首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >如何用solve_ivp解谱方法求解偏微分方程?

如何用solve_ivp解谱方法求解偏微分方程?
EN

Stack Overflow用户
提问于 2020-04-01 23:12:04
回答 1查看 552关注 0票数 0

我想用谱方法来解偏微分方程。这样的方程,初始条件是u(t=0,x)=(a^2)* formula (X),u'_t (t=0)=0。

为了解决这个问题,我使用了带有光谱方法的python。以下是代码:

代码语言:javascript
复制
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()

以下是我的预期结果:

expected result

我的python代码可以运行,但我无法获得预期的结果,也找不到我的python代码my result的结果是problem.Following

我也尝试了一个新的代码,但仍然不能解决

代码语言:javascript
复制
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()
EN

回答 1

Stack Overflow用户

发布于 2020-04-03 14:24:51

您在设计状态向量并在ODE函数中使用状态向量时遇到了一些小问题。总体意图是u[:N]是波函数,u[N:]是它的时间导数。现在您需要波函数的第二个空间导数,因此您需要使用

代码语言:javascript
复制
uxx= psdiff(u[:N],period=L,order=2)

此时,您使用时间导数,使其成为方程中不存在的混合三阶导数。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/60974226

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档