在这个网站和我的参考书上搜索后,我发现我不知道为什么我的代码不起作用。
正如我的教授在课堂上向我们展示的那样,我为质量弹簧系统(具有摊销)做了一个四阶 Runge-Kutta 实现。但是,如您所见,生成的图形非常奇怪。
我最终写的代码是这样的:
#! /usr/bin/env python3
#-*- coding: utf-8 -*-
def f(data, t, x1, v1):
from math import cos
F = data["F"]
c = data["c"]
k = data["k"]
m = data["m"]
omega = data["omega"]
return( [v1, (F*cos(omega*t) - c*v1 - k*x1)/m] )
def run(data = {}):
xi, vi, ti = [data["u1"]], [data["v1"]], [data["t_ini"]]
dt = data["dt"]
while ti[-1] <= data["t_end"]:
xn = xi[-1]
vn = vi[-1]
tn = ti[-1]
K1 = f(data, t = tn, x1 = xn, v1 = vn)
K1 = [dt*K1[i] for i in range(len(K1))]
K2 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K1[0], v1 = vn + 0.5*K1[1])
K2 = [dt*K2[i] for i in range(len(K2))]
K3 = f(data, t = tn + 0.5*dt, x1 = xn + 0.5*K2[0], v1 = vn + 0.5*K2[1])
K3 = [dt*K3[i] for i in range(len(K3))]
K4 = f(data, t = tn + dt, x1 = xn + K3[0], v1 = vn + K3[1])
K4 = [dt*K4[i] for i in range(len(K4))]
xn = xn + (K1[0] + 2*K2[0] + 2*K3[0] + K4[0])/6
vn = xn + (K1[1] + 2*K2[1] + 2*K3[1] + K4[1])/6
ti.append(tn+dt)
xi.append(xn)
vi.append(vn)
return(ti, xi, vi)
这是由 main.py 文件导入的,该文件仅包含程序的 GUI 和绘图部分,并且该函数是在课堂上推导出来的,所以我相信错误出在 Runge-Kutta 本身。(可能是我搞砸了一些愚蠢的事情。)
我尝试在“xn”和“vn”中切换 K,在 f() 中强制“F”和“c”值,重写所有内容并手动编写每个 K 的每个元素(如 K11、K12、K21 ,等等),但它只给出指数结果。此外,将 f() 的返回切换为 numpy 数组并没有任何帮助。
我发现了一些关于RK4方法的问题,但是我无法解决这个问题,也不明白哪里出了问题。我对这个方法有一些了解,但这实际上是我第一次实现它,所以请大家帮忙。
如果重要的话,我正在为 python3 使用 Anaconda 发行版。