尝试使用四阶龙格-库塔法求解非线性摆二阶微分方程，未得到预期结果

Question

我是 Python 新手，所以我只将过程限制为 matplotlib，而不是进入 NumPy。

我按照 Abhijit Kar Gupta 的书《Python 中的科学计算》编写了 Python 脚本，使用四阶 Runge-Kutta 方法求解二阶微分方程。

书中，作者给出的例子很简单。但对于我的问题，非线性摆（其中小角度近似 sin(x)~x 不准确），我尝试按照我的预期修改函数。

我正在尝试绘制角度 (x) 与时间 (t) 的变化。但我认为我没有得到正确的输出。以 x=170 度角度释放的摆锤的性质如下：https://en.wikipedia.org/wiki/Pendulum_(mechanics)#/media/File:Pendulum_170deg.gif。\

虽然我发现 x = 160 度 = 8.0*pi/9.0 弧度，但该图仍然应该看起来有点像正弦振荡，周期不均匀。但我得到的只是像突然停止的钟摆一样。

请告诉我哪里做错了。另外，如果我不想查看绘图，而是想要获得类似列表的输出（例如包含时间和角度的两列的表格），我该怎么办？另外，应该做什么才能获得特定时间（即任何外部输入）的角度？

这是我写的代码：

Python 3.12.0 (tags/v3.12.0:0fb18b0, Oct  2 2023, 13:03:39) [MSC v.1935 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license()" for more information.

# Solving a non-linear pendulum equation of motion using RK4 method:
# Equation of Motion: d^x/dt^2 + k.sin(x)=0, where 'x' is angle, k=g/l. We take g=10m/s^2, l=10m, so k=1(s^{-2}).
# We modify the single 2nd order differential equation as a system of linear equations, consisting of two 1st order differential equation: (1) dx/dt=y=f1(x,y,t), (2) dy/dt=-k.sin(x)=f2(x,y,t).

from math import sin, pi

# Defining two functions-
def f1(t,x,y): return y

def f2(t,x,y): return -k*sin(x)

k=1.0                      #parameter
t,x,y=0,8.0*pi/9.0,0             #initial values (t: second, x: radian, y: (initial angular velocity) radian/second)
h=0.01                   #increment in t
T,X,Y=[t],[x],[y]          #lists to store data

# Loop starts with RK4 steps:
for i in range(200):
    a1=h*f1(t,x,y)                      
    b1=h*f2(t,x,y)                      
    a2=h*f1(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    b2=h*f2(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    a3=h*f1(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    b3=h*f2(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    a4=h*f1(t+h,x+a3,y+b3)              
    b4=h*f2(t+h,x+a3,y+b3)              
    x=x+(1/6)*(a1+2*a2+2*a3+a4)         # apprx. value of x1 for t+h
    y=y+(1/6)*(b1+2*b2+2*b3+b4)         # apprx. value of y1 for t+h
    t=t+h                               # current value of independent variable t
    T.append(t)
    X.append(x)
    Y.append(y)

    
# for plotting
import matplotlib.pyplot as plt
plt.plot(T,X,lw=3)
[<matplotlib.lines.Line2D object at 0x00000190D0CD5910>]
plt.xlim(0,50)    # xlim is just time range
(0.0, 50.0)
plt.ylim(-pi,pi)  # ylim is the range of angle (in radians)
(-3.141592653589793, 3.141592653589793)
plt.grid()
plt.show()

我得到的情节：

如果我放大蓝线：

这与一般的钟摆相差甚远。为什么解决方案突然停止？

Answer 1

正如评论中所指出的，问题是您只选择了 20 个样本进行绘图。您可以将样本数量参数化为时间的函数，从而得到一个图表，其中 t 的范围从 0 到 50，如您所愿：

# Solving a non-linear pendulum equation of motion using RK4 method:
# Equation of Motion: d^x/dt^2 + k.sin(x)=0, where 'x' is angle, k=g/l. We take g=10m/s^2, l=10m, so k=1(s^{-2}).
# We modify the single 2nd order differential equation as a system of linear equations, consisting of two 1st order differential equation: (1) dx/dt=y=f1(x,y,t), (2) dy/dt=-k.sin(x)=f2(x,y,t).

from math import sin, pi

# Defining two functions-
def f1(t,x,y): return y

def f2(t,x,y): return -k*sin(x)

k=1.0                      #parameter
t,x,y=0,8.0*pi/9.0,0             #initial values (t: second, x: radian, y: (initial angular velocity) radian/second)
h=0.01                   #increment in t
T,X,Y=[t],[x],[y]          #lists to store data
tmax=50                  #max value of t (duration of the plot in units of time)
n=int(tmax/h)            #number of samples

# Loop starts with RK4 steps:
for i in range(n):
    a1=h*f1(t,x,y)                      
    b1=h*f2(t,x,y)                      
    a2=h*f1(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    b2=h*f2(t+0.5*h,x+0.5*a1,y+0.5*b1)  
    a3=h*f1(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    b3=h*f2(t+0.5*h,x+0.5*a2,y+0.5*b2)  
    a4=h*f1(t+h,x+a3,y+b3)              
    b4=h*f2(t+h,x+a3,y+b3)              
    x=x+(1/6)*(a1+2*a2+2*a3+a4)         # apprx. value of x1 for t+h
    y=y+(1/6)*(b1+2*b2+2*b3+b4)         # apprx. value of y1 for t+h
    t=t+h                               # current value of independent variable t
    T.append(t)
    X.append(x)
    Y.append(y)

    
# for plotting
import matplotlib.pyplot as plt
plt.plot(T,X,lw=3)
plt.xlim(0,tmax)    # xlim is just time range
plt.ylim(-pi,pi)  # ylim is the range of angle (in radians)
plt.grid()
plt.show()