python - 使用 RK4 显式求解 ODE

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)

如何求解作为时间函数的 theta 并从 t=0 到 t=40 绘制。系统从静止开始theta = 0 and d(theta)/dt = 0。

Answer 1

似乎您要模拟具有振荡外力的强制物理钟摆。

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

正如您已经正确认识到的那样，这需要转换为 1 阶系统，然后打包为向量值 ODE func

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 积分器的标准格式。