CODEWITHSUNDEEP
python
J Wright
J Wright

Reputation: 21

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. Below is the python code for the 4th order Runge-Kutta that evaluates the following system of two 2nd order ODE:

enter image description here

I need help fixing it. When I run my code it plots a flat line, but I know it should oscillate because I want to plot a position versus time graph. Something is probably wrong in my last few lines for the plot commands, but I am not sure what. Thank you for any help.

import numpy as np
import matplotlib.pyplot as plt


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

h=.01

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

y=np.zeros((len(t),4))

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

Upvotes: 1

Views: 1143

Answers (1)

Warren Weckesser
Warren Weckesser

Reputation: 114976

You'll need a nonzero initial condition, e.g,

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

and you'll also have to fix the derivatives calculated in f:

    return np.array([x[1], (-k1*x[0]-k2*(x[0] - x[2]))/m1, x[3], (-k2*(x[2]-x[0])/m2)])

With those changes, the plot is

plot

P.S. You might be interested in this: Coupled spring-mass system

Upvotes: 1

Related Questions