Numerical solutions to differential equation

Question

I have been given the following question:

---

The differential equation describing the angular position $\theta$ of a mechanical arm: $$\ddot{\theta}=\frac{a(b-\theta)-\theta\dot{\theta}^2}{1-\theta^2}$$ where $a=100ms^{-2}$, b=15. If $\theta(0)=2\pi, \dot{\theta}=0$. Compute $\theta$ & $\dot{\theta}$ when t=0.5s with tolerance of $10^{-3}$.

---

I am required to create a python program to solve this. I have no problem with making the program. However, normally these problems require me to split the differential equation into a system of first order differential equations (by introducing a new variable).

The reason for posting is to see if anyone can help me to split this into a system of DE's, as I have no idea where to start.

Answer 1

(The task in the question is ambiguous and gives a not solvable problem in the apparent interpretations. However, the original task that can be found in the exercise 14 on page 271 of Numerical Methods in Engineering with MATLAB® By Jaan Kiusalaas, is well solvable.)

The usual method to transform a higher order ODE into a first order system via the derivatives array $\dot \theta = \omega$, $\dot\omega=rhs(\theta,\omega)$ can be implemented in a python function as

a, b = 100.0, 15.0

def ODEsystem(t,u):
    theta, omega = u
    return [ omega, (a*(b-theta)-theta*omega**2)/(1+theta**2) ]

theta0 = 2*pi; omega0 = 0
u0 = [ theta0, omega0 ]

---

This then can be fed to an ODE solver and evaluate for instance as in the following lines where Richardson extrapolation is used to estimate the error of the second order method.

steps = [2**k for k in range(10)]
sol = [ RK2(ODEsystem, u0, [0.0,0.5], subdiv = N)[-1,0] for N in steps]
for k in range(1,len(steps)):
    print steps[k], sol[k],(sol[k-1]-sol[k])/3, (4*sol[k]-sol[k-1])/3 

with the result

steps    theta(0.5)     error estimate      extrapolation

  2    8.62612817225    0.116288833908       8.50983933834
  4    8.44007543611    0.0620175787123      8.3780578574
  8    8.39252457314    0.0158502876557      8.37667428549
 16    8.38067290584    0.00395055576866     8.37672235007
 32    8.37773561225    0.000979097860883    8.37675651439
 64    8.37700612603    0.000243162075971    8.37676296395
128    8.37682447081    6.05517382795e-05    8.37676391907
256    8.37677915382    1.5105664442e-05     8.37676404815
512    8.37676783713    3.7722280369e-06     8.37676406491

---

For demonstration purposes a second-order fixed-step Runge-Kutta method was used, the explicit midpoint or improved Euler method.

def RK2(odefunc, u, times, subdiv = 1):
    f = lambda u,t: np.array(odefunc(t,u))
    u = np.array(u); times = np.array(times)
    uout = np.zeros((len(times),)+u.shape)
    uout[0] = u;
    for k in range(len(times)-1):
        t = times[k]
        h = (times[k+1]-times[k])/subdiv
        for j in range(subdiv):
            k1 = f(u,t)*h
            k2 = f(u+0.5*k1, t+0.5*h)*h
            u, t = u+k2, t+h
        uout[k+1]=u
    return uout

Answer 2

The new variable is usually the derivative of the variable you have, so define $x=\dot \theta$ and your equation becomes $$\dot x=\frac{a(b-\theta)-\theta x^2}{1-\theta^2}\\ \dot {\theta}=x$$ which is first order in both variables.