I have two masses, and two springs which are atattched to a wall on one end. For context I attached the system of equations. enter image description here

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):
    return np.array([x[1], (-k1*x[0]-k2*(x[0] - x[2]))/m1, x[3], (-k2*(x[2]-x[0])/m2)])



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


  • 1
    $\begingroup$ 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. $\endgroup$ Nov 23, 2020 at 14:55
  • $\begingroup$ 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. $\endgroup$
    – Laurent90
    Nov 23, 2020 at 15:02

1 Answer 1


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

A = np.array([[0,1,0,0],
              [-(k1+k2)/m1, 0, k2/m1, 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
for t in t_exp[1:]:
    sol_exp.append( exp_dtA.dot( sol_exp[-1] ) )
sol_exp = np.array(sol_exp).T

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)')

Here is the result:

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.