python两个耦合二阶ODE Runge Kutta第四顺序

Question

我是python的初学者。我有2个二阶ODE的问题，它们是耦合的。我想用Runge Kutta第四顺序解决它。我制作了2个矩阵。[A]和[B]表示V'= A * C + B.但我无法得到答案。谁能帮我？

这是方程式

方程

这是我的代码：

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

Answer 1

切除未使用和矛盾的变量和语句，减少和纠正代码

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] + 0.5*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()

如果这是以任何方式，您需要确定正确的结果。