In [1]:
%reload_ext autoreload
%autoreload 2
In [2]:
from PDESolver import *
1D Tests in periodic conditions¶
1D stationnary equation with psiOp()¶
In [3]:
# Define symbolic variables
x, xi = symbols('x xi')
u = Function('u')(x)
# Define the PDE: psiOp[xi^2 + 1, u](x) = 0
equation = Eq(psiOp(xi**2 + 1, u), sin(x)) # rhs will be overwritten
# Define exact solution
def u_exact(x):
return np.sin(x)/2
# Create the solver
solver = PDESolver(equation)
# Setup parameters
Lx = 2 * np.pi
Nx = 256
solver.setup(Lx=Lx, Nx=Nx, initial_condition=None)
# Solve the stationary problem
u_num = solver.solve_stationary_psiOp(order=1)
# Exact solution on grid
u_ref = u_exact(solver.X)
# Compute relative L2 error
error = np.linalg.norm(np.real(u_num) - u_ref) / np.linalg.norm(u_ref)
solver.test(u_exact=u_exact, threshold=5e-3, component='abs')
solver.show_stationary_solution(u=u_num, component='real')
*********************************
* Partial differential equation *
*********************************
⎛ 2 ⎞
psiOp⎝ξ + 1, u(x)⎠ = sin(x)
********************
* Equation parsing *
********************
Equation rewritten in standard form: -sin(x) + psiOp(xi**2 + 1, u(x))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -sin(x) + psiOp(xi**2 + 1, u(x))
Analyzing term: -sin(x)
--> Classified as source term
Analyzing term: psiOp(xi**2 + 1, u(x))
--> Classified as pseudo linear term (psiOp)
Final linear terms: {}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, xi**2 + 1)]
Source terms: [-sin(x)]
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
symbol =
2
ξ + 1
⚠️ For psiOp, use interactive_symbol_analysis.
*******************************
* Solving the stationnary PDE *
*******************************
boundary condition: periodic
symbol =
2
ξ + 1
✅ Elliptic pseudo-differential symbol: inversion allowed.
Right inverse asymptotic symbol:
1
──────
2
ξ + 1
⚡ Optimization: symbol independent of x — direct product in Fourier.
Testing a stationary solution.
Test error : 2.151e-16
1D stationnary equation with psiOp() depending on $x$¶
In [4]:
# Define symbolic variables
x, xi = symbols('x xi', real=True)
u = Function('u')(x)
# Define the symbol
p_symbol = xi**2 + x**2 + 1
# Define the right-hand side f(x) = P[u_exact]
f_expr = 3 * (1 - x**2) * exp(-x**2)
# Define the full equation: P[u] = f
equation = Eq(psiOp(p_symbol, u), f_expr)
# Exact and initial
def u_exact(x): return np.exp(-x**2)
solver = PDESolver(equation)
solver.setup(Lx=30, Nx=1024, initial_condition=None)
u_num = solver.solve_stationary_psiOp(order=1)
solver.test(u_exact=u_exact, threshold=2, component='real')
********************************* * Partial differential equation * *********************************
2
⎛ 2 2 ⎞ ⎛ 2⎞ -x
psiOp⎝x + ξ + 1, u(x)⎠ = ⎝3 - 3⋅x ⎠⋅ℯ
********************
* Equation parsing *
********************
Equation rewritten in standard form: -(3 - 3*x**2)*exp(-x**2) + psiOp(x**2 + xi**2 + 1, u(x))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -(3 - 3*x**2)*exp(-x**2) + psiOp(x**2 + xi**2 + 1, u(x))
Analyzing term: -(3 - 3*x**2)*exp(-x**2)
--> Classified as source term
Analyzing term: psiOp(x**2 + xi**2 + 1, u(x))
--> Classified as pseudo linear term (psiOp)
Final linear terms: {}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, x**2 + xi**2 + 1)]
Source terms: [-(3 - 3*x**2)*exp(-x**2)]
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
symbol =
2 2
x + ξ + 1
⚠️ For psiOp, use interactive_symbol_analysis.
*******************************
* Solving the stationnary PDE *
*******************************
boundary condition: periodic
symbol =
2 2
x + ξ + 1
✅ Elliptic pseudo-differential symbol: inversion allowed.
Right inverse asymptotic symbol:
4.0⋅ⅈ⋅x⋅ξ 1
- ────────────── + ───────────
3 2 2
⎛ 2 2 ⎞ x + ξ + 1
⎝x + ξ + 1⎠
⚙️ 1D Kohn-Nirenberg Quantification
Testing a stationary solution. Test error : 1.827e-01
1D integro-differential equation¶
In [5]:
# Define the symbols and the integro-differential equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eps = 0.001
# Define the equation with an integral term using Op
integral_term = Op(1 / (I * (kx + eps)), u)
eq = Eq(diff(u, t), diff(u, x) + integral_term - u)
# Create the solver
# solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=2/3)
solver = PDESolver(eq, time_scheme='ETD-RK4')
# Simulation parameters
Lt = 2.0 # Total simulation time
Lx = 20 # Spatial domain size
Nx = 1024 # Number of grid points
Nt = 200 # Number of time steps
# Initial condition
initial_condition = lambda x: np.sin(x)
# Setup the solver
solver.setup(
Lx=Lx, Nx=Nx, Lt=Lt, Nt=Nt,
initial_condition=initial_condition
)
# Solve the equation
solver.solve()
# Exact solution at any time t
def u_exact(x, t):
return np.exp(-t) * np.sin(x)
# Automatic testing at multiple times
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=0.5, # Adjust threshold if necessary
component='real'
)
for i in range(n_test + 1)
]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ -ⅈ ⎞ ∂
──(u(t, x)) = -u(t, x) + Op⎜──────────, u(t, x)⎟ + ──(u(t, x))
∂t ⎝kx + 0.001 ⎠ ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: u(t, x) - Op(-I/(kx + 0.001), u(t, x)) + Derivative(u(t, x), t) - Derivative(u(t, x), x)
Expanded equation: u(t, x) - Op(-I/(kx + 0.001), u(t, x)) + Derivative(u(t, x), t) - Derivative(u(t, x), x)
Analyzing term: u(t, x)
--> Classified as linear
Analyzing term: -Op(-I/(kx + 0.001), u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x), x)
Derivative found: Derivative(u(t, x), x)
--> Classified as linear
Final linear terms: {u(t, x): 1, Derivative(u(t, x), t): 1, Derivative(u(t, x), x): -1}
Final nonlinear terms: []
Symbol terms: [(-1, -I/(kx + 0.001))]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅kx - ⅈ⋅ω + 1
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: I/(kx + 0.001)
Raw characteristic equation:
ⅈ
-ⅈ⋅kx - ⅈ⋅ω + 1 + ──────────
kx + 0.001
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2 ⎤
⎢- 1000.0⋅kx - kx - 1000.0⋅ⅈ⋅kx + 1000.0 - ⅈ⎥
⎢────────────────────────────────────────────⎥
⎣ 1000.0⋅kx + 1.0 ⎦
--- Final linear operator ---
⎛ 2 ⎞
-ⅈ⋅⎝- 1000.0⋅kx - kx - 1000.0⋅ⅈ⋅kx + 1000.0 - ⅈ⎠
──────────────────────────────────────────────────
1000.0⋅kx + 1.0
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.01, max allowed dt = 9.765625e-06
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.749e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.484e-01
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 3.472e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 4.210e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 4.840e-01
1D convolution equation¶
In [6]:
# Define symbols
t, x, kx, xi, omega, lam = symbols('t x kx xi omega lambda')
u_func = Function('u')
u = u_func(t, x)
lam = 1.0
# Define the kernel f(x) = exp(-lambda * |x|)
f_kernel = exp(-lam * Abs(x))
# Convolution equation: ∂t u = - OpConv(f_kernel, u)
eq = Eq(diff(u, t), -Op(fourier_transform(f_kernel, x, kx/(2*pi)), u))
# Create solver with lambda=1
solver = PDESolver(eq.subs(lam, 1), time_scheme='ETD-RK4')
# Simulation parameters
Lt = 2.0
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=100,
initial_condition=lambda x: np.cos(x)
)
# Solve
solver.solve()
# Exact solution: u(x, t) = cos(x) * exp(- (2*lambda / (lambda^2 + k^2)) * t )
def u_exact(x, t):
lam_val = 1.0
k_val = 1.0 # mode k=1 (cos(x))
decay = 2 * lam_val / (lam_val**2 + k_val**2)
return np.cos(x) * np.exp(-decay * t)
# Automatic testing at multiple times
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-3, component='abs')
for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 2.0 ⎞
──(u(t, x)) = -Op⎜───────, u(t, x)⎟
∂t ⎜ 2 ⎟
⎝kx + 1 ⎠
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op(2.0/(kx**2 + 1), u(t, x)) + Derivative(u(t, x), t)
Expanded equation: Op(2.0/(kx**2 + 1), u(t, x)) + Derivative(u(t, x), t)
Analyzing term: Op(2.0/(kx**2 + 1), u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1}
Final nonlinear terms: []
Symbol terms: [(1, 2.0/(kx**2 + 1))]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 2.0/(kx**2 + 1)
Raw characteristic equation:
2.0
-ⅈ⋅ω + ───────
2
kx + 1
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ -2.0⋅ⅈ ⎤
⎢─────────⎥
⎢ 2 ⎥
⎣kx + 1.0⎦
--- Final linear operator ---
-2.0
─────────
2
kx + 1.0
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
⚠️ Insufficient high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * ******************* Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 1.065e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.055e-04
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.048e-04
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.042e-04
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.040e-04
1D transport equation¶
In [7]:
# Define the symbols and the 1D transport equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eq = Eq(diff(u, t), -diff(u, x)) # Rightward transport at speed +1
# Create the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
Lt=5.0
# Domain setup and initial condition
solver.setup(
Lx=10, Nx=256, Lt=Lt, Nt=2000,
initial_condition=lambda x: np.exp(-x**2)
)
# Solving
solver.solve()
# Exact solution at t = 2.0 (with periodic wrapping)
def u_exact(x, t):
L = 10 # Domain length
return np.exp(-((x - t + L/2) % L - L/2)**2)
# Automatic testing
n_test = 5
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ∂
──(u(t, x)) = -──(u(t, x))
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), t) + Derivative(u(t, x), x)
Expanded equation: Derivative(u(t, x), t) + Derivative(u(t, x), x)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: Derivative(u(t, x), x)
Derivative found: Derivative(u(t, x), x)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), x): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
ⅈ⋅(kx - ω)
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
ⅈ⋅(kx - ω)
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
[kx]
--- Final linear operator ---
-ⅈ⋅kx
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.0025, max allowed dt = 0.00024285117048934224
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
⚠️ Insufficient high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
/home/fifi/pdesolver_prod/PDESolver.py:4234: RuntimeWarning: invalid value encountered in divide return np.where(np.abs(z) > 1e-12, (np.exp(z) - 1) / z, 1.0) /home/fifi/pdesolver_prod/PDESolver.py:4237: RuntimeWarning: invalid value encountered in divide return np.where(np.abs(z) > 1e-12, (np.exp(z) - 1 - z) / z**2, 0.5)
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 2.871e-12
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.749e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 4.749e-02
Closest available time to t_eval=3.0: 3.0 Test error t = 3.0: 4.749e-02
Closest available time to t_eval=4.0: 4.0 Test error t = 4.0: 4.761e-02
Closest available time to t_eval=5.0: 5.0 Test error t = 5.0: 4.788e-02
1D Fractional diffusion equation¶
In [8]:
# Define the symbols and the equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Fractional diffusion equation with alpha = 1.5 and added dissipation
alpha = 1.5
nu = 0.01 # Diffusion coefficient
# eq = Eq(diff(u, t), Op(abs(kx)**alpha, u) + nu * diff(u, x, 2))
eq = Eq(diff(u, t), -Op(abs(kx)**alpha, u))
# Create the solver
solver = PDESolver(eq, dealiasing_ratio=0.5, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=2000, # Reduced dt by increasing Nt
initial_condition=lambda x: np.sin(x)
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, t):
return np.sin(x) * np.exp(-t * abs(1)**alpha)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 1.5 ⎞
──(u(t, x)) = -Op⎝│kx│ , u(t, x)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op(Abs(kx)**1.5, u(t, x)) + Derivative(u(t, x), t)
Expanded equation: Op(Abs(kx)**1.5, u(t, x)) + Derivative(u(t, x), t)
Analyzing term: Op(Abs(kx)**1.5, u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1}
Final nonlinear terms: []
Symbol terms: [(1, Abs(kx)**1.5)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: Abs(kx)**1.5
Raw characteristic equation:
1.5
-ⅈ⋅ω + │kx│
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 3/2⎤
⎣-ⅈ⋅│kx│ ⎦
--- Final linear operator ---
3/2
-│kx│
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.888e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.771e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.651e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.530e-02
1D transport equation with Op¶
In [9]:
# Define symbols and equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Parameters
c = 1.0
beta = 1.3
# Equation
eq = Eq(diff(u, t) + c * diff(u, x), -Op(abs(kx)**beta, u))
# Create the solver
solver = PDESolver(eq, dealiasing_ratio=0.5, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=2000,
initial_condition=lambda x: np.cos(x)
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, t):
return np.cos(x - c * t) * np.exp(-t * abs(1)**beta)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=7e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ∂ ⎛ 1.3 ⎞
──(u(t, x)) + 1.0⋅──(u(t, x)) = -Op⎝│kx│ , u(t, x)⎠
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op(Abs(kx)**1.3, u(t, x)) + Derivative(u(t, x), t) + 1.0*Derivative(u(t, x), x)
Expanded equation: Op(Abs(kx)**1.3, u(t, x)) + Derivative(u(t, x), t) + 1.0*Derivative(u(t, x), x)
Analyzing term: Op(Abs(kx)**1.3, u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: 1.0*Derivative(u(t, x), x)
Derivative found: Derivative(u(t, x), x)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), x): 1.00000000000000}
Final nonlinear terms: []
Symbol terms: [(1, Abs(kx)**1.3)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
ⅈ⋅(1.0⋅kx - ω)
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: Abs(kx)**1.3
Raw characteristic equation:
1.3
ⅈ⋅(1.0⋅kx - ω) + │kx│
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found --- ⎡ 13⎤ ⎢ ──⎥ ⎢ 10⎥ ⎣kx - ⅈ⋅│kx│ ⎦ --- Final linear operator --- ⎛ 13⎞ ⎜ ──⎟ ⎜ 10⎟ -ⅈ⋅⎝kx - ⅈ⋅│kx│ ⎠ ***************** * CFL condition * ***************** CFL condition violated: dt = 0.001, max allowed dt = 9.587379924285257e-05 ******************** * Symbol condition * ******************** ✅ Spectral stability satisfied: Re(a(k)) ≤ 0 ✅ Proper high-frequency dissipation ✅ Reasonable spectral growth ✔ Symbol analysis completed. ******************* * Symbol plotting * *******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 1.065e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.474e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 2.282e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 3.492e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 6.492e-02
1D heat equation with periodic initial condition¶
In [10]:
# Define the symbols and the equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# 1D heat equation
eq = Eq(diff(u, t), diff(u, x, x))
# Create the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=100,
initial_condition=lambda x: np.sin(x)
)
# Solve the equation
solver.solve()
# Exact solution function
def u_exact(x, t):
return np.sin(x) * np.exp(-t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂ ∂
──(u(t, x)) = ───(u(t, x))
∂t 2
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))
Expanded equation: Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
kx - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2
kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-ⅈ⋅kx ⎦
--- Final linear operator ---
2
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * ******************* Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 6.131e-03
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.145e-03
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 6.167e-03
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 6.197e-03
1D heat equation with compact support initial condition¶
In [11]:
# Definition of symbols and equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u = Function('u')(t, x)
# 1D heat equation
eq = Eq(diff(u, t), diff(u, x, x))
# Creation of the solver with ETD-RK4 scheme
solver = PDESolver(eq, time_scheme='ETD-RK4')
# Simulation parameters
Lt = 2.0 # Total time
Lx = 5 * np.pi # Spatial domain
Nx = 256 # Grid points
Nt = 100 # Time steps
# Solver setup with Gaussian initial condition
solver.setup(
Lx=Lx,
Nx=Nx,
Lt=Lt,
Nt=Nt,
initial_condition=lambda x: np.exp(-2 * x**2) # Gaussian σ=0.5
)
# Solving the equation
solver.solve()
# Exact solution for the Gaussian with diffusion
def u_exact(x, t):
sigma = 0.5
variance = sigma**2 + 2 * t # Variance increases linearly
return (np.exp(-x**2 / (2 * variance)) /
np.sqrt(1 + 2 * t / sigma**2))
# Automated tests at different times
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=5e-3,
component='real'
)
for i in range(n_test + 1)
]
# Animate the solution
ani = solver.animate()
HTML(ani.to_jshtml())
*********************************
* Partial differential equation *
*********************************
2
∂ ∂
──(u(t, x)) = ───(u(t, x))
∂t 2
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))
Expanded equation: Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
kx - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2
kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-ⅈ⋅kx ⎦
--- Final linear operator ---
2
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 3.111e-53
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 7.390e-12
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.035e-07
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 2.867e-05
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.799e-04
********************* * Solution plotting * *********************
Out[11]:
1D heat equation with Op¶
In [12]:
# Define the symbols and the equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# 1D heat equation
#eq = Eq(diff(u, t), diff(u, x, 2))
eq = Eq(diff(u, t), Op((I*kx)**2, u)/2 + diff(u, x, 2)/2)
#eq = Eq(Op(-(I*omega),u), Op(-(I*kx)**2, u)/2 + diff(u, x, 2)/2)
# Create the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=2.0, Nt=100,
initial_condition=lambda x: np.sin(x)
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, t):
return np.sin(x) * np.exp(-t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂
───(u(t, x))
⎛ 2 ⎞ 2
∂ Op⎝-kx , u(t, x)⎠ ∂x
──(u(t, x)) = ───────────────── + ────────────
∂t 2 2
********************
* Equation parsing *
********************
Equation rewritten in standard form: -Op(-kx**2, u(t, x))/2 + Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))/2
Expanded equation: -Op(-kx**2, u(t, x))/2 + Derivative(u(t, x), t) - Derivative(u(t, x), (x, 2))/2
Analyzing term: -Op(-kx**2, u(t, x))/2
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))/2
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): -1/2}
Final nonlinear terms: []
Symbol terms: [(-1/2, -kx**2)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
kx
─── - ⅈ⋅ω
2
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: kx**2/2
Raw characteristic equation:
2
kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-ⅈ⋅kx ⎦
--- Final linear operator ---
2
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * ******************* Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
/home/fifi/pdesolver_prod/PDESolver.py:4234: RuntimeWarning: invalid value encountered in divide return np.where(np.abs(z) > 1e-12, (np.exp(z) - 1) / z, 1.0) /home/fifi/pdesolver_prod/PDESolver.py:4237: RuntimeWarning: invalid value encountered in divide return np.where(np.abs(z) > 1e-12, (np.exp(z) - 1 - z) / z**2, 0.5)
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 6.131e-03
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.145e-03
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 6.167e-03
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 6.197e-03
1D heat equation in Op formulation¶
In [13]:
# Define the symbols and the equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# 1D heat equation in Op formulation
# eq = Eq(Op((-I*omega), u), diff(u, x, x))
eq = Eq(Op((-I*omega), u), Op((-I*kx)**2, u))
# Create the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=512, Lt=Lt, Nt=400,
initial_condition=lambda x: np.sin(x),
initial_velocity=lambda x: np.zeros_like(x) # ∂u/∂t(x,0) = 0
)
# Solve the equation
solver.solve()
# Exact solution function
def u_exact(x, t):
return np.sin(x) * np.exp(-t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=2e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
⎛ 2 ⎞
Op(-ⅈ⋅ω, u(t, x)) = Op⎝-kx , u(t, x)⎠
********************
* Equation parsing *
********************
Equation rewritten in standard form: -Op(-kx**2, u(t, x)) + Op(-I*omega, u(t, x))
Expanded equation: -Op(-kx**2, u(t, x)) + Op(-I*omega, u(t, x))
Analyzing term: -Op(-kx**2, u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Op(-I*omega, u(t, x))
--> Classified as symbolic linear term (Op)
Final linear terms: {}
Final nonlinear terms: []
Symbol terms: [(-1, -kx**2), (1, -I*omega)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
0
--- Symbolic symbol analysis ---
symb_omega: -I*omega
symb_k: kx**2
Raw characteristic equation:
2
kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-ⅈ⋅kx ⎦
--- Final linear operator ---
2
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 2.169e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.516e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.504e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.492e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.480e-02
1D diffusion equation with psiOp¶
In [14]:
# Definition of the PDE: ∂u/∂t = - ∂²u/∂x² (diffusion equation)
u = Function('u')
t, x, xi = symbols('t x xi')
eq = Eq(diff(u(t,x), t), psiOp(-xi**2, u(t,x)))
# Solver configuration
solver = PDESolver(eq)
Lt=2.0
solver.setup(Lx=2*np.pi, Nx=128, Lt=Lt, Nt=200,
initial_condition=lambda x: np.sin(x))
# Solving and visualization
solver.solve()
# Exact solution function
def u_exact(x, t):
return np.sin(x) * np.exp(-t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 2 ⎞
──(u(t, x)) = psiOp⎝-ξ , u(t, x)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: -psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), t)
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), t)
Analyzing term: -psiOp(-xi**2, u(t, x))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(-1, -xi**2)]
Source terms: []
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
-ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = [(-1, -xi**2)]
✅ Time derivative coefficient detected: 1
symbol =
2
-ξ
⚠️ For psiOp, use interactive_symbol_analysis.
*******************
* Solving the PDE *
*******************
Closest available time to t_eval=0.0: 0.0
Test error t = 0.0: 1.734e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.852e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.771e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.745e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.786e-02
1D heat equation with Op¶
In [15]:
# Define symbols and equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Parameters
alpha = 1.8
nu = 0.1
# Equation
eq = Eq(diff(u, t), -Op(abs(kx)**alpha, u) + nu * diff(u, x, 2))
# Create the solver
solver = PDESolver(eq, dealiasing_ratio=0.5, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=2000,
initial_condition=lambda x: np.sin(x)
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, t):
return np.sin(x) * np.exp(-t * (abs(1)**alpha + nu * 1**2))
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=2e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂ ⎛ 1.8 ⎞ ∂
──(u(t, x)) = - Op⎝│kx│ , u(t, x)⎠ + 0.1⋅───(u(t, x))
∂t 2
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op(Abs(kx)**1.8, u(t, x)) + Derivative(u(t, x), t) - 0.1*Derivative(u(t, x), (x, 2))
Expanded equation: Op(Abs(kx)**1.8, u(t, x)) + Derivative(u(t, x), t) - 0.1*Derivative(u(t, x), (x, 2))
Analyzing term: Op(Abs(kx)**1.8, u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -0.1*Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): -0.100000000000000}
Final nonlinear terms: []
Symbol terms: [(1, Abs(kx)**1.8)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
0.1⋅kx - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: Abs(kx)**1.8
Raw characteristic equation:
2 1.8
0.1⋅kx - ⅈ⋅ω + │kx│
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ ⎛ 2 9/5⎞⎤
⎣ⅈ⋅⎝- 0.1⋅kx - │kx│ ⎠⎦
--- Final linear operator ---
2 9/5
- 0.1⋅kx - │kx│
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.958e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.695e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.434e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.192e-02
1D Burgers equation¶
In [16]:
# Definition of symbols and equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u = Function('u')(t, x)
# 1D Burgers equation
nu = 0.1
eq = Eq(diff(u, t), -u * diff(u, x) + nu * diff(u, x, x))
# Corrected exact solution
def phi(x, t):
return 2 + np.sin(x) * np.exp(-nu * t)
def dphi_dx(x, t):
return np.cos(x) * np.exp(-nu * t)
def u_exact(x, t):
return -2 * nu * dphi_dx(x, t) / phi(x, t)
# Creation of the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
Lt = 1.0
# Configuration with corrected initial condition
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=200,
initial_condition=lambda x: u_exact(x, 0)
)
# Solving
solver.solve()
# Automated tests
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=5e-2, # Adjusted threshold
component='real'
)
for i in range(n_test + 1)
]
*********************************
* Partial differential equation *
*********************************
2
∂ ∂ ∂
──(u(t, x)) = - u(t, x)⋅──(u(t, x)) + 0.1⋅───(u(t, x))
∂t ∂x 2
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: u(t, x)*Derivative(u(t, x), x) + Derivative(u(t, x), t) - 0.1*Derivative(u(t, x), (x, 2))
Expanded equation: u(t, x)*Derivative(u(t, x), x) + Derivative(u(t, x), t) - 0.1*Derivative(u(t, x), (x, 2))
Analyzing term: u(t, x)*Derivative(u(t, x), x)
--> Classified as nonlinear
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -0.1*Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): -0.100000000000000}
Final nonlinear terms: [u(t, x)*Derivative(u(t, x), x)]
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
0.1⋅kx - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2
0.1⋅kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-0.1⋅ⅈ⋅kx ⎦
--- Final linear operator ---
2
-0.1⋅kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 2.790e-03
Closest available time to t_eval=0.25: 0.25 Test error t = 0.25: 1.455e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.692e-02
Closest available time to t_eval=0.75: 0.75 Test error t = 0.75: 3.759e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.685e-02
1D equation with psiOp depending on spatial variable but not on $\xi$¶
In [17]:
# === Symbols ===
t, x, xi = symbols('t x xi', real=True)
u = Function('u')
# === Symbol dependent on x, periodic
symbol_expr = 1 + sin(x)
# === Equation ∂t u = -ψOp(symbol_expr, u) = - (1 + sin x) * u
eq = Eq(diff(u(t,x), t), -psiOp(symbol_expr, u(t,x)))
# === Solver creation
solver = PDESolver(eq)
# === Parameters
Lx = 2 * np.pi
Nx = 128
Lt = 2.0
Nt = 200
# === Initial condition
initial_condition = lambda x: np.cos(x)
# === Exact solution
def u_exact(x, t):
return np.cos(x) * np.exp(-t * (1 + np.sin(x)))
# === Setup
solver.setup(
Lx=Lx,
Nx=Nx,
Lt=Lt,
Nt=Nt,
boundary_condition='periodic',
initial_condition=initial_condition
)
# === Solving
solver.solve()
# === Tests
n_test = 4
for i in range(n_test + 1):
t_eval = i * Lt / n_test
solver.test(
u_exact=u_exact,
t_eval=t_eval,
threshold=5e-2,
component='real'
)
*********************************
* Partial differential equation *
*********************************
∂
──(u(t, x)) = -psiOp(sin(x) + 1, u(t, x))
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: psiOp(sin(x) + 1, u(t, x)) + Derivative(u(t, x), t)
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: psiOp(sin(x) + 1, u(t, x)) + Derivative(u(t, x), t)
Analyzing term: psiOp(sin(x) + 1, u(t, x))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, sin(x) + 1)]
Source terms: []
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
-ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = [(1, sin(x) + 1)]
✅ Time derivative coefficient detected: 1
symbol =
sin(x) + 1
⚠️ For psiOp, use interactive_symbol_analysis.
*******************
* Solving the PDE *
*******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.019e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 9.440e-03
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.165e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.267e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.203e-02
1D wave equation with periodic initial condition¶
In [18]:
# Define the symbols and the equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Wave equation
eq = Eq(diff(u, t, t), diff(u, x, x))
# Create the solver
solver = PDESolver(eq, time_scheme="ETD-RK4")
Lt=2.0
# Domain and initial conditions setup
solver.setup(
Lx=2 * np.pi, Nx=512, Lt=Lt, Nt=500,
initial_condition=lambda x: np.sin(x),
initial_velocity=lambda x: np.zeros_like(x) # ∂u/∂t(x,0) = 0
)
# Solving
solver.solve()
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution
def u_exact(x, t):
return np.sin(x) * np.cos(t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=2e-1, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂
───(u(t, x)) = ───(u(t, x))
2 2
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Expanded equation: Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), (t, 2)): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx - ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx - ω
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
[-kx, kx]
--- Final linear operator --- 2 -kx ***************** * CFL condition * ***************** ******************** * Symbol condition * ******************** ✅ Spectral stability satisfied: Re(a(k)) ≤ 0 ✅ Proper high-frequency dissipation ✅ Reasonable spectral growth ✔ Symbol analysis completed. ******************* * Symbol plotting * *******************
***************************** * Wave propagation analysis * *****************************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 2.169e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.022e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 2.173e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.530e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.070e-02
1D wave equation with compact support initial condition¶
In [19]:
# Define the symbols and the wave equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u = Function('u')(t, x)
# 1D wave equation: ∂²u/∂t² = ∂²u/∂x²
eq = Eq(diff(u, t, t), diff(u, x, x))
# Create the PDE solver with ETD-RK4 time integration
solver = PDESolver(eq, time_scheme="ETD-RK4")
# Simulation parameters
Lt = 10.0 # Total simulation time
Lx = 30.0 # Spatial domain size (increased for better Gaussian containment)
Nx = 1024 # Number of spatial grid points
Nt = 1000 # Number of time steps
# Gaussian initial condition parameters
sigma = 2.0 # Standard deviation of the Gaussian
x0 = 0.0 # Center position
# Initial condition: Gaussian pulse
initial_condition = lambda x: np.exp(-(x - x0)**2 / (2*sigma**2))
# Exact solution using d'Alembert's formula for wave equation
def u_exact(x, t):
return 0.5 * (np.exp(-(x - t - x0)**2 / (2*sigma**2)) +
np.exp(-(x + t - x0)**2 / (2*sigma**2)))
# Setup the solver with increased domain size
solver.setup(
Lx=Lx, Nx=Nx, Lt=Lt, Nt=Nt,
initial_condition=initial_condition,
initial_velocity=lambda x: np.zeros_like(x) # Zero initial velocity
)
# Solve the equation
solver.solve()
# Visualize energy conservation
solver.plot_energy()
solver.plot_energy(log=True)
# Perform automatic validation tests
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=1e-1, # Adjusted threshold for Gaussian decay
component='real'
)
for i in range(n_test + 1)
]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂
───(u(t, x)) = ───(u(t, x))
2 2
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Expanded equation: Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), (t, 2)): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx - ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx - ω
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
[-kx, kx]
--- Final linear operator ---
2
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
***************************** * Wave propagation analysis * *****************************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 1.360e-14
Closest available time to t_eval=2.5: 2.5 Test error t = 2.5: 3.380e-02
Closest available time to t_eval=5.0: 5.0 Test error t = 5.0: 3.143e-02
Closest available time to t_eval=7.5: 7.5 Test error t = 7.5: 3.108e-02
Closest available time to t_eval=10.0: 10.0 Test error t = 10.0: 3.247e-02
1D wave equation with a source term¶
In [20]:
# Define the symbols and the wave equation with a source term
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Source term
source_term = cos(x) * (3 * cos(sqrt(2) * t) - sqrt(2) * sin(sqrt(2) * t))
# Equation with source term
eq = Eq(diff(u, t, t), diff(u, x, x) - u + source_term)
# Create the solver with ETD-RK4 scheme
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=1)
# Simulation parameters
Lt = 2.0 # Total simulation time
solver.setup(
Lx=20, Nx=1024, Lt=Lt, Nt=400,
initial_condition=lambda x: np.cos(x),
initial_velocity=lambda x: np.zeros_like(x)
)
# Solving
solver.solve()
# Plot energy evolution
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution at t = 2: harmonic oscillation with frequency sqrt(2)
def u_exact(x, t):
return np.cos(x) * np.cos(np.sqrt(2) * t)
# Automatic testing
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=1, # Adjust threshold if necessary
component='real'
)
for i in range(n_test + 1)
]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂
───(u(t, x)) = (-√2⋅sin(√2⋅t) + 3⋅cos(√2⋅t))⋅cos(x) - u(t, x) + ───(u(t, x))
2 2
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: -(-sqrt(2)*sin(sqrt(2)*t) + 3*cos(sqrt(2)*t))*cos(x) + u(t, x) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Expanded equation: u(t, x) + sqrt(2)*sin(sqrt(2)*t)*cos(x) - 3*cos(x)*cos(sqrt(2)*t) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Analyzing term: u(t, x)
--> Classified as linear
Analyzing term: sqrt(2)*sin(sqrt(2)*t)*cos(x)
--> Classified as source term
Analyzing term: -3*cos(x)*cos(sqrt(2)*t)
--> Classified as source term
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {u(t, x): 1, Derivative(u(t, x), (t, 2)): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: [sqrt(2)*sin(sqrt(2)*t)*cos(x), -3*cos(x)*cos(sqrt(2)*t)]
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx - ω + 1
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx - ω + 1
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ _________ _________⎤
⎢ ╱ 2 ╱ 2 ⎥
⎣-╲╱ kx + 1 , ╲╱ kx + 1 ⎦
--- Final linear operator ---
2
- kx - 1
***************** * CFL condition * ***************** ******************** * Symbol condition * ******************** ✅ Spectral stability satisfied: Re(a(k)) ≤ 0 ✅ Proper high-frequency dissipation ✅ Reasonable spectral growth ✔ Symbol analysis completed. ******************* * Symbol plotting * *******************
***************************** * Wave propagation analysis * *****************************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 8.905e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 6.417e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.890e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 2.464e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.127e-01
1D wave equation with Op¶
In [21]:
# Define symbols and equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
# Parameters
gamma = 1.4
# Equation
eq = Eq(diff(u, t, 2) - diff(u, x, 2), -Op(abs(kx)**gamma, u))
# Create the solver
solver = PDESolver(eq, dealiasing_ratio=0.5, time_scheme='ETD-RK4')
Lt=2.0
# Domain and initial condition setup
solver.setup(
Lx=2 * np.pi, Nx=256, Lt=Lt, Nt=2000,
initial_condition=lambda x: np.sin(x),
initial_velocity=lambda x: np.zeros_like(x) # Correction ici
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, t):
return np.sin(x) * np.cos(np.sqrt(1 + abs(1)**gamma) * t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=0.2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂ ⎛ 1.4 ⎞
───(u(t, x)) - ───(u(t, x)) = -Op⎝│kx│ , u(t, x)⎠
2 2
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op(Abs(kx)**1.4, u(t, x)) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Expanded equation: Op(Abs(kx)**1.4, u(t, x)) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Analyzing term: Op(Abs(kx)**1.4, u(t, x))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), (t, 2)): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: [(1, Abs(kx)**1.4)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx - ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: Abs(kx)**1.4
Raw characteristic equation:
2 2 1.4
kx - ω + │kx│
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ _______________ _______________⎤
⎢ ╱ 2 1.4 ╱ 2 1.4 ⎥
⎣-╲╱ kx + │kx│ , ╲╱ kx + │kx│ ⎦
--- Final linear operator ---
2 1.4
- kx - │kx│
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
***************************** * Wave propagation analysis * *****************************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.523e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.416e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 3.260e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.805e-02
1D wave equation with psiOp¶
In [22]:
# Define the symbols
t, x, xi = symbols('t x xi', real=True)
u_func = Function('u')
u = u_func(t, x)
# Wave equation with psiOp
eq = Eq(diff(u, t, t), psiOp((-I*xi)**2, u))
# Create the solver
solver = PDESolver(eq)
Lt = 5.0
# Domain and initial conditions setup
solver.setup(
Lx=2 * np.pi, Nx=512, Lt=Lt, Nt=1000,
initial_condition=lambda x: np.sin(x),
initial_velocity=lambda x: np.zeros_like(x) # ∂u/∂t(x,0) = 0
)
# Solve
solver.solve()
# Plot energy
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution
def u_exact(x, t):
return np.sin(x) * np.cos(t)
# Automatic testing
n_test = 10
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=1, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂ ⎛ 2 ⎞
───(u(t, x)) = psiOp⎝-ξ , u(t, x)⎠
2
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: -psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), (t, 2))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), (t, 2))
Analyzing term: -psiOp(-xi**2, u(t, x))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), (t, 2)): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(-1, -xi**2)]
Source terms: []
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
-ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2
-ω
Temporal order from dispersion relation: 2
self.pseudo_terms = [(-1, -xi**2)]
✅ Time derivative coefficient detected: 1
symbol =
2
-ξ
⚠️ For psiOp, use interactive_symbol_analysis. ******************* * Solving the PDE * *******************
No energy data recorded. Call compute_energy() within solve(). No energy data recorded. Call compute_energy() within solve(). Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 2.169e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.526e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.972e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 6.340e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 9.970e-02
Closest available time to t_eval=2.5: 2.5 Test error t = 2.5: 3.567e-02
Closest available time to t_eval=3.0: 3.0 Test error t = 3.0: 1.146e-02
Closest available time to t_eval=3.5: 3.5 Test error t = 3.5: 1.818e-02
Closest available time to t_eval=4.0: 4.0 Test error t = 4.0: 5.193e-02
Closest available time to t_eval=4.5: 4.5 Test error t = 4.5: 2.080e-01
Closest available time to t_eval=5.0: 5.0 Test error t = 5.0: 1.534e-01
1D equation with psiOp depending on spatial potential¶
In [23]:
# Symbols
t, x, xi = symbols('t x xi', real=True)
u = Function('u')
# Symbol for the pseudo-differential operator (regular, x-dependent)
p_expr = xi**2 + cos(x)
# Define the PDE
equation = Eq(diff(u(t,x), t, t), -psiOp(p_expr, u(t,x)))
solver = PDESolver(equation)
# Setup the solver
Lx = 2 * np.pi
Nx = 512
Lt = 2.0
Nt = 1000
k = 5
omega = np.sqrt(k**2 + 0) # approximation si V(x) ≈ 0
initial_condition = lambda x: np.cos(k * x)
initial_velocity = lambda x: 0 * x
# Exact solution
u_exact = lambda x, t: np.cos(omega * t) * np.cos(k * x)
solver.setup(
Lx=Lx,
Nx=Nx,
Lt=Lt,
Nt=Nt,
boundary_condition='periodic',
initial_condition=initial_condition,
initial_velocity=initial_velocity
)
# Solve the PDE
solver.solve()
# Perform tests
for i in range(5):
solver.test(u_exact=u_exact, t_eval=i * Lt / 4, threshold=3e-1, component='real')
*********************************
* Partial differential equation *
*********************************
2
∂ ⎛ 2 ⎞
───(u(t, x)) = -psiOp⎝ξ + cos(x), u(t, x)⎠
2
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: psiOp(xi**2 + cos(x), u(t, x)) + Derivative(u(t, x), (t, 2))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: psiOp(xi**2 + cos(x), u(t, x)) + Derivative(u(t, x), (t, 2))
Analyzing term: psiOp(xi**2 + cos(x), u(t, x))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), (t, 2)): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, xi**2 + cos(x))]
Source terms: []
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2
-ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2
-ω
Temporal order from dispersion relation: 2
self.pseudo_terms = [(1, xi**2 + cos(x))]
✅ Time derivative coefficient detected: 1
symbol =
2
ξ + cos(x)
⚠️ For psiOp, use interactive_symbol_analysis.
*******************
* Solving the PDE *
*******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 4.700e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 7.342e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 2.809e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 2.035e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 7.893e-02
1D KdV equation¶
In [24]:
# Define the symbols and the KdV equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eq = Eq(diff(u, t) + 6 * u * diff(u, x) - diff(u, x, x, x), 0) # Standard KdV form
# Soliton parameters
c = 0.5 # Speed and amplitude of the soliton
x0 = 0.0 # Initial position
# Initial condition: centered soliton
initial_condition = lambda x: c / 2 * (1 / np.cosh(np.sqrt(c)/2 * (x - x0)))**2
# Create the solver with periodic boundary conditions and adapted dealiasing
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=2/3)
# Simulation parameters
Lx = 40
Nx = 2048
Lt = 5.0
Nt = 1000
# Setup
solver.setup(Lx=Lx, Nx=Nx, Lt=Lt, Nt=Nt, initial_condition=initial_condition)
# Solve
solver.solve()
# Exact solution at t = Lt
def u_exact(x, t):
return c / 2 * (1 / np.cosh(np.sqrt(c)/2 * (x - c*t - x0)))**2
# Automatic testing
n_test = 6
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=0.5, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
3
∂ ∂ ∂
6⋅u(t, x)⋅──(u(t, x)) + ──(u(t, x)) - ───(u(t, x)) = 0
∂x ∂t 3
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: 6*u(t, x)*Derivative(u(t, x), x) + Derivative(u(t, x), t) - Derivative(u(t, x), (x, 3))
Expanded equation: 6*u(t, x)*Derivative(u(t, x), x) + Derivative(u(t, x), t) - Derivative(u(t, x), (x, 3))
Analyzing term: 6*u(t, x)*Derivative(u(t, x), x)
--> Classified as nonlinear
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 3))
Derivative found: Derivative(u(t, x), (x, 3))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 3)): -1}
Final nonlinear terms: [6*u(t, x)*Derivative(u(t, x), x)]
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
⎛ 3 ⎞
ⅈ⋅⎝kx - ω⎠
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
⎛ 3 ⎞
ⅈ⋅⎝kx - ω⎠
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 3⎤
⎣kx ⎦
--- Final linear operator ---
3
-ⅈ⋅kx
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.005, max allowed dt = 2.346607975283368e-09
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
⚠️ Insufficient high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 5.816e-09
Closest available time to t_eval=0.8333333333333334: 0.8500000000000001 Test error t = 0.8500000000000001: 9.495e-02
Closest available time to t_eval=1.6666666666666667: 1.6500000000000001 Test error t = 1.6500000000000001: 1.760e-01
Closest available time to t_eval=2.5: 2.5 Test error t = 2.5: 2.556e-01
Closest available time to t_eval=3.3333333333333335: 3.35 Test error t = 3.35: 3.318e-01
Closest available time to t_eval=4.166666666666667: 4.15 Test error t = 4.15: 4.014e-01
Closest available time to t_eval=5.0: 5.0 Test error t = 5.0: 4.741e-01
1D Schrodinger equation¶
In [25]:
# Define the symbols and the 1D Schrödinger equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eq = Eq(diff(u, t), -I * diff(u, x, x))
# Create the solver with ETD-RK4 time scheme
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=1/2)
Lt=2.0
# Configure the domain and initial condition: modulated Gaussian packet
solver.setup(
Lx=20, Nx=1024, Lt=Lt, Nt=500,
initial_condition=lambda x: np.exp(-x**2) * np.exp(1j * x)
)
# Solve
solver.solve()
# Exact solution (one of the two propagated packets, here moving to the left)
def u_exact(x, t):
return 1 / np.sqrt(1 - 4j * t) * np.exp(1j * (x + t)) * np.exp(-((x + 2*t)**2) / (1 - 4j * t))
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=0.5, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂ ∂
──(u(t, x)) = -ⅈ⋅───(u(t, x))
∂t 2
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), t) + I*Derivative(u(t, x), (x, 2))
Expanded equation: Derivative(u(t, x), t) + I*Derivative(u(t, x), (x, 2))
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: I*Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 2)): I}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
⎛ 2 ⎞
ⅈ⋅⎝- kx - ω⎠
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
⎛ 2 ⎞
ⅈ⋅⎝- kx - ω⎠
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2⎤
⎣-kx ⎦
--- Final linear operator ---
2
ⅈ⋅kx
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.004, max allowed dt = 3.7745082245147894e-07
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
⚠️ Insufficient high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 5.367e-45
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 4.410e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.590e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 9.757e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.648e-01
1D Schrodinger equation with psiOp¶
In [26]:
# Define the symbols and the 1D Schrödinger equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eq = Eq(diff(u, t), -I * psiOp(-xi**2, u))
# Create the solver with ETD-RK4 time scheme
solver = PDESolver(eq)
Lt=2.0
# Configure the domain and initial condition: modulated Gaussian packet
solver.setup(
Lx=20, Nx=1024, Lt=Lt, Nt=500,
initial_condition=lambda x: np.exp(-x**2) * np.exp(1j * x)
)
# Solve
solver.solve()
# Exact solution (one of the two propagated packets, here moving to the left)
def u_exact(x, t):
return 1 / np.sqrt(1 - 4j * t) * np.exp(1j * (x + t)) * np.exp(-((x + 2*t)**2) / (1 - 4j * t))
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=0.5, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 2 ⎞
──(u(t, x)) = -ⅈ⋅psiOp⎝-ξ , u(t, x)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: I*psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), t)
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: I*psiOp(-xi**2, u(t, x)) + Derivative(u(t, x), t)
Analyzing term: I*psiOp(-xi**2, u(t, x))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(I, -xi**2)]
Source terms: []
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
-ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = [(I, -xi**2)]
✅ Time derivative coefficient detected: 1
symbol =
2
-ξ
⚠️ For psiOp, use interactive_symbol_analysis. ******************* * Solving the PDE * ******************* Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 5.367e-45
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 4.410e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.589e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 9.680e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.597e-01
1D Klein-Gordon equation¶
In [27]:
# Define the symbols and the Klein-Gordon equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
eq = Eq(diff(u, t, t), diff(u, x, x) - u)
# Create the solver with ETD-RK4 scheme
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=1/4)
Lt=2.0
# Simulation parameters
solver.setup(
Lx=20, Nx=1024, Lt=Lt, Nt=200,
initial_condition=lambda x: np.cos(x),
initial_velocity=lambda x: np.zeros_like(x)
)
# Solving
solver.solve()
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution at t = 2: harmonic oscillation with frequency sqrt(2)
def u_exact(x, t):
return np.cos(np.sqrt(2) * t) * np.cos(x)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂
───(u(t, x)) = -u(t, x) + ───(u(t, x))
2 2
∂t ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: u(t, x) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Expanded equation: u(t, x) + Derivative(u(t, x), (t, 2)) - Derivative(u(t, x), (x, 2))
Analyzing term: u(t, x)
--> Classified as linear
Analyzing term: Derivative(u(t, x), (t, 2))
Derivative found: Derivative(u(t, x), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x), (x, 2))
Derivative found: Derivative(u(t, x), (x, 2))
--> Classified as linear
Final linear terms: {u(t, x): 1, Derivative(u(t, x), (t, 2)): 1, Derivative(u(t, x), (x, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx - ω + 1
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx - ω + 1
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ _________ _________⎤
⎢ ╱ 2 ╱ 2 ⎥
⎣-╲╱ kx + 1 , ╲╱ kx + 1 ⎦
--- Final linear operator ---
2
- kx - 1
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.01, max allowed dt = 0.00975625926088619
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
***************************** * Wave propagation analysis * *****************************
******************* * Solving the PDE * *******************
Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 8.905e-04
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 6.145e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.619e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 2.391e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.125e-01
1D biharmonic equation¶
In [28]:
# Definition of the 1D biharmonic equation
t, x, kx, xi, omega = symbols('t x kx xi omega')
u_func = Function('u')
u = u_func(t, x)
biharmonic_eq = Eq(diff(u, t), -diff(u, x, x, x, x))
# Creation of the solver with periodic boundary conditions
solver = PDESolver(biharmonic_eq)
# Configuration of the domain and initial condition
Lx = 2 * np.pi
Nx = 256
Lt = 1.0
Nt = 200
solver.setup(
Lx=Lx, Nx=Nx,
Lt=Lt, Nt=Nt,
initial_condition=lambda x: np.sin(x)
)
# Solving
solver.solve()
# Exact solution
def u_exact(x, t):
return np.sin(x) * np.exp(-t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5e-2, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
4
∂ ∂
──(u(t, x)) = -───(u(t, x))
∂t 4
∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x), t) + Derivative(u(t, x), (x, 4))
Expanded equation: Derivative(u(t, x), t) + Derivative(u(t, x), (x, 4))
Analyzing term: Derivative(u(t, x), t)
Derivative found: Derivative(u(t, x), t)
--> Classified as linear
Analyzing term: Derivative(u(t, x), (x, 4))
Derivative found: Derivative(u(t, x), (x, 4))
--> Classified as linear
Final linear terms: {Derivative(u(t, x), t): 1, Derivative(u(t, x), (x, 4)): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
4
kx - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
4
kx - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 4⎤
⎣-ⅈ⋅kx ⎦
--- Final linear operator ---
4
-kx
*****************
* CFL condition *
*****************
********************
* Symbol condition *
********************
✅ Spectral stability satisfied: Re(a(k)) ≤ 0
✅ Proper high-frequency dissipation
✅ Reasonable spectral growth
✔ Symbol analysis completed.
*******************
* Symbol plotting *
*******************
******************* * Solving the PDE * ******************* Closest available time to t_eval=0.0: 0.0 Test error t = 0.0: 6.134e-03
Closest available time to t_eval=0.25: 0.25 Test error t = 0.25: 7.655e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 7.383e-03
Closest available time to t_eval=0.75: 0.75 Test error t = 0.75: 7.131e-03
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.903e-03
In [ ]: