我想用谱法求解偏微分方程。像这样的方程,,初始条件是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()
您的状态向量设计以及在 ODE 函数中使用它时遇到一些小问题。总体意图是
u[:N]
是波函数,u[N:]
是它的时间导数。现在你想要波函数的二阶空间导数,因此你需要使用
uxx= psdiff(u[:N],period=L,order=2)
此时您使用时间导数,使其成为方程中未出现的混合三阶导数。
您通过
uxx
而不是 u[N:]
计算出 u[N:]
。 du1dt
应等于 u[:N]
而不是 u[N:]
。我还改变了计算u
的二阶导数的公式。我匹配了你提供的结果。以下是代码:
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
#RHS of equations
def f(t, u, a, deltax):
N = len(u) // 2
uxx = np.gradient(np.gradient(u[:N]) / deltax) / deltax
du1dt=u[N:]
du2dt =a**2 * uxx
dudt=np.append(du1dt, du2dt)
return dudt
a=1
amin=-40; bmax=40
N=1001
deltax = (bmax - amin) / (N - 1)
x = np.linspace(amin, bmax, N)
t_eval = [0, 5, 10, 20]
u01 = 2*np.cosh(x)**(-1)
u02 = np.zeros(N)
# y0
inital=np.append(u01, u02)
sola1 = solve_ivp(f, t_span=[0, 40], t_eval=t_eval, y0=inital, args=(a, deltax))
ts = [0, 5, 10, 20]
for i, t in enumerate(ts):
fig, ax = plt.subplots(figsize=(10, 2))
ax.plot(x,sola1.y[:N, i])
plt.title(f't={t}')
plt.ylim([0, 2])