python two coupled second order ODEs Runge Kutta 4th order

Question

I am a beginner in python. I have a problem with 2 ODEs that are second order and they are coupled. I want to solve it with Runge Kutta 4th order. I have made 2 matrices.[A] and [B] that V' = A*C + B . But I couldn't get the answer. Can anyone help me?

This is equations

Equations

this is my code:

import numpy as np
from math import sqrt

from sympy import *
import matplotlib.pyplot as plt

x,v1, v2, v3, v4, dv1, dv2, dv3, dv4, t = symbols('x v1 v2 v3 v4 dv1 dv2 dv3 dv4 t')
# Parameters
mr = 
m =
c = 
k = 
F =  

# Equations:


A = Matrix([[0, 0, 1, 0], [0, 0, 0, 1], [-k / mr, k / mr, -c / mr, c / mr], [k / m, -k / m, c / m, -c / m]])

B = Matrix([0, 0, 0, -F])

dv1 = diff(v1, t)
dv2 = diff(v2, t)
dv3 = diff(v3, t)
dv4 = diff(v4, t)


C = Matrix([v1,v2,v3,v4])


def V(v1,v2,v3,v4,t):
    V = A*C+B

    return V

def rk4(f, v1, v2, v3 ,v4 , x0 , x1 , n):
    k11 = [0] * (n + 1)
    k21 = [0] * (n + 1)
    k31 = [0] * (n + 1)
    k41 = [0] * (n + 1)
    k12 = [0] * (n + 1)
    k22 = [0] * (n + 1)
    k32 = [0] * (n + 1)
    k42 = [0] * (n + 1)
    k13 = [0] * (n + 1)
    k23 = [0] * (n + 1)
    k33 = [0] * (n + 1)
    k43 = [0] * (n + 1)
    k14 = [0] * (n + 1)
    k24 = [0] * (n + 1)
    k34 = [0] * (n + 1)
    k44 = [0] * (n + 1)


    h = (x1 - x0) / n
    for i in range(1, n + 1):

        k11 = h * f(v1, v2, v3, v4, t)
        k21 = h * f(v1, v2, v3, v4, t)
        k31 = h * f(v1, v2, v3, v4, t)
        k41 = h * f(v1, v2, v3, v4, t)

        k12 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k11)
        k22 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k21)
        k32 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k31)
        k42 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k41)

        k13 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k12)
        k23 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k22)
        k33 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k32)
        k43 = h * f(v1 + 0.5 * h, v2 + 0.5 * h, v3 + 0.5 * h, v4 + 0.5 * h, t + 0.5 * k42)

        k14 = h * f(v1 + h, v2 + h, v3 + h, v4 + h, t + k13)
        k24 = h * f(v1 + h, v2 + h, v3 + h, v4 + h, t + k23)
        k34 = h * f(v1 + h, v2 + h, v3 + h, v4 + h, t + k33)
        k44 = h * f(v1 + h, v2 + h, v3 + h, v4 + h, t + k43)

        t[i] = t + i * h
        jv1[i] = v1 + (k11 + k12 + k12 + k13 + k13 + k14) / 6
        jv2[i] = v2 + (k21 + k22 + k22 + k23 + k23 + k24) / 6
        jv3[i] = v3 + (k31 + k32 + k32 + k33 + k33 + k34) / 6
        jv4[i] = v4 + (k41 + k42 + k42 + k43 + k43 + k44) / 6
    return jv1, jv2 , jv3, jv4


jv1, jv2 , jv3 , jv4  = rk4(V, 0, 0, 1, 1 , 0 , 1 , 10)

plt.plot(jv1, jv2 , jv3 , jv4 , t)
plt.grid('on')
plt.show()

Can you clarify what your question is? You say you "couldn't get the answer", but you don't actually say anywhere what does happen. Are you getting an exception? The wrong output? No output? What do you expect the output to be like? — Blckknght, Sep 16 '18 at 17:54
Please add your reasoning for using sympy in a numerical code. Check again how RK4 works. Decide if you want a vector based approach or a coordinate wise approach. At the moment it is a mix that just can not work. — Lutz Lehmann, Sep 16 '18 at 18:27
I got error in code , it could not run. get like 10errors after using rk4 function.I need to get 4 variables in list — faraz parvaresh, Sep 16 '18 at 19:47
I used sympy because I have 4 coupled variables it cann't solved dependently so I have a matrix for solving it. — faraz parvaresh, Sep 16 '18 at 19:49

Lutz Lehmann · Accepted Answer · 2020-10-03T20:33:11.930

Cutting away unused and contradicting variables and statements a running reduction and correction of your code is

import numpy as np
from math import sqrt

import matplotlib.pyplot as plt

# Parameters
mr = 2
m = 1
c = 5
k = 36
F = 3

# Equations:
   
def V(u,t):
    x1, x2, v1, v2 = u
    dx = x2-x1;
    dv = v2-v1;

    return np.array([ v1, v2, (c*dv+k*dx)/mr, -(F+c*dv+k*dx)/m ])

def rk4(f, u0, t0, tf , n):
    t = np.linspace(t0, tf, n+1)
    u = np.array((n+1)*[u0])
    h = t[1]-t[0]
    for i in range(n):
        k1 = h * f(u[i], t[i])    
        k2 = h * f(u[i] + 0.5 * k1, t[i] + 0.5*h)
        k3 = h * f(u[i] + 0.5 * k2, t[i] + 0.5*h)
        k4 = h * f(u[i] + k3, t[i] + h)
        u[i+1] = u[i] + (k1 + 2*(k2 + k3) + k4) / 6
    return u, t


u, t  = rk4(V, np.array([0., 0., 1., 1.]) , 0. , 1. , 10)
x1, x2, v1, v2 = u.T
plt.plot(t, x1, t, x2)
plt.grid('on')
plt.show()

If that is in any way a correct result is up to you to determine.

python two coupled second order ODEs Runge Kutta 4th order

1 Answers1

Linked