# System of second order ODEs Runge Kutta 4th order

I have two masses, and two springs which are atattched to a wall on one end. For context I attached the system of equations. How can I add a command to see my analytic solution for my system of equations? I want to be able to compare the numeric vs analytic. I think it would be something like (t, answer, t, answer). Thank you.

import numpy as np
import matplotlib.pyplot as plt

def f(x,t):
k1=20
k2=30
m1=3
m2=5
return np.array([x, (-k1*x-k2*(x - x))/m1, x, (-k2*(x-x)/m2)])

h=.01

t=np.arange(0,15+h,h)

y = np.zeros((len(t), 4))
y[0, :] = [1, 0, 0, 0]

for i in range(0,len(t)-1):
k1 = h * f( y[i,:], t[i] )
k2 = h * f( y[i,:] + 0.5 * k1, t[i] + 0.5 * h )
k3 = h * f( y[i,:] + 0.5 * k2, t[i] + 0.5 * h )
k4 = h * f( y[i,:] + k3, t[i+1] )
y[i+1,:] = y[i,:] + ( k1 + 2.0 * ( k2 + k3 ) + k4 ) / 6.0

plt.plot(t,y[:,0],t,y[:,2])

plt.gca().legend(('x_1','x_2'))
plt.show()
$$$$

• The symbolic solutions of sympy are not yet reliable, there was progress in the recent versions, but also problems reported with linear systems. You could compare your RK4 solution against the likewise numerical scipy.integrate.odeint solution with sufficiently low error tolerances. – Lutz Lehmann Nov 23 '20 at 14:55
• Agreed. It's almost never a good idea to reimplement a time integrator yourself when there already exists much better packages. I would rather suggest you use the solve_ivp module instead of odeint. It is more modern and is IMO clearer in terms of inputs and outputs. As for the analytical solution, your system is linear, so you should be able to find it by trasnforming your equations to a 1st-order system $\dot{x} = Ax$, computing the eigenvectors and eigenvalues of $A$, and solving analytically in the space of eigenvector before going back to the original space. – Laurent90 Nov 23 '20 at 15:02

Here is a quick rewrite of your code, with the addition of the analytical solution using a matrix exponential:

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg
import scipy.integrate

k1=20
k2=30
m1=3
m2=5
A = np.array([[0,1,0,0],
[-(k1+k2)/m1, 0, k2/m1, 0],
[0,0,0,1],
[k2/m2,0,-k2/m2,0]])
def f(t,x): # /!\ changed the order of the arguments for solve_ivp
return A.dot(x)

x0 = np.array([1, 0, 0, 0])
tf = 2.

# compute numerical solution with adaptive time stepping
sol = scipy.integrate.solve_ivp(fun=f, y0=x0, t_span=(0,tf), method='DOP853', atol=1e-12, rtol=1e-12)

# compute analytical solution for the linear system
# dx/dt = Ax  --> x(t) = exp(t*A) * x(0)
dt_exp = 0.05
t_exp = np.arange(0,tf,dt_exp) # time points
exp_dtA = scipy.linalg.expm(dt_exp * A ) # used to compute the solution during one time step
sol_exp=[x0]
for t in t_exp[1:]:
sol_exp.append( exp_dtA.dot( sol_exp[-1] ) )
sol_exp = np.array(sol_exp).T

plt.figure()
plt.plot(sol.t, sol.y[0,:], label=r'$$x_1$$', color='r', linestyle='', marker='.')
plt.plot(sol.t, sol.y[2,:], label=r'$$x_2$$', color='b', linestyle='', marker='.')
plt.plot(t_exp, sol_exp[0,:], label=r'$$x_{1,th}$$', color='r', linestyle='-', marker=None)
plt.plot(t_exp, sol_exp[2,:], label=r'$$x_{2,th}$$', color='b', linestyle='-', marker=None)
plt.ylabel('position (m)')
plt.xlabel('t (s)')
plt.grid()
plt.legend()
`

Here is the result: • Thank you for your help, I found this very useful. – J Wright Nov 23 '20 at 16:32
• Has this solved your problem ? Or do you need further help ? – Laurent90 Nov 23 '20 at 17:57
• This was very helpful for me, so it answers it. Thank you very much. – J Wright Nov 23 '20 at 18:35
• You're welcome ! Can you accept the answer so that its status is updated ? – Laurent90 Nov 23 '20 at 18:45