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

Question

我尝试使用odeint 解决达芬方程：

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();

精确匹配