求解非线性一阶微分方程并在图中打断

Question

我正在尝试使用 python 中的四阶 runge-kutta 方法求解椭圆微分方程。

执行后，我得到了应该获得的实际情节的一小部分，并伴随着错误提示：

"RuntimeWarning: invalid value encountered in double_scalars"

import numpy as np
import matplotlib.pyplot as plt

#Define constants
g=9.8
L=1.04

#Define the differential Function
def fun(y,x):
    return-(2*(g/L)*(np.cos(y)-np.cos(np.pi/6)))**(1/2)

#Define variable arrays
x=np.zeros(1000)
y=np.zeros(1000)
y[0]=np.pi/6
dx=0.5

#Runge-Kutta Method
for i in range(len(y)-1):
    k1=fun(x[i],y[i])
    k2=fun(x[i]+dx/2, y[i]+dx*k1/2)
    k3=fun(x[i]+dx/2, y[i]+dx*k2/2)
    k4=fun(x[i]+dx, y[i]+dx*k3)
    
    
    y[i+1]=y[i]+dx/6*(k1+2*k2+2*k3+k4)
    x[i+1]=x[i]+dx
    
#print(y)
#print(x)
plt.plot(x,y)
plt.xlabel('Time')
plt.ylabel('Theta')
plt.grid()

我得到的图表是这样的，

我的问题是为什么我会收到错误消息？感谢您的帮助！

Answer 1

导致此行为的几点。首先，您切换了 ODE 函数中参数的顺序，可能是为了使其与 odeint 兼容。使用 tfirst=True 可选参数来避免这种情况，并始终首先使用自变量。

错误的实际来源是术语

(np.cos(y)-np.cos(np.pi/6)))**(1/2)

请记住，在您的版本中 y 的值为 x[i]，因此在某些时候根下的表达式变为负数。

如果你纠正了第一点，你可能仍然会遇到第二个错误，因为精确解是抛物线式地向不动点移动，所以RK4的阶段很可能会超调。可以通过提供足够安全的平方根函数来解决这个问题，

def softroot(x): return x/max(1e-12,abs(x))**0.5

#Define the differential Function
def fun(x,y):
    return -(2*(g/L)*softroot(np.cos(y)-np.cos(np.pi/6)))

#Define variable arrays
dx=0.01
x=np.arange(0,1,dx)
y=np.zeros(x.shape)
y[0]=np.pi/6
...

产生一个情节

因为解决方案已经在固定点开始。将初始点稍微向下移动到 y[0]=np.pi/6-1e-8 会跳转到下面的固定点。

求解非线性一阶微分方程并在图中打断

solving a Non- linear first order differential equation and getting a break in the plot

python

differential-equations