用 python 求解非线性非齐次微分方程（Duffing 振荡器）

Question

我尝试使用 odeint:

来求解 Duffing 方程

def func(z, t):
    q, p = z
    return [p, F*np.cos(OMEGA * t) + Fm*np.cos(3*OMEGA * t) + 2*gamma*omega*p - omega**2 * q - beta * q**3]
OMEGA = 1.4
T = 1 / OMEGA
F = 0.2
gamma = 0.1
omega = -1.0
beta = 0.0
Fm = 0

z0 = [1, 0]

psi = 1 / ((omega**2 - OMEGA**2)**2 + 4*gamma**2*OMEGA**2)
wf = np.sqrt(omega**2 - gamma**2)

t = np.linspace(0, 100, 1000)
sol1 = odeintw(func, z0, t, atol=1e-13, rtol=1e-13, mxstep=1000)

当F = gamma = beta = 0时，我们有一个由两个线性齐次方程组成的系统。很简单！

但是当F not equal 0 系统变得不均匀。问题是数值解和解析解不重合：

^图一

数值解没有考虑方程的非齐次解

我还没弄清楚是否可以在这里使用 solve_bvp 功能。你能帮帮我吗？

Answer 1

插入常数，等式变为

x'' + 2*c*x' + x = F*cos(W*t)

一般解的形式是

x(t)=A*cos(W*t)+B*sin(W*t)+exp(-c*t)*(C*cos(w*t)+D*sin(w*t))

w^2=1-c^2

对于得到的特定解决方案

  -W^2*(A*cos(W*t)+B*sin(W*t))
+2*c*W*(B*cos(W*t)-A*sin(W*t))
+     (A*cos(W*t)+B*sin(W*t))
=F*cos(W*t)

 (1-W^2)*A + 2*c*W*B = F
-2*c*W*A + (1-W^2)*B = 0

初始条件需要

A+C = 1
W*B-c*C+w*D=0

在python中这样编码

F=0.2; c=0.1; W=1.4
w=(1-c*c)**0.5

A,B = np.linalg.solve([[1-W*W, 2*c*W],[-2*c*W,1-W*W]], [F,0])
C = 1-A; D = (c*C-W*B)/w

print(f"w={w}, A={A}, B={B}, C={C}, D={D}")

与输出

w=0.99498743710662, A=-0.192, B=0.056, C=1.192, D=0.04100554286257586

继续

t = np.linspace(0, 100, 1000)
u=odeint(lambda u,t:[u[1], F*np.cos(W*t)-2*c*u[1]-u[0]], [1,0],t, atol=1e-13, rtol=1e-13)
plt.plot(t,u[:,0],label="odeint", lw=3); 
plt.plot(t,A*np.cos(W*t)+B*np.sin(W*t)+np.exp(-c*t)*(C*np.cos(w*t)+D*np.sin(w*t)), label="exact"); 
plt.legend(); plt.grid(); plt.show();

结果完全吻合

用 python 求解非线性非齐次微分方程（Duffing 振荡器）

Solve non-linear non homogeneous differential equation with python (Duffing oscillator)

python

numerical-methods

differential-equations

numerical-integration