CODEWITHSUNDEEP
pythonnumpyscipynumerical-methodspde
user499244
user499244

Reputation: 21

how to use solve_ivp solve Partial Differential Equations with spectral method?

I want to use the spectral method to solve partial differential equations. The equations like that, formula,the initial condition is u(t=0,x)=(a^2)*sech(x),u'_t (t=0)=0.

To solve it, I use the python with the spectral method. Following is the code,

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()

Following is my expected result,

expected result.

My python code can run,but I can't get the expected result,and can't find the problem.Following is the result from my python code, my result

-----------------------------Update---------------------------------------------- I also try a new code,but still can't solve

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()

Upvotes: 2

Views: 1810

Answers (2)

Tarik
Tarik

Reputation: 11209

You got uxx calculated out of u[N:] instead of u[:N]. du1dt should be equal to u[:N] instead of u[N:]. I also changed the formula to calculate the second derivative of u. I matched the results you provided. Following is the code:

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

pde_def = r'$\frac{\partial^2u}{\partial t}=a^2\frac{\partial^2u}{\partial x^2}$' + '\n' + r'$u(t=0, x) = a^2 sech(x)$' + '\n' + r'$\frac{\partial u}{\partial t}_{(t=0)} = 0$'
        
#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=1000
assert N % 2 == 0
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))

for i in range(len(sola1.t)):
    fig, ax = plt.subplots(figsize=(20, 4))
    ax.plot(x,sola1.y[:N, i])
    plt.title(f'{pde_def}\na = {a}\nSolution at t = {sola1.t[i]}', fontsize=30)
    plt.ylim([0, 2])

Upvotes: 1

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

You have some slight problem with the design of your state vector and using this in the ODE function. The overall intent is that u[:N] is the wave function and u[N:] its time derivative. Now you want the second space derivative of the wave function, thus you need to use

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

at the moment you use the time derivative, making this a mixed third derivative that does not occur in the equation.

Upvotes: 1

Related Questions