二阶 ODE 系统 Runge Kutta 四阶答案

【问题标题】：System of second order ODEs Runge Kutta 4th order二阶 ODE 系统 Runge Kutta 四阶
【发布时间】：2020-11-23 14:49:07
【问题描述】：

我有两个质量，两个弹簧，一端固定在墙上。对于上下文，我附上了方程组。下面是四阶 Runge-Kutta 的 Python 代码，它评估以下两个二阶 ODE 的系统：

我需要帮助修复它。当我运行我的代码时，它会绘制一条平线，但我知道它应该会振荡，因为我想绘制位置与时间的关系图。我的最后几行 plot 命令可能有问题，但我不确定是什么。感谢您的帮助。

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

【问题讨论】：

标签： python

【解决方案1】：

你需要一个非零的初始条件，例如，

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

您还必须修正在f 中计算的导数：

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

随着这些变化，情节是

附：您可能对此感兴趣：Coupled spring-mass system

【讨论】：