explicitly solving ODEs using RK4

Question

我的两个一阶微分如下

y1' = -sin(y0) + (gamma)*cos(y0)sin(beta * x)

和

y0' = y1

其中

(theta)'' = y, (theta)' = y1, theta = y0

我最初的方程式是

(((d^2)*theta)/dt^2)=-sin(theta)+(gamma)cos(theta)sin(Bx)

如何求解θ作为时间的函数，并绘制从t=0到t=40的曲线。系统以theta = 0 and d(theta)/dt = 0静止启动。

Answer 1

似乎要模拟带有振荡外力的强迫物理摆。

theta''+sin(theta) = gamma * cos(theta)*sin(beta*t)

正如您正确认识到的，这需要转换为1阶系统，然后打包为向量值ODE函数

def odefunc(t,y):
    y0, y1 = y
    return y1, -sin(y0)+gamma*cos(y0)*sin(beta*t)

如果常数是全局变量，或者作为参数

def odefunc(t,y,params):
    y0, y1 = y
    beta, gamma = params
    return y1, -sin(y0)+gamma*cos(y0)*sin(beta*t)

然后必须使用lambda表达式将参数减少为ODE积分器的标准格式。