In [1]:
%reload_ext autoreload
%autoreload 2
In [2]:
from PDESolver import *
2D tests in periodic conditions¶
2D stationnary equation with psiOp()¶
In [3]:
# Symbolic variables
x, y, xi, eta = symbols('x y xi eta', real=True)
u = Function('u')(x, y)
# Symbolic equation
symbol = xi**2 + eta**2 + 1
equation = Eq(psiOp(symbol, u), -sin(x) * sin(y))
# Define exact solution
def u_exact(x, y):
return -np.sin(x) * np.sin(y) / 3
# Solver creation
solver = PDESolver(equation)
# Domain
Lx, Ly = 2 * np.pi, 2 * np.pi
Nx, Ny = 128, 128
# Setup with exact solution as source
solver.setup(Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, initial_condition=None)
# Stationary solution
u_num = solver.solve_stationary_psiOp()
# Accuracy test
solver.test(u_exact=u_exact, threshold=20, component='real')
solver.show_stationary_solution(u=u_num, component='real')
********************************* * Partial differential equation * *********************************
⎛ 2 2 ⎞
psiOp⎝η + ξ + 1, u(x, y)⎠ = -sin(x)⋅sin(y)
********************
* Equation parsing *
********************
Equation rewritten in standard form: sin(x)*sin(y) + psiOp(eta**2 + xi**2 + 1, u(x, y))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: sin(x)*sin(y) + psiOp(eta**2 + xi**2 + 1, u(x, y))
Analyzing term: sin(x)*sin(y)
--> Classified as source term
Analyzing term: psiOp(eta**2 + xi**2 + 1, u(x, y))
--> Classified as pseudo linear term (psiOp)
Final linear terms: {}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, eta**2 + xi**2 + 1)]
Source terms: [sin(x)*sin(y)]
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
symbol =
2 2
η + ξ + 1
⚠️ For psiOp, use interactive_symbol_analysis.
*******************************
* Solving the stationnary PDE *
*******************************
boundary condition: periodic
symbol =
2 2
η + ξ + 1
✅ Elliptic pseudo-differential symbol: inversion allowed.
Right inverse asymptotic symbol:
1
───────────
2 2
η + ξ + 1
⚡ Optimization: Symbol independent of x and y — direct product in 2D Fourier.
Testing a stationary solution.
Test error : 3.321e-16
2D stationnary equation with psiOp() depending on spatial variables¶
In [4]:
# 1) Define symbolic variables
x, y, xi, eta = symbols('x y xi eta', real=True)
u = Function('u')(x, y)
# 2) Define the 2D symbol p(x, y, xi, eta)
p_symbol = xi**2 + eta**2 + x**2 + y**2 + 1
# 3) Choose an exact solution u_exact(x, y)
u_exact_expr = exp(-x**2 - y**2)
def u_exact(x, y):
return np.exp(-x**2 - y**2)
# 4) Compute the right-hand side f_expr = P[u_exact] symbolically
f_expr = (-3*x**2 - 3*y**2 + 5)*u_exact_expr
# 5) Build the stationary equation: psiOp(p_symbol, u) = f_expr
equation_2d = Eq(psiOp(p_symbol, u), f_expr)
# 6) Initialize and setup the solver
solver2d = PDESolver(equation_2d)
solver2d.setup(Lx=10.0, Ly=10.0, Nx=64, Ny=64, initial_condition=None)
# 7) Solve stationary psiOp problem (first-order asymptotic inverse)
u_num2d = solver2d.solve_stationary_psiOp(order=1)
# 8) Test numerical vs exact
solver2d.test(u_exact=u_exact, threshold=3e-1, component='real')
*********************************
* Partial differential equation *
*********************************
2 2
⎛ 2 2 2 2 ⎞ ⎛ 2 2 ⎞ - x - y
psiOp⎝η + x + ξ + y + 1, u(x, y)⎠ = ⎝- 3⋅x - 3⋅y + 5⎠⋅ℯ
********************
* Equation parsing *
********************
Equation rewritten in standard form: -(-3*x**2 - 3*y**2 + 5)*exp(-x**2 - y**2) + psiOp(eta**2 + x**2 + xi**2 + y**2 + 1, u(x, y))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -(-3*x**2 - 3*y**2 + 5)*exp(-x**2 - y**2) + psiOp(eta**2 + x**2 + xi**2 + y**2 + 1, u(x, y))
Analyzing term: -(-3*x**2 - 3*y**2 + 5)*exp(-x**2 - y**2)
--> Classified as source term
Analyzing term: psiOp(eta**2 + x**2 + xi**2 + y**2 + 1, u(x, y))
--> Classified as pseudo linear term (psiOp)
Final linear terms: {}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(1, eta**2 + x**2 + xi**2 + y**2 + 1)]
Source terms: [-(-3*x**2 - 3*y**2 + 5)*exp(-x**2 - y**2)]
⚠️ Pseudo‑differential operator detected: all other linear terms have been rejected.
symbol =
2 2 2 2
η + x + ξ + y + 1
⚠️ For psiOp, use interactive_symbol_analysis.
*******************************
* Solving the stationnary PDE *
*******************************
boundary condition: periodic
symbol =
2 2 2 2
η + x + ξ + y + 1
✅ Elliptic pseudo-differential symbol: inversion allowed.
Right inverse asymptotic symbol:
4.0⋅ⅈ⋅η⋅y 4.0⋅ⅈ⋅x⋅ξ
──────────────────────── + ────────────────────────
2 2
⎛ 2 2 2 2 ⎞ ⎛ 2 2 2 2 ⎞
⎝η + x + ξ + y + 1⎠ ⎝η + x + ξ + y + 1⎠ 1
- ─────────────────────────────────────────────────── + ─────────────────────
2 2 2 2 2 2 2 2
η + x + ξ + y + 1 η + x + ξ + y + 1
⚙️ 2D Kohn-Nirenberg Quantification
Testing a stationary solution. Test error : 2.038e-01
2D transport equation¶
In [5]:
# Definition of the 2D transport equation
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
transport_eq = Eq(diff(u, t), -diff(u, x) - diff(u, y)) # Diagonal transport
# Creation of the solver with periodic conditions
solver = PDESolver(transport_eq)
# Configuration of the 2D domain and initial condition
L = 10.0 # Square domain [-5, 5] × [-5, 5]
N = 512 # Spatial resolution
Lt=3.0
solver.setup(
Lx=L, Ly=L,
Nx=N, Ny=N,
Lt=Lt, Nt=2000,
initial_condition=lambda x, y: np.exp(-(x**2 + y**2)), # 2D Gaussian
n_frames=100
)
# Exact solution with periodic wrapping
def u_exact(x, y, t):
Lx = Ly = L # Square domain
x_shift = (x - t + Lx/2) % Lx - Lx/2 # Periodic wrapping in x
y_shift = (y - t + Ly/2) % Ly - Ly/2 # Periodic wrapping in y
return np.exp(-(x_shift**2 + y_shift**2))
# Solving and validation
solver.solve()
# 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 *
*********************************
∂ ∂ ∂
──(u(t, x, y)) = - ──(u(t, x, y)) - ──(u(t, x, y))
∂t ∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x, y), t) + Derivative(u(t, x, y), x) + Derivative(u(t, x, y), y)
Expanded equation: Derivative(u(t, x, y), t) + Derivative(u(t, x, y), x) + Derivative(u(t, x, y), y)
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Analyzing term: Derivative(u(t, x, y), x)
Derivative found: Derivative(u(t, x, y), x)
--> Classified as linear
Analyzing term: Derivative(u(t, x, y), y)
Derivative found: Derivative(u(t, x, y), y)
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1, Derivative(u(t, x, y), x): 1, Derivative(u(t, x, y), y): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment: ⅈ⋅(kx + ky - ω) --- Symbolic symbol analysis --- symb_omega: 0 symb_k: 0 Raw characteristic equation: ⅈ⋅(kx + ky - ω) Temporal order from dispersion relation: 1 self.pseudo_terms = [] --- Solutions found --- [kx + ky] --- Final linear operator --- -ⅈ⋅(kx + ky) ***************** * CFL condition * ***************** CFL condition violated: dt = 0.0015, max allowed dt = 3.035639631116778e-05 ******************** * 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.166e-12
Closest available time to t_eval=0.75: 0.75 Test error t = 0.75: 4.030e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 4.030e-02
Closest available time to t_eval=2.25: 2.25 Test error t = 2.25: 4.030e-02
Closest available time to t_eval=3.0: 3.0 Test error t = 3.0: 4.029e-02
2D heat equation¶
In [6]:
# Define the symbols and the unknown function
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
# Define the 2D heat equation
eq = Eq(diff(u, t), diff(u, x, x) + diff(u, y, y))
# Initialize the solver
solver = PDESolver(eq)
Lt=2.0
# Configure the solver with the problem parameters
solver.setup(
Lx=2 * np.pi, Ly=2 * np.pi, # Spatial domain: [0, 2π] × [0, 2π]
Nx=512, Ny=512, # Spatial resolution: 512×512 points
Lt=Lt, Nt=100, # Total time: 1.0, number of time steps: 100
initial_condition=lambda x, y: np.sin(x) * np.sin(y), # Initial condition
n_frames=50
)
# Solve the equation
solver.solve()
# Exact solution at t = 1
def u_exact(x, y, t):
"""
Analytical solution of the 2D Poisson equation.
The solution is given by u(x, y, t) = sin(x) * sin(y) * exp(-2t).
"""
return np.sin(x) * np.sin(y) * np.exp(-2.0 * 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)]
# Animate the solution
ani = solver.animate(overlay='contour')
HTML(ani.to_jshtml())
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂ ∂
──(u(t, x, y)) = ───(u(t, x, y)) + ───(u(t, x, y))
∂t 2 2
∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Expanded equation: Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1, Derivative(u(t, x, y), (x, 2)): -1, Derivative(u(t, x, y), (y, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx + ky - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx + ky - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ ⎛ 2 2⎞⎤
⎣ⅈ⋅⎝- kx - ky ⎠⎦
--- Final linear operator ---
2 2
- kx - ky
***************** * 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.068e-03
Closest available time to t_eval=0.5: 0.48 Test error t = 0.48: 4.088e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.081e-02
Closest available time to t_eval=1.5: 1.48 Test error t = 1.48: 4.075e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 4.069e-02
********************* * Solution plotting * *********************
Out[6]:
2D heat equation with psiOp¶
In [7]:
# Define symbols
t, x, y, xi, eta = symbols('t x y xi eta', real=True)
u_func = Function('u')
u = u_func(t, x, y)
# Define the equation using psiOp
eq = Eq(diff(u, t), psiOp(-(xi**2 + eta**2), u))
# Create the solver instance
solver = PDESolver(eq)
# Total simulation time
Lt = 2.0
# Setup the solver
solver.setup(
Lx=2 * np.pi, Ly=2 * np.pi,
Nx=512, Ny=512,
Lt=Lt, Nt=100,
initial_condition=lambda x, y: np.sin(x) * np.sin(y),
n_frames=50
)
# Solve the PDE
solver.solve()
# Define the exact solution
def u_exact(x, y, t):
return np.sin(x) * np.sin(y) * np.exp(-2.0 * t)
# Test the result
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)
]
# Animate the solution
ani = solver.animate(overlay='front')
HTML(ani.to_jshtml())
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 2 2 ⎞
──(u(t, x, y)) = psiOp⎝- η - ξ , u(t, x, y)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: -psiOp(-eta**2 - xi**2, u(t, x, y)) + Derivative(u(t, x, y), t)
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -psiOp(-eta**2 - xi**2, u(t, x, y)) + Derivative(u(t, x, y), t)
Analyzing term: -psiOp(-eta**2 - xi**2, u(t, x, y))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(-1, -eta**2 - 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, -eta**2 - xi**2)]
✅ Time derivative coefficient detected: 1
symbol =
2 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: 3.068e-03
Closest available time to t_eval=0.5: 0.48 Test error t = 0.48: 4.088e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 4.081e-02
Closest available time to t_eval=1.5: 1.48 Test error t = 1.48: 4.075e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 4.069e-02
********************* * Solution plotting * *********************
Out[7]:
2D Schrödinger equation¶
In [8]:
# Define the symbols and the unknown function
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
# 2D Schrödinger equation: i ∂t u = ∂xx u + ∂yy u
eq = Eq(I * diff(u, t), (diff(u, x, x) + diff(u, y, y)))
# Create the solver
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=1)
# Domain
Lx = Ly = 20
Nx = Ny = 256
Lt = 2.0
Nt = 200
# Initial condition: Gaussian wave packet modulated e^{-x^2 - y^2} e^{i(x + y)}
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=lambda x, y: np.exp(-x**2 - y**2) * np.exp(1j * (x + y))
)
# Solve the equation
solver.solve()
# Exact solution
def u_exact(x, y, t):
"""
Analytical solution of the 2D Schrödinger equation.
The solution is given by:
u(x, y, t) = 1 / sqrt(1 + 4j * t)**2 * exp(i * (x + y - 2 * t)) * exp(-((x + 2t)^2 + (y + 2t)^2) / (1 + 4j * t)).
"""
return 1 / np.sqrt(1 + 4j * t)**2 * np.exp(1j * (x + y - 2 * t)) * np.exp(-((x + 2 * t)**2 + (y + 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.3, component='abs') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2
∂ ∂ ∂
ⅈ⋅──(u(t, x, y)) = ───(u(t, x, y)) + ───(u(t, x, y))
∂t 2 2
∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: I*Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Expanded equation: I*Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Analyzing term: I*Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): I, Derivative(u(t, x, y), (x, 2)): -1, Derivative(u(t, x, y), (y, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx + ky + ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx + ky + ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 2 2⎤
⎣- kx - ky ⎦
--- Final linear operator ---
⎛ 2 2⎞
-ⅈ⋅⎝- kx - ky ⎠
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.01, max allowed dt = 1.2078426318447326e-05
********************
* 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.786e-43
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.053e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.433e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.099e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.575e-01
2D Schrödinger equation with psiOp¶
In [9]:
# Define the symbols and the unknown function
t, x, y, xi, eta = symbols('t x y xi eta', real=True)
u_func = Function('u')
u = u_func(t, x, y)
# 2D Schrödinger equation: i ∂t u = ∂xx u + ∂yy u
eq = Eq(I * diff(u, t), psiOp(-(xi**2 + eta**2), u))
# Create the solver
solver = PDESolver(eq)
# Domain
Lx = Ly = 20
Nx = Ny = 256
Lt = 2.0
Nt = 200
# Initial condition: Gaussian wave packet modulated e^{-x^2 - y^2} e^{i(x + y)}
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=lambda x, y: np.exp(-x**2 - y**2) * np.exp(1j * (x + y))
)
# Solve the equation
solver.solve()
# Exact solution
def u_exact(x, y, t):
"""
Analytical solution of the 2D Schrödinger equation.
The solution is given by:
u(x, y, t) = 1 / sqrt(1 + 4j * t)**2 * exp(i * (x + y - 2 * t)) * exp(-((x + 2t)^2 + (y + 2t)^2) / (1 + 4j * t)).
"""
return 1 / np.sqrt(1 + 4j * t)**2 * np.exp(1j * (x + y - 2 * t)) * np.exp(-((x + 2 * t)**2 + (y + 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.3, component='abs') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
∂ ⎛ 2 2 ⎞
ⅈ⋅──(u(t, x, y)) = psiOp⎝- η - ξ , u(t, x, y)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: -psiOp(-eta**2 - xi**2, u(t, x, y)) + I*Derivative(u(t, x, y), t)
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -psiOp(-eta**2 - xi**2, u(t, x, y)) + I*Derivative(u(t, x, y), t)
Analyzing term: -psiOp(-eta**2 - xi**2, u(t, x, y))
--> Classified as pseudo linear term (psiOp)
Analyzing term: I*Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): I}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(-1, -eta**2 - 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, -eta**2 - xi**2)]
✅ Time derivative coefficient detected: I
symbol =
⎛ 2 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: 2.786e-43
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.053e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.379e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.037e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 2.277e-01
2D fractional diffusion equation¶
In [10]:
# Define the symbols and the equation
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
# Parameters
alpha = 1.5 # Fractional order
nu = 0.01 # Classical diffusion coefficient
# 2D fractional diffusion equation
eq = Eq(diff(u, t), -Op((kx**2 + ky**2)**(alpha/2), u))
# Create the solver
solver = PDESolver(eq, dealiasing_ratio=0.5)
# Domain and initial condition setup
Lx, Ly = 2 * np.pi, 2 * np.pi # Domain size
Nx, Ny = 128, 128 # Grid resolution
Lt = 1.0 # Total simulation time
Nt = 100 # Number of time steps
# Initial condition: sinusoidal wave in 2D
initial_condition = lambda x, y: np.sin(x) * np.sin(y)
# Setup the solver
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=initial_condition
)
# Solve the equation
solver.solve()
# Exact solution function at t = 1
def u_exact(x, y, t):
return np.sin(x) * np.sin(y) * np.exp(-t * (1**2 + 1**2)**(alpha/2))
# 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 *
*********************************
⎛ 0.75 ⎞
∂ ⎜⎛ 2 2⎞ ⎟
──(u(t, x, y)) = -Op⎝⎝kx + ky ⎠ , u(t, x, y)⎠
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: Op((kx**2 + ky**2)**0.75, u(t, x, y)) + Derivative(u(t, x, y), t)
Expanded equation: Op((kx**2 + ky**2)**0.75, u(t, x, y)) + Derivative(u(t, x, y), t)
Analyzing term: Op((kx**2 + ky**2)**0.75, u(t, x, y))
--> Classified as symbolic linear term (Op)
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1}
Final nonlinear terms: []
Symbol terms: [(1, (kx**2 + ky**2)**0.75)]
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
-ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 1.68179283050743*(kx**2)**0.75
Raw characteristic equation:
0.75
⎛ 2⎞
-ⅈ⋅ω + 1.68179283050743⋅⎝kx ⎠
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ 3/4⎤
⎢ ⎛ 2⎞ ⎥
⎣-1.68179283050743⋅ⅈ⋅⎝kx ⎠ ⎦
--- Final linear operator ---
3/4
⎛ 2⎞
-1.68179283050743⋅⎝kx ⎠
*****************
* 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: 2.452e-02
Closest available time to t_eval=0.25: 0.25 Test error t = 0.25: 2.175e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.203e-02
Closest available time to t_eval=0.75: 0.75 Test error t = 0.75: 2.235e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 2.274e-02
2D non-linear equation¶
In [11]:
# Symbolic variables
t, x, y = symbols('t x y')
u_func = Function('u')
u = u_func(t, x, y)
# Definition of the equation ∂u/∂t = Δu + u²
eq = Eq(diff(u, t), diff(u, x, 2) + diff(u, y, 2) + u**2)
# Creation of the solver
solver = PDESolver(eq)
# Square domain centered around 0
L = 10.0
N = 256
Lt = 1 # Reduced time for improved stability (nonlinear term)
Nt = 1000 # More time steps for precision
# Initial condition: t = 0
def initial_condition(x, y):
return np.exp(-(x**2 + y**2)) # Centered Gaussian
# ❌ No simple exact solution for this nonlinear equation.
# For validation: use numerical convergence
# But if you absolutely want an "exact solution", here is an approximate case:
# This is a heuristic approximation, valid for short durations:
def u_exact(x, y, t):
"""Empirical approximation for short duration"""
denom = 4 * t + 1
return 1 / denom * np.exp(-(x**2 + y**2) / denom)
# Solver setup
solver.setup(
Lx=L, Ly=L,
Nx=N, Ny=N,
Lt=Lt, Nt=Nt,
initial_condition=initial_condition
)
# Solving
solver.solve()
# Validation at different times
n_test = 4
test_results = [
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
∂ 2 ∂ ∂
──(u(t, x, y)) = u (t, x, y) + ───(u(t, x, y)) + ───(u(t, x, y))
∂t 2 2
∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: -u(t, x, y)**2 + Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Expanded equation: -u(t, x, y)**2 + Derivative(u(t, x, y), t) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Analyzing term: -u(t, x, y)**2
--> Classified as nonlinear
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1, Derivative(u(t, x, y), (x, 2)): -1, Derivative(u(t, x, y), (y, 2)): -1}
Final nonlinear terms: [-u(t, x, y)**2]
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2
kx + ky - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2
kx + ky - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ ⎛ 2 2⎞⎤
⎣ⅈ⋅⎝- kx - ky ⎠⎦
--- Final linear operator ---
2 2
- kx - ky
*****************
* 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: 4.060e-12
Closest available time to t_eval=0.25: 0.25 Test error t = 0.25: 8.510e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.323e-01
Closest available time to t_eval=0.75: 0.75 Test error t = 0.75: 1.612e-01
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.815e-01
2D wave equation¶
In [12]:
# Define the symbols and the unknown function
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
# 2D wave equation
eq = Eq(diff(u, t, t), diff(u, x, x) + diff(u, y, y))
# Initialize the solver
solver = PDESolver(eq, time_scheme='ETD-RK4')
# Parameters
Lx = Ly = 2 * np.pi
Nx = Ny = 512
Lt = 2.0
Nt = 400
# Configure the solver with the problem parameters
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=lambda x, y: np.sin(x) * np.sin(y),
initial_velocity=lambda x, y: np.zeros_like(x),
)
# Solve the equation
solver.solve()
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution: sin(x) sin(y) cos(ωt), with ω = sqrt(2)
def u_exact(x, y, t):
"""
Analytical solution of the 2D wave equation.
The solution is given by u(x, y, t) = sin(x) * sin(y) * cos(ωt),
where ω = sqrt(kx^2 + ky^2) = sqrt(1^2 + 1^2).
"""
omega = np.sqrt(1**2 + 1**2)
return np.sin(x) * np.sin(y) * np.cos(omega * t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=5, component='abs') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2 2
∂ ∂ ∂
───(u(t, x, y)) = ───(u(t, x, y)) + ───(u(t, x, y))
2 2 2
∂t ∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x, y), (t, 2)) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Expanded equation: Derivative(u(t, x, y), (t, 2)) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Analyzing term: Derivative(u(t, x, y), (t, 2))
Derivative found: Derivative(u(t, x, y), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), (t, 2)): 1, Derivative(u(t, x, y), (x, 2)): -1, Derivative(u(t, x, y), (y, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2 2
kx + ky - ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2 2
kx + ky - ω
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ ___________ ___________⎤
⎢ ╱ 2 2 ╱ 2 2 ⎥
⎣-╲╱ kx + ky , ╲╱ kx + ky ⎦
--- Final linear operator ---
2 2
- kx - ky
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.005, max allowed dt = 0.0030619460109133708
********************
* 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.741e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 9.994e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 2.534e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.466e-02
2D wave equation with a source term¶
In [13]:
# Define the symbols and the wave equation with a source term
t, x, y, omega = symbols('t x y omega')
u_func = Function('u')
u = u_func(t, x, y)
# Source term
source_term = cos(x) * cos(y) * (5 * cos(sqrt(3) * t) - sqrt(3) * sin(sqrt(3) * t))
# Equation with source term
eq = Eq(diff(u, t, t), diff(u, x, x) + diff(u, y, y) - u + source_term)
# Create the solver with ETD-RK4 scheme
solver = PDESolver(eq, time_scheme='ETD-RK4', dealiasing_ratio=2/3)
# Simulation parameters
Lt = 2.0 # Total simulation time
Lx, Ly = 20, 20 # Spatial domain size
Nx, Ny = 128, 128 # Number of grid points
Nt = 200 # Number of time steps
# Initial conditions
initial_condition = lambda x, y: np.cos(x) * np.cos(y)
initial_velocity = lambda x, y: np.zeros_like(x)
# Setup the solver
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=initial_condition,
initial_velocity=initial_velocity
)
# Solve the equation
solver.solve()
# Plot energy evolution
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution at any time t
def u_exact(x, y, t):
return np.cos(x) * np.cos(y) * np.cos(np.sqrt(3) * t)
# Automatic testing at multiple times
n_test = 4
test_set = [
solver.test(
u_exact=u_exact,
t_eval=i * Lt / n_test,
threshold=5, # Adjust threshold if necessary
component='real'
)
for i in range(n_test + 1)
]
*********************************
* Partial differential equation *
*********************************
2 2 2
∂ ∂ ∂
───(u(t, x, y)) = (-√3⋅sin(√3⋅t) + 5⋅cos(√3⋅t))⋅cos(x)⋅cos(y) - u(t, x, y) + ───(u(t, x, y)) + ───(u(t, x, y))
2 2 2
∂t ∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: -(-sqrt(3)*sin(sqrt(3)*t) + 5*cos(sqrt(3)*t))*cos(x)*cos(y) + u(t, x, y) + Derivative(u(t, x, y), (t, 2)) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Expanded equation: u(t, x, y) + sqrt(3)*sin(sqrt(3)*t)*cos(x)*cos(y) - 5*cos(x)*cos(y)*cos(sqrt(3)*t) + Derivative(u(t, x, y), (t, 2)) - Derivative(u(t, x, y), (x, 2)) - Derivative(u(t, x, y), (y, 2))
Analyzing term: u(t, x, y)
--> Classified as linear
Analyzing term: sqrt(3)*sin(sqrt(3)*t)*cos(x)*cos(y)
--> Classified as source term
Analyzing term: -5*cos(x)*cos(y)*cos(sqrt(3)*t)
--> Classified as source term
Analyzing term: Derivative(u(t, x, y), (t, 2))
Derivative found: Derivative(u(t, x, y), (t, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {u(t, x, y): 1, Derivative(u(t, x, y), (t, 2)): 1, Derivative(u(t, x, y), (x, 2)): -1, Derivative(u(t, x, y), (y, 2)): -1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: [sqrt(3)*sin(sqrt(3)*t)*cos(x)*cos(y), -5*cos(x)*cos(y)*cos(sqrt(3)*t)]
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2 2
kx + ky - ω + 1
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2 2
kx + ky - ω + 1
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ _______________ _______________⎤
⎢ ╱ 2 2 ╱ 2 2 ⎥
⎣-╲╱ kx + ky + 1 , ╲╱ kx + ky + 1 ⎦
--- Final linear operator ---
2 2
- kx - ky - 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: 2.204e-02
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 1.098e-01
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.931e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.080e-01
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.210e-01
2D wave equation with psiOp¶
In [14]:
# Define the symbols and the unknown function
t, x, y, xi, eta = symbols('t x y xi eta', real=True)
u_func = Function('u')
u = u_func(t, x, y)
# 2D wave equation
eq = Eq(diff(u, t, t), psiOp(-(xi**2 + eta**2), u))
# Initialize the solver
solver = PDESolver(eq)
# Parameters
Lx = Ly = 2 * np.pi
Nx = Ny = 512
Lt = 2.0
Nt = 400
# Configure the solver with the problem parameters
solver.setup(
Lx=Lx, Ly=Ly, Nx=Nx, Ny=Ny, Lt=Lt, Nt=Nt,
initial_condition=lambda x, y: np.sin(x) * np.sin(y),
initial_velocity=lambda x, y: np.zeros_like(x),
)
# Solve the equation
solver.solve()
solver.plot_energy()
solver.plot_energy(log=True)
# Exact solution: sin(x) sin(y) cos(ωt), with ω = sqrt(2)
def u_exact(x, y, t):
"""
Analytical solution of the 2D wave equation.
The solution is given by u(x, y, t) = sin(x) * sin(y) * cos(ωt),
where ω = sqrt(kx^2 + ky^2) = sqrt(1^2 + 1^2).
"""
omega = np.sqrt(1**2 + 1**2)
return np.sin(x) * np.sin(y) * np.cos(omega * t)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=0.2, component='abs') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2
∂ ⎛ 2 2 ⎞
───(u(t, x, y)) = psiOp⎝- η - ξ , u(t, x, y)⎠
2
∂t
********************
* Equation parsing *
********************
Equation rewritten in standard form: -psiOp(-eta**2 - xi**2, u(t, x, y)) + Derivative(u(t, x, y), (t, 2))
⚠️ psiOp detected: skipping expansion for safety
Expanded equation: -psiOp(-eta**2 - xi**2, u(t, x, y)) + Derivative(u(t, x, y), (t, 2))
Analyzing term: -psiOp(-eta**2 - xi**2, u(t, x, y))
--> Classified as pseudo linear term (psiOp)
Analyzing term: Derivative(u(t, x, y), (t, 2))
Derivative found: Derivative(u(t, x, y), (t, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), (t, 2)): 1}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: [(-1, -eta**2 - 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, -eta**2 - xi**2)]
✅ Time derivative coefficient detected: 1
symbol =
2 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.141e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 1.356e-01
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 3.670e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.339e-02
2D Klein-Gordon equation¶
In [15]:
# Physical parameters
c = 1.0 # Wave speed
m = 1.0 # Field mass
# Definition of the 2D Klein-Gordon equation
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
klein_gordon_eq = Eq(diff(u, t, t), c**2*(diff(u, x, x) + diff(u, y, y)) - m**2*u)
# Creation of the solver
solver = PDESolver(klein_gordon_eq)
# Domain configuration
L = 2 * np.pi # Square domain [0, 2π] × [0, 2π]
N = 512 # Grid points per dimension
T_final = 2.0 # Final time
Nt = 200 # Time steps
# Initial conditions
kx = 1
ky = 1
omega = np.sqrt(c**2*(kx**2 + ky**2) + m**2) # Direct numerical calculation
solver.setup(
Lx=L, Ly=L,
Nx=N, Ny=N,
Lt=T_final, Nt=Nt,
initial_condition=lambda x, y: np.sin(x) * np.sin(y),
initial_velocity=lambda x, y: np.zeros_like(x) # Initial time derivative is zero
)
# Replace the definition of the exact solution with:
omega_val = float(np.sqrt(c**2*(kx**2 + ky**2) + m**2)) # Convert to float
def u_exact(x, y, t):
return np.sin(x) * np.sin(y) * np.cos(omega_val * t)
# Solving and validation
solver.solve()
solver.plot_energy()
solver.plot_energy(log=True)
# Automatic testing
n_test = 4
test_set = [solver.test(u_exact=u_exact, t_eval=i * Lt / n_test, threshold=1e-1, component='real') for i in range(n_test + 1)]
*********************************
* Partial differential equation *
*********************************
2 2 2
∂ ∂ ∂
───(u(t, x, y)) = -1.0⋅u(t, x, y) + 1.0⋅───(u(t, x, y)) + 1.0⋅───(u(t, x, y))
2 2 2
∂t ∂x ∂y
********************
* Equation parsing *
********************
Equation rewritten in standard form: 1.0*u(t, x, y) + Derivative(u(t, x, y), (t, 2)) - 1.0*Derivative(u(t, x, y), (x, 2)) - 1.0*Derivative(u(t, x, y), (y, 2))
Expanded equation: 1.0*u(t, x, y) + Derivative(u(t, x, y), (t, 2)) - 1.0*Derivative(u(t, x, y), (x, 2)) - 1.0*Derivative(u(t, x, y), (y, 2))
Analyzing term: 1.0*u(t, x, y)
--> Classified as linear
Analyzing term: Derivative(u(t, x, y), (t, 2))
Derivative found: Derivative(u(t, x, y), (t, 2))
--> Classified as linear
Analyzing term: -1.0*Derivative(u(t, x, y), (x, 2))
Derivative found: Derivative(u(t, x, y), (x, 2))
--> Classified as linear
Analyzing term: -1.0*Derivative(u(t, x, y), (y, 2))
Derivative found: Derivative(u(t, x, y), (y, 2))
--> Classified as linear
Final linear terms: {u(t, x, y): 1.00000000000000, Derivative(u(t, x, y), (t, 2)): 1, Derivative(u(t, x, y), (x, 2)): -1.00000000000000, Derivative(u(t, x, y), (y, 2)): -1.00000000000000}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
2 2 2
1.0⋅kx + 1.0⋅ky - 1.0⋅ω + 1.0
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
2 2 2
1.0⋅kx + 1.0⋅ky - 1.0⋅ω + 1.0
Temporal order from dispersion relation: 2
self.pseudo_terms = []
--- Solutions found ---
⎡ _________________ _________________⎤
⎢ ╱ 2 2 ╱ 2 2 ⎥
⎣-╲╱ kx + ky + 1.0 , ╲╱ kx + ky + 1.0 ⎦
--- Final linear operator ---
2 2
- kx - ky - 1.0
*****************
* CFL condition *
*****************
CFL condition violated: dt = 0.01, max allowed dt = 0.0030619696478911772
********************
* 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: 3.068e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 2.144e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 9.016e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 1.206e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 1.349e-02
2D biharmonic equation¶
In [16]:
# Definition of the 2D biharmonic equation
t, x, y, kx, ky, omega = symbols('t x y kx ky omega')
u_func = Function('u')
u = u_func(t, x, y)
biharmonic_eq = Eq(diff(u, t), -(diff(u, x, 4) + 2*diff(u, x, 2, y, 2) + diff(u, y, 4)))
# Creation of the solver with periodic boundary conditions
solver = PDESolver(biharmonic_eq)
# Configuration of the domain
L = 2 * np.pi # Square domain [0, 2π] × [0, 2π]
N = 512 # Grid points per dimension
Lt = 2.0 # Final time
Nt = 400 # Time steps
# Initial sinusoidal 2D condition
initial_condition = lambda x, y: np.sin(x) * np.sin(y)
solver.setup(
Lx=L, Ly=L,
Nx=N, Ny=N,
Lt=Lt, Nt=Nt,
initial_condition=initial_condition
)
# Corresponding exact solution
def u_exact(x, y, t):
return np.sin(x) * np.sin(y) * np.exp(-4*t)
# Solving
solver.solve()
# 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 *
*********************************
4 4 4
∂ ∂ ∂ ∂
──(u(t, x, y)) = - ───(u(t, x, y)) - ───(u(t, x, y)) - 2⋅───────(u(t, x, y))
∂t 4 4 2 2
∂x ∂y ∂y ∂x
********************
* Equation parsing *
********************
Equation rewritten in standard form: Derivative(u(t, x, y), t) + Derivative(u(t, x, y), (x, 4)) + Derivative(u(t, x, y), (y, 4)) + 2*Derivative(u(t, x, y), (x, 2), (y, 2))
Expanded equation: Derivative(u(t, x, y), t) + Derivative(u(t, x, y), (x, 4)) + Derivative(u(t, x, y), (y, 4)) + 2*Derivative(u(t, x, y), (x, 2), (y, 2))
Analyzing term: Derivative(u(t, x, y), t)
Derivative found: Derivative(u(t, x, y), t)
--> Classified as linear
Analyzing term: Derivative(u(t, x, y), (x, 4))
Derivative found: Derivative(u(t, x, y), (x, 4))
--> Classified as linear
Analyzing term: Derivative(u(t, x, y), (y, 4))
Derivative found: Derivative(u(t, x, y), (y, 4))
--> Classified as linear
Analyzing term: 2*Derivative(u(t, x, y), (x, 2), (y, 2))
Derivative found: Derivative(u(t, x, y), (x, 2), (y, 2))
--> Classified as linear
Final linear terms: {Derivative(u(t, x, y), t): 1, Derivative(u(t, x, y), (x, 4)): 1, Derivative(u(t, x, y), (y, 4)): 1, Derivative(u(t, x, y), (x, 2), (y, 2)): 2}
Final nonlinear terms: []
Symbol terms: []
Pseudo terms: []
Source terms: []
*******************************
* Linear operator computation *
*******************************
Characteristic equation before symbol treatment:
4 2 2 4
kx + 2⋅kx ⋅ky + ky - ⅈ⋅ω
--- Symbolic symbol analysis ---
symb_omega: 0
symb_k: 0
Raw characteristic equation:
4 2 2 4
kx + 2⋅kx ⋅ky + ky - ⅈ⋅ω
Temporal order from dispersion relation: 1
self.pseudo_terms = []
--- Solutions found ---
⎡ ⎛ 4 2 2 4⎞⎤
⎣ⅈ⋅⎝- kx - 2⋅kx ⋅ky - ky ⎠⎦
--- Final linear operator ---
4 2 2 4
- kx - 2⋅kx ⋅ky - ky
*****************
* 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.068e-03
Closest available time to t_eval=0.5: 0.5 Test error t = 0.5: 6.168e-02
Closest available time to t_eval=1.0: 1.0 Test error t = 1.0: 6.143e-02
Closest available time to t_eval=1.5: 1.5 Test error t = 1.5: 6.118e-02
Closest available time to t_eval=2.0: 2.0 Test error t = 2.0: 6.093e-02
In [ ]: