как использовать solve_ivp решить дифференциальные уравнения в частных производных спектральным методом? - PullRequest
Новая
       73

как использовать solve_ivp решить дифференциальные уравнения в частных производных спектральным методом?

0 голосов
/ 01 апреля 2020

Я хочу использовать спектральный метод для решения уравнений в частных производных. Уравнения, подобные этому, формула , исходное условие: u (t = 0, x) = (a ^ 2) * sech (x), u'_t (t = 0) = 0.

Для решения этой проблемы я использую python спектральным методом. Ниже приведен код, 

import numpy as np
from scipy.integrate import solve_ivp
from scipy.fftpack import diff as psdiff

#RHS of equations
def f(t,u):
    uxx= psdiff(u[N:],period=L,order=2)
    du1dt=u[:N]
    du2dt =a**2*uxx
    dudt=np.append(du1dt,du2dt)
    return dudt

a=1
amin=-40;bmax=40
L=bmax-amin;N=256
deltax=L/N
x=np.arange(amin,bmax,deltax)

u01 = 2*np.cosh(x)**(-1)
u02=np.zeros(N)
# y0
inital=np.append(u01,u02)

sola1 = solve_ivp(f, t_span=[0,40],y0=inital,args=(a,))

import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.plot(x,sola1.y[:N,5])
plt.show()

Ниже приведен мой ожидаемый результат,

Ожидаемый результат .

Мой python код может выполняться, но Я не могу получить ожидаемый результат и не могу найти проблему. Ниже приведен результат из моего python кода, моего результата

-------- ---------------------Обновить---------------------------- ------------------ Я также пробую новый код, но все еще не могу решить 

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
from scipy.fftpack import diff as psdiff
from itertools import chain 
def lambdifide_odes(y,t,a):
    # uxx =- (1j)**2*k**2*u[:N]
    u1=y[::2]
    u2=y[1::2]
    dudt=np.empty_like(y)
    du1dt=dudt[::2]
    du2dt=dudt[1::2]

    du1dt=u2
    uxx=psdiff(u1,order=2,period=L)
    du2dt=a**2*uxx
    return dudt
a=1
amin=-40;bmax=40
L=bmax-amin;N=256
deltax=L/N
x=np.arange(amin,bmax,deltax)
u01 = 2*np.cosh(x)**(-1)
u02=np.zeros(N)
initial=np.array(list(chain.from_iterable(zip(u01,u02))))
t0=np.linspace(0,40,100)
sola1 = odeint(lambdifide_odes,y0=initial,t=t0,args=(a,))
fig, ax = plt.subplots()
ax.plot(x,sola1[20,::2])
plt.show()

1 Ответ

0 голосов
/ 03 апреля 2020

У вас есть небольшая проблема с дизайном вектора состояния и использованием этого в функции ODE. Общее намерение состоит в том, что u[:N] - это волновая функция, а u[N:] - ее производная по времени. Теперь вам нужна вторая пространственная производная волновой функции, поэтому вам нужно использовать 

uxx= psdiff(u[:N],period=L,order=2)

в тот момент, когда вы используете производную по времени, делая эту смешанную третью производную, которая не встречается в уравнении.

...