Ода не работает для уравнения Стюарта Ландау?

Question

Я не уверен, что проблема в этом куске кода. Я пытаюсь решить следующее ODE. Я решил это с 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

Считая первый компонент более правильным и предполагая, что система является реальной и мнимой частью сложного ODE, я прочитал для этого 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. численное решение точно соответствует ожидаемому от точного решения.

Answer 2

Вы вызываете .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.