Runge–Kutta methods for ODE integration in Python with additional constraints

Question

I have a question on solving second order differential equations using RK4, considering additional constraints on the first derivative. I am doing the example shown here with some modifications

θ′′(t) + b θ′(t) + c sin(θ(t)) = 0

The additional constraint is:

when θ_i θ_(i+1)<0, then θ′_(i+1)=0.9θ′_i,

where i is the steps of t, i+1 is one step after i. In the real world, it says when the direction of displacement reverses, its velocity is reduced to 90%.

Vectorially if y(t) = (θ(t), ω(t)), then the equation is ẏ = f(t,y), where f(t,y) = (y₂(t), −by₂(t) − cos(y₁(t))).

In this problem, how should I modify the code if I want to add constraints on ω or θ′(t) (which are the same thing)? Here is my code which didn't work. The additional condition makes θ′ non-continuous. The following "homemade" solution cannot update θ′ properly. I am new to Python and this is my first StackOverflow post. Any guidance is much appreciated.

def rungekutta4(f, y0, t, args=()):
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    for i in range(n - 1):
        h = t[i+1] - t[i]
        if y[i][0]*y[i+1][0]<0:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)*0.9
        else:
            k1 = f(y[i], t[i], *args)
            k2 = f(y[i] + k1 * h / 2., t[i] + h / 2., *args)
            k3 = f(y[i] + k2 * h / 2., t[i] + h / 2., *args)
            k4 = f(y[i] + k3 * h, t[i] + h, *args)
            y[i+1] = y[i] + (h / 6.) * (k1 + 2*k2 + 2*k3 + k4)        
    return y

Answer 1

Unless I completely misunderstand you, what you want is a simple case distinction in f: Mathematically, you have f₂(t,y) = −by₂(t) − cos(y₁(t)) if θ_i²−1 = y₁(t)²−1 < 0 and 0.9·(y₂−1) otherwise. This is all still only a dependency of f on y, which can simply implemented programming-wise.

For example, you could define:

def f(y):
    θ = y[0]
    ω = y[1]
    return [
        θ,
        -b*ω-cos(θ) if θ**2<1 else 0.9*(ω-1)
    ]

While this may work as it is, you may face problems due to f not being continuous (or having a continuous derivative), in particular if you want to use more advanced integrators with step-size control instead of your homebrewed one. To avoid this, replace the if–else construct with a (sharp) sigmoid. For more details on this, see this answer of mine.

Answer 2

In the current formulation and taking the idea that each time the pendulum passes the vertical position its velocity is reduced by 10%, this can be approximately arranged as

    for i in range(n - 1):
        h = t[i+1] - t[i]
        y[i+1] = RK4step(f,t[i],y[i],h, args)
        if y[i+1,0]*y[i,0] < 0: y[i+1,1] *= 0.9
    return y

that is, first compute the new value and then apply the condition. The time step should be small enough that the angle only changes a few degrees. For larger time steps you would have to split the step containing the zero crossing, using some root finding method like the secant method to find a more accurate time of the root.