常微分方程不适用于斯图亚特·兰道方程吗？

Question

我不确定这段代码中有什么问题。我正在试着解决下面这首颂歌。我用ODEint解决了这个问题，由于某种原因，它只给出了一个简单的答案，这肯定是不正确的。我已经纠正了雅可比矩阵中的错误，并更改了参数，但都不起作用。我用FreeFem++用同样的参数解决了这段代码的一个版本，所以这是令人惊讶的。

from scipy.integrate import ode
import numpy as np

A0= [10**-4, 0]
t0 = 0
B0 = 10**-4 + 0j


coeff = np.complex(0.145,-0.088)
coeffr = np.real(coeff)
coeffi = np.imag(coeff)


def f(t,A):
    return [(A[0]**2 + A[1]**2)*(- coeffr*A[0] + coeffi*A[1]), (A[0]**2 + A[1]**2)*(-coeffr*A[0] - coeffi*A[1])]

def jac(t,A):
    return [[ (A[0]**2 + A[1]**2)*(- coeffr) + (2*A[0])*(- coeffr*A[0] + coeffi*A[1]), (A[0]**2 + A[1]**2)*(coeffi) + (2*A[1])*(- coeffr*A[0] + coeffi*A[1])],
            [ (A[0]**2 + A[1]**2)*(-coeffr) + (2*A[0])*(-coeffr*A[1] - coeffi*A[0]) , (A[0]**2 + A[1]**2)*(- coeffi) +  (2*A[1])*(-coeffr*A[0] - coeffi*A[1])]]

def fcomp(t,B):
    return -coeff*abs(B)**2*B

s = ode(fcomp).set_integrator('zvode', method='bdf')
s.set_initial_value(B0, t0)
t1 = 100
dt = 0.01

import numpy as np

solution = np.empty((0,100))
time = np.empty((0,100))

while s.successful() and s.t < t1:
    sol = s.integrate(s.t+dt)[0]
    solution = np.append(solution,sol)
    tim = s.t+dt
    time = np.append(time,tim)
    print(s.t+dt, s.integrate(s.t+dt))


import matplotlib.pyplot as plt


plt.plot(time,np.abs(solution))
plt.show()

致以最亲切的问候，

凯瑟琳

Answer 1

您调用了.set_f_params(2.0).set_jac_params(2.0)，但是f和jac都没有额外的参数。变化

r.set_initial_value(A0, t0).set_f_params(2.0).set_jac_params(2.0)

只是为了

r.set_initial_value(A0, t0)

您将遇到的另一个问题是f和jac中的表达式缺少一些乘法符号。例如，在f中，术语

(A[0]**2 + A[1]**2)(- coeffr*A[0] + coeffi*A[1])

应该是

(A[0]**2 + A[1]**2) * (- coeffr*A[0] + coeffi*A[1])

同样的问题发生在第二个学期，在jac中的几个地方。

Answer 2

假设第一个部分可能更正确，并假设系统是复杂颂歌的真实和虚构部分，我读取了那个z'=-c*z*|z|^2。

请注意，实部和虚部是

A2*(-coeffr*A[0]+coeffi*A[1]),  A2*(-coeffr*A[1]-coeffi*A[0])

请注意虚数部分中的切换索引。

使用新的接口solve_ivp，整个任务可以写得更短

z0 = 1e-4+0j
coeff = 0.145-0.088j
t0, tf, dt = 0, 100, 0.01

sol = solve_ivp(lambda t,z: -coeff*z*abs(z)**2, [t0,tf], [z0], t_eval=np.arange(t0,tf,dt))

plt.plot(sol.t, abs(sol.y[0]))

这确实给出了几乎恒定的解决方案，其中1e-11的范围变化在误差容限之内。

这是否合理呢？绝对值的方程式为

d/dt abs(z)^2 = 2*Re(z.conj*dz/dt) = -2*coeffr*abs(z)^4

它有一个解决方案

abs(z)^2 = abs(z(0))^2/(1+2*coeffr*abs(z(0))^2*t)

实际上，在t=100，分母是1+2*0.145*1e-8*100=1+2.9e-7，它非常接近于1，数值解与精确解的预期完全一致。