mathematicskit.pde #

Returns:

PDESolution – extra["stability"] is check_cfl() for this step.

step(u, dt, scheme='upwind')[source]#

Advance grid values u by one step of scheme.

Parameters:

u (ndarray)
dt (float)
scheme (str)

Return type:

Returns:

ndarray, shape (n,)

class mathematicskit.pde.BurgersConservationLaw1D(u0, x_range=(0.0, 1.0), n=200, bc='outflow')[source]#

Bases: object

Inviscid Burgers’ equation \(u_t + (u^2/2)_x = 0\) on \([a, b]\), finite-volume schemes.

Parameters:

u0 (Callable | ndarray) – Initial cell values, sampled at the cell centers.
x_range (tuple of float)
n (int) – Number of cells.
bc (str) – "outflow" copies the edge cells into ghost cells (waves leave freely).

Examples

>>> import numpy as np
>>> law = BurgersConservationLaw1D(lambda x: np.where(x < 0.0, 1.0, 0.0), x_range=(-1.0, 1.0), n=200)
>>> sol = law.solve(0.8, dt=0.005, scheme="godunov")
>>> exact = burgers_riemann_solution(1.0, 0.0, sol.x, 0.8)
>>> bool(np.mean(np.abs(sol.final - exact)) < 0.01)
True

solve(t_final, dt, scheme='godunov', save_every=1)[source]#

Step from t = 0 to t_final.

Parameters:

t_final (float) – dt is shrunk if needed so a whole number of steps lands exactly on t_final.
dt (float) – dt is shrunk if needed so a whole number of steps lands exactly on t_final.
scheme (str)
save_every (int)

Return type:

Returns:

PDESolution – extra["stability"] is the CFL check with speed max |u0|.

step(u, dt, scheme='godunov')[source]#

Advance cell values u by one conservative step of scheme.

Parameters:

u (ndarray)
dt (float)
scheme (str)

Return type:

Returns:

ndarray, shape (n,)

class mathematicskit.pde.BurgersEquation1D(u0, length=6.283185307179586, n=64, nu=0.05)[source]#

Viscous Burgers’ equation \(u_t + u\,u_x = \nu\,u_{xx}\), periodic, Fourier pseudo-spectral.

Derivatives come from fourier_differentiation_matrix(); the nonlinear product \(u\,u_x\) is formed pointwise in physical space (the pseudo-spectral, or collocation, approach of Orszag 1969-1971 and Kreiss & Oliger 1972), without dealiasing. The resulting ODE system is integrated by mathematicskit.integrators through an @njit right-hand side that receives the flattened matrices in params.

Parameters:

u0 (Callable | ndarray)
length (float) – Period.
n (int) – Even number of grid points.
nu (float) – Viscosity.

Examples

>>> import numpy as np
>>> burgers = BurgersEquation1D(lambda x: np.sin(x), n=32, nu=0.1)
>>> sol = burgers.solve(1.0, dt=1e-3, save_every=100)
>>> bool(abs(sol.final.mean() - sol.u[0].mean()) < 1e-10)  # the mean is conserved
True

max_stable_dt(method='rk4')[source]#

Linearized RK4 estimate: the smaller of the diffusive and advective limits.

Diffusion contributes eigenvalues down to \(-\nu k_{\max}^2\) and advection (frozen at \(\max|u_0|\)) up to \(\pm i \max|u_0| k_{\max}\), with \(k_{\max} = \pi n / L\).

Return type:: float
Parameters:: method (str)

class mathematicskit.pde.EllipticSolution(x, u, y=None, method='', residual_norm=0.0, solver_result=None)[source]#

Bases: object

Output of a steady boundary-value solve (Laplace, Poisson, u'' = f).

u[i] (1D) or u[i, j] (2D, "ij" indexing) is the discrete solution at x[i] (and y[j]), boundary points included.

Parameters:

x (ndarray)
u (ndarray)
y (ndarray | None)
method (str)
residual_norm (float)
solver_result (object | None)

method: str = ''#

"direct", "cg", "jacobi", "gauss_seidel", "sor", or "chebyshev".

Type:: str

residual_norm: float = 0.0#

2-norm of the residual of the assembled linear system at the solution.

Type:: float

solver_result: object | None = None#

The full IterativeSolveResult (residual history, iteration count) when an iterative mathematicskit.linalg solver was used.

Type:: IterativeSolveResult, optional

u: ndarray#

Discrete solution.

Type:: ndarray, shape (n_x,) or (n_x, n_y)

x: ndarray#

Spatial grid along x.

Type:: ndarray, shape (n_x,)

y: ndarray | None = None#

Spatial grid along y (2D problems only).

Type:: ndarray, shape (n_y,), optional

class mathematicskit.pde.HeatEquation1D(u0, length=1.0, n=51, alpha=1.0, bc='dirichlet', boundary_values=(0.0, 0.0))[source]#

The 1D heat equation \(u_t = \alpha u_{xx}\) on \([0, L]\).

Parameters:

u0 (Callable | ndarray) – Initial temperature, u0(x) or its grid values.
length (float) – Rod length \(L\).
n (int) – Number of grid points (both endpoints included for Dirichlet; the right endpoint omitted for periodic).
alpha (float) – Diffusivity.
bc (str)
boundary_values (tuple of float) – Fixed end temperatures (u(0), u(L)) for bc="dirichlet".

Examples

>>> import numpy as np
>>> heat = HeatEquation1D(lambda x: np.sin(np.pi * x), n=41)
>>> sol = heat.solve(0.1, dt=1e-4)
>>> exact = np.sin(np.pi * sol.x) * np.exp(-np.pi**2 * 0.1)
>>> bool(np.max(np.abs(sol.final - exact)) < 1e-3)
True
>>> cn = heat.solve_theta(0.1, dt=1e-3, theta=0.5)
>>> bool(np.max(np.abs(cn.final - exact)) < 1e-3)
True

max_stable_dt(method='rk4')[source]#

Largest stable step: 2.785 dx^2 / (4 alpha) for "rk4", dx^2 / (2 alpha) for "ftcs", inf for implicit schemes.

Return type:: float
Parameters:: method (str)

solve_theta(t_final, dt, theta=0.5, save_every=1)[source]#

March the theta-method: FTCS (theta=0), Crank-Nicolson (0.5), backward Euler (1).

Parameters:

t_final (float)
dt (float)
theta (float) – Implicitness in [0, 1].
save_every (int) – Keep every save_every-th step (the last one is always kept).

Return type:

Returns:

PDESolution – method is "ftcs", "crank_nicolson", "backward_euler", or "theta=<value>"; extra["stability"] is check_diffusion_stability() for this step.

class mathematicskit.pde.HeatEquation2D(u0, lengths=(1.0, 1.0), n=(41, 41), alpha=1.0, boundary_value=0.0)[source]#

The 2D heat equation \(u_t = \alpha (u_{xx} + u_{yy})\) on a rectangle, fixed boundary temperature.

Semi-discretized with the five-point Laplacian and integrated by the method of lines through mathematicskit.integrators (solve()), or by the theta-method (solve_theta()).

Parameters:

u0 (Callable | ndarray) – Initial temperature u0(X, Y) (called on "ij" meshgrids) or its grid values.
lengths (tuple of float) – Rectangle side lengths (Lx, Ly); the domain is [0, Lx] x [0, Ly].
n (tuple of int) – Grid points (nx, ny), boundaries included.
alpha (float) – Diffusivity.
boundary_value (float) – Temperature held on the whole boundary.

Examples

>>> import numpy as np
>>> heat = HeatEquation2D(lambda X, Y: np.sin(np.pi * X) * np.sin(np.pi * Y), n=(21, 21))
>>> sol = heat.solve(0.02, dt=2e-4, save_every=10)
>>> sol.u.shape
(11, 21, 21)
>>> exact = np.exp(-2 * np.pi**2 * 0.02)
>>> bool(abs(sol.final[10, 10] - exact) < 5e-3)
True

max_stable_dt(method='rk4')[source]#

Largest stable step for "rk4" or "ftcs"; inf for implicit schemes.

Return type:: float
Parameters:: method (str)

solve_theta(t_final, dt, theta=0.5, save_every=1)[source]#

March the theta-method on the five-point Laplacian; see HeatEquation1D.solve_theta().

Return type:

Returns:

PDESolution – extra["stability"] is the step-size check from stability().

Parameters:

t_final (float)
dt (float)
theta (float)
save_every (int)

class mathematicskit.pde.MethodOfLinesPDE[source]#

Bases: ABC

Common base for a time-dependent PDE semi-discretized in space.

Concrete subclasses set self._rhs_njit (an @njit function with signature (state, t, params) -> dstate), self.params, and self.state0 (the flattened initial ODE state), and implement _to_field() to turn a batch of ODE states back into grid values (re-attaching boundary points, reshaping 2D grids).

max_stable_dt(method='rk4')[source]#

Largest step size for which method is linearly stable on this grid.

Subclasses override this; the default inf means “no known limit”.

Return type:: float
Parameters:: method (str)

params: ndarray = array([], dtype=float64)#

rhs(state, t=0.0)[source]#

Evaluate the semi-discrete right-hand side dU/dt at state.

Return type:

Parameters:

state (ndarray)
t (float)

solve(t_final, dt=None, method='rk4', save_every=1, **kwargs)[source]#

Integrate the semi-discrete system from t = 0 to t_final.

Parameters:

t_final (float)
dt (Optional[float]) – Fixed step (required for "rk4"), shrunk if needed so a whole number of steps lands exactly on t_final; initial step attempt for "dopri5".
method (str) – "rk4" or "dopri5" (subclasses may add more, e.g. "leapfrog").
save_every (int) – Keep every save_every-th time sample (the last one is always kept).
**kwargs – Forwarded to mathematicskit.integrators.dopri5_integrate().

Return type:

Returns:

PDESolution – extra["stability"] holds stability() for the fixed step dt (omitted for "dopri5", which controls its own step).

stability(dt, method='rk4')[source]#

Check dt against max_stable_dt() for method.

Parameters:

dt (float)
method (str)

Return type:

StabilityResult

Returns:

StabilityResult – number and limit are expressed as step sizes here.

state0: ndarray = array([], dtype=float64)#

x: ndarray = array([], dtype=float64)#

y: ndarray | None = None#

class mathematicskit.pde.PDESolution(t, x, u, y=None, method='', extra=<factory>)[source]#

Bases: object

Output of a time-dependent PDE solve, sampled on a fixed spatial grid.

For a 1D problem u[k, i] is the solution at time t[k] and position x[i]; for a 2D problem u[k, i, j] is the solution at (x[i], y[j]) ("ij" indexing). Boundary points are included.

Parameters:

t (ndarray)
x (ndarray)
u (ndarray)
y (ndarray | None)
method (str)
extra (dict)

extra: dict#

Free-form diagnostics, e.g. "stability" (a StabilityResult for the step size used) or "v" (the velocity field of a wave-equation solve).

Type:: dict

property final: ndarray#

The solution at the last saved time, u[-1].

Type:: ndarray

method: str = ''#

Time-stepping scheme used, e.g. "rk4", "dopri5", "leapfrog", "crank_nicolson", or "upwind".

Type:: str

snapshot(t)[source]#

Return the saved solution at the time sample nearest to t.

Parameters:: t (float)
Return type:: ndarray
Returns:: ndarray, shape (n_x,) or (n_x, n_y)

t: ndarray#

Saved time samples.

Type:: ndarray, shape (n_t,)

u: ndarray#

Solution at each saved time.

Type:: ndarray, shape (n_t, n_x) or (n_t, n_x, n_y)

x: ndarray#

Spatial grid along x.

Type:: ndarray, shape (n_x,)

y: ndarray | None = None#

Spatial grid along y (2D problems only).

Type:: ndarray, shape (n_y,), optional

class mathematicskit.pde.StabilityResult(number, limit, stable, scheme='', quantity='')[source]#

Bases: object

Outcome of a CFL / diffusion-number stability check for one step size.

Parameters:

number (float)
limit (float)
stable (bool)
scheme (str)
quantity (str)

limit: float#

The largest stable value of number for this scheme (inf for an unconditionally stable scheme, 0 for an unconditionally unstable one).

Type:: float

number: float#

The dimensionless number checked – the Courant number |c| dt / dx or the diffusion number alpha dt / dx^2.

Type:: float

quantity: str = ''#

"courant" or "diffusion".

Type:: str

scheme: str = ''#

The scheme the limit applies to.

Type:: str

stable: bool#

Whether number <= limit.

Type:: bool

class mathematicskit.pde.WaveEquation1D(u0, v0=None, length=1.0, n=101, c=1.0)[source]#

A vibrating string \(u_{tt} = c^2 u_{xx}\) on \([0, L]\) with \(u(0) = u(L) = 0\).

Parameters:

u0 (Callable | ndarray) – Initial displacement.
v0 (Union[Callable, ndarray, None]) – Initial velocity (default zero).
length (float)
n (int) – Grid points, both ends included.
c (float) – Wave speed.

Examples

>>> import numpy as np
>>> string = WaveEquation1D(lambda x: np.sin(np.pi * x), n=101)
>>> sol = string.solve(1.0, dt=0.01, method="leapfrog")  # Courant number 1
>>> sol.extra["stability"].stable
True
>>> bool(np.allclose(sol.final, -np.sin(np.pi * sol.x), atol=1e-3))  # half a period later
True

energy(solution)[source]#

Discrete energy \(E = \tfrac12 \sum_i v_i^2\,\Delta x + \tfrac12 c^2 \sum_i \big((u_{i+1} - u_i)/\Delta x\big)^2 \Delta x\) at each saved time.

Conserved by the exact PDE; leapfrog keeps it bounded (oscillating near its initial value), RK4 slowly dissipates it.

Parameters:: solution (PDESolution) – From solve().
Return type:: ndarray
Returns:: ndarray, shape (n_t,)

max_stable_dt(method='leapfrog')[source]#

dx / c for "leapfrog" (the CFL limit), \(\sqrt2\,dx/c\) for "rk4", inf otherwise.

Return type:: float
Parameters:: method (str)

solve(t_final, dt=None, method='leapfrog', save_every=1, **kwargs)[source]#

Integrate from t = 0 to t_final with "leapfrog" (default), "rk4", or "dopri5".

See MethodOfLinesPDE.solve(). The velocity field is returned in extra["v"] (same shape as u).

Return type:

Parameters:

t_final (float)
dt (float | None)
method (str)
save_every (int)

mathematicskit.pde.amplification_factor(scheme, number, xi)[source]#

Von Neumann amplification factor \(G(\xi)\) of a linear scheme.

With \(r\) the diffusion number and \(\nu\) the Courant number (advection speed \(c > 0\)), and \(s = \sin^2(\xi/2)\):

`"ftcs_heat"`	\(1 - 4 r s\)
`"backward_euler_heat"`	\(1 / (1 + 4 r s)\)
`"crank_nicolson_heat"`	\((1 - 2 r s) / (1 + 2 r s)\)
`"upwind"`	\(1 - \nu (1 - e^{-i\xi})\)
`"lax_friedrichs"`	\(\cos\xi - i \nu \sin\xi\)
`"lax_wendroff"`	\(1 - i \nu \sin\xi - \nu^2 (1 - \cos\xi)\)
`"ftcs_advection"`	\(1 - i \nu \sin\xi\)

See Strikwerda, Ch. 2.2, and LeVeque, Ch. 9.6 and 10.5.

Parameters:

scheme (str) – A key of AMPLIFICATION_SCHEMES.
number (float) – The scheme’s diffusion number \(r\) or Courant number \(\nu\).
xi (float or array-like) – Phase angle(s) \(\xi = k\,\Delta x\), usually in [0, pi].

Return type:

Returns:

ndarray of complex

Examples

>>> import numpy as np
>>> float(abs(amplification_factor("upwind", 1.0, np.pi)))
1.0
>>> float(abs(amplification_factor("ftcs_heat", 0.5, np.pi)))
1.0

mathematicskit.pde.burgers_cole_hopf_solution(u0, t, nu, length=6.283185307179586, n=256)[source]#

Exact periodic solution of viscous Burgers’ equation by the Hopf-Cole transformation.

Hopf (1950) and Cole (1951) showed that \(u = -2\nu\,\varphi_x/\varphi\) turns \(u_t + u u_x = \nu u_{xx}\) into the heat equation \(\varphi_t = \nu \varphi_{xx}\), with \(\varphi(x, 0) = \exp\!\big(-\tfrac{1}{2\nu}\int_0^x u_0\big)\). Here the antiderivative, the heat evolution (each Fourier mode times \(e^{-\nu k^2 t}\)), and \(\varphi_x\) are all computed spectrally with numpy.fft, so the result is exact up to the resolution of the n-point grid. \(\varphi\) spans roughly \(e^{\pm \max|\int u_0| / (2\nu)}\), so very small nu exhausts double precision.

Parameters:

u0 (callable or ndarray, shape (n,)) – Periodic initial data with zero mean (so that \(\varphi\) is periodic).
t (float)
nu (float)
length (float) – Period.
n (int) – Even number of grid points.

Return type:

tuple[ndarray, ndarray]

Returns:

x (ndarray, shape (n,))
u (ndarray, shape (n,))

Examples

>>> import numpy as np
>>> x, u = burgers_cole_hopf_solution(np.sin, t=1.0, nu=0.1, n=128)
>>> sol = BurgersEquation1D(np.sin, n=128, nu=0.1).solve(1.0, dt=1e-3)
>>> bool(np.max(np.abs(sol.final - u)) < 1e-10)
True

mathematicskit.pde.burgers_riemann_solution(u_left, u_right, x, t, x0=0.0)[source]#

Exact entropy solution of the Burgers Riemann problem with a jump at x0.

A shock moving at the Rankine-Hugoniot speed \(s = (u_L + u_R)/2\) if \(u_L > u_R\); otherwise a rarefaction fan \(u = (x - x_0)/t\) between \(u_L\) and \(u_R\).

Parameters:

u_left (float)
u_right (float)
x (array-like)
t (float) – Positive time.
x0 (float)

Return type:

Returns:

ndarray

Examples

>>> burgers_riemann_solution(1.0, 0.0, [0.4, 0.6], t=1.0).tolist()  # shock at x = 0.5
[1.0, 0.0]
>>> burgers_riemann_solution(0.0, 1.0, [-0.5, 0.25, 2.0], t=1.0).tolist()  # rarefaction
[0.0, 0.25, 1.0]

mathematicskit.pde.chebyshev_differentiation_matrix(n, a=-1.0, b=1.0)[source]#

Chebyshev collocation differentiation matrix on n Chebyshev points of [a, b].

\[D_{ij} = \frac{c_i}{c_j} \frac{(-1)^{i+j}}{x_i - x_j}\ (i \ne j), \qquad D_{ii} = -\sum_{j \ne i} D_{ij},\]

with \(c_0 = c_{n-1} = 2\) and \(c_i = 1\) otherwise; the “negative sum trick” for the diagonal makes \(D\) exact on constants and improves rounding (Trefethen 2000, Ch. 6, cheb.m). Points come from mathematicskit.numerical_analysis.chebyshev_nodes(), in increasing order.

Parameters:

n (int) – Number of points (polynomial degree n - 1).
a (float)
b (float)

Return type:

tuple[ndarray, ndarray]

Returns:

D (ndarray, shape (n, n))
x (ndarray, shape (n,))

Examples

>>> import numpy as np
>>> D, x = chebyshev_differentiation_matrix(12)
>>> bool(np.allclose(D @ x**5, 5 * x**4))  # exact for polynomials of degree < n
True

mathematicskit.pde.chebyshev_poisson_1d(f, n=32, a=-1.0, b=1.0, ua=0.0, ub=0.0)[source]#

Solve \(u'' = f\) on \([a, b]\) with \(u(a) = u_a\), \(u(b) = u_b\) by Chebyshev collocation.

Collocates \(D^2 u = f\) at the interior Chebyshev points and moves the known boundary columns to the right-hand side (Trefethen 2000, Ch. 7). For analytic f the error decays geometrically in n.

Parameters:

f (float | Callable)
n (int)
a (float)
b (float)
ua (float)
ub (float)

Return type:

Returns:

EllipticSolution

Examples

>>> import numpy as np
>>> sol = chebyshev_poisson_1d(lambda x: np.exp(x), n=16, ua=np.exp(-1), ub=np.exp(1))
>>> bool(np.max(np.abs(sol.u - np.exp(sol.x))) < 1e-12)
True

mathematicskit.pde.chebyshev_poisson_2d(f, n=24, a=-1.0, b=1.0)[source]#

Solve \(u_{xx} + u_{yy} = f\) on the square \([a, b]^2\) with \(u = 0\) on the boundary.

Tensor-product Chebyshev collocation: the interior operator is \(\tilde D^2 \otimes I + I \otimes \tilde D^2\) (Trefethen 2000, Ch. 7, program p16), solved densely with mathematicskit.linalg.lu_solve_system().

Parameters:

f (float | Callable) – Source term, f(X, Y) on "ij" meshgrids.
n (int) – Chebyshev points per direction.
a (float)
b (float)

Return type:

Returns:

EllipticSolution

Examples

>>> import numpy as np
>>> f = lambda X, Y: -2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
>>> sol = chebyshev_poisson_2d(f, n=20)
>>> X, Y = np.meshgrid(sol.x, sol.y, indexing="ij")
>>> bool(np.max(np.abs(sol.u - np.sin(np.pi * X) * np.sin(np.pi * Y))) < 1e-8)
True

mathematicskit.pde.check_cfl(c, dt, dx, scheme='upwind')[source]#

Check the Courant-Friedrichs-Lewy condition for an advection/wave scheme.

The limits (LeVeque 2007, Ch. 10) are \(\nu \le 1\) for "upwind", "godunov", "lax_friedrichs", "lax_wendroff" and the "leapfrog" (velocity-Verlet) wave scheme; forward-time centered-space "ftcs" advection is unstable for every \(\nu > 0\).

Parameters:

c (float)
dt (float)
dx (float)
scheme (str)

Return type:

StabilityResult

Returns:

StabilityResult

Examples

>>> check_cfl(1.0, 0.009, 0.01).stable
True
>>> check_cfl(1.0, 0.011, 0.01).stable
False

mathematicskit.pde.check_diffusion_stability(alpha, dt, dx, theta=0.0, ndim=1)[source]#

Check the diffusion-number limit of the theta-method for the heat equation.

The theta-method (\(\theta = 0\) explicit FTCS, \(\tfrac12\) Crank-Nicolson, 1 backward Euler) on the standard second-difference Laplacian in ndim dimensions (equal spacing) is stable iff

\[r \le \frac{1}{2\,\text{ndim}\,(1 - 2\theta)} \quad (\theta < \tfrac12),\]

and unconditionally stable for \(\theta \ge \tfrac12\) (Strikwerda, Ch. 6.3; LeVeque, Ch. 9.6).

Parameters:

alpha (float)
dt (float)
dx (float)
theta (float) – Implicitness, in [0, 1].
ndim (int) – Number of space dimensions.

Return type:

StabilityResult

Returns:

StabilityResult

Examples

>>> check_diffusion_stability(1.0, 0.004, 0.1).stable  # r = 0.4
True
>>> check_diffusion_stability(1.0, 0.006, 0.1).stable  # r = 0.6
False
>>> check_diffusion_stability(1.0, 1.0, 0.1, theta=0.5).limit
inf

mathematicskit.pde.courant_number(c, dt, dx)[source]#

Courant number \(\nu = |c|\,\Delta t / \Delta x\).

Parameters:

c (float) – Wave / advection speed.
dt (float) – Time step and grid spacing.
dx (float) – Time step and grid spacing.

Return type:

float

Returns:

float

Examples

>>> courant_number(2.0, 0.01, 0.04)
0.5

mathematicskit.pde.dalembert_solution(f, x, t, c=1.0, g_antiderivative=None, length=None)[source]#

D’Alembert’s solution of the wave equation (d’Alembert, 1747).

\[u(x, t) = \frac{f(x - ct) + f(x + ct)}{2} + \frac{G(x + ct) - G(x - ct)}{2c}\]

with \(f\) the initial displacement and \(G' = g\) the initial velocity. On the whole line this is exact as written. Passing length solves the fixed-end string on \([0, L]\) instead, by replacing \(f\) (and \(g\)) with their odd \(2L\)-periodic extensions (method of images; Strauss, Ch. 4.1).

Parameters:

f (Callable) – Initial displacement, vectorized over NumPy arrays.
x (array-like)
t (float)
c (float)
g_antiderivative (Optional[Callable]) – An antiderivative \(G\) of the initial velocity (on \([0, L]\) when length is given); zero velocity if omitted.
length (Optional[float]) – String length for the fixed-end problem.

Return type:

Returns:

ndarray, same shape as x

Examples

>>> import numpy as np
>>> bump = lambda s: np.exp(-100 * s**2)
>>> u = dalembert_solution(bump, np.array([-1.0, 0.0, 1.0]), t=1.0)
>>> np.round(u, 6).tolist()  # the bump splits into two half-height copies
[0.5, 0.0, 0.5]

mathematicskit.pde.diffusion_number(alpha, dt, dx)[source]#

Diffusion number \(r = \alpha\,\Delta t / \Delta x^2\).

Parameters:

alpha (float) – Diffusivity.
dt (float) – Time step and grid spacing.
dx (float) – Time step and grid spacing.

Return type:

float

Returns:

float

Examples

>>> round(diffusion_number(1.0, 0.001, 0.1), 12)
0.1

mathematicskit.pde.fem_poisson_1d(f, a=0.0, b=1.0, n_elements=10, ua=0.0, ub=0.0, nodes=None, quad_points=3)[source]#

Solve \(u'' = f\) on \([a, b]\), \(u(a) = u_a\), \(u(b) = u_b\), with P1 finite elements.

Parameters:

f (float | Callable) – Source term, vectorized over NumPy arrays if callable.
a (float) – Interval, ignored when nodes is given.
b (float) – Interval, ignored when nodes is given.
n_elements (int) – Number of equal elements, ignored when nodes is given.
ua (float) – Boundary values.
ub (float) – Boundary values.
nodes (Optional[ndarray]) – Strictly increasing mesh nodes (a nonuniform mesh), endpoints included.
quad_points (int) – Gauss-Legendre points per element for the load integrals (exact for f polynomial of degree <= 2 * quad_points - 2).

Return type:

Returns:

EllipticSolution – x are the mesh nodes and u the nodal values.

Examples

>>> import numpy as np
>>> nodes = np.array([0.0, 0.1, 0.35, 0.4, 0.8, 1.0])  # deliberately uneven
>>> sol = fem_poisson_1d(lambda x: 6 * x, nodes=nodes)  # exact u = x^3 - x
>>> bool(np.allclose(sol.u, nodes**3 - nodes))  # exact at the nodes
True

mathematicskit.pde.fourier_derivative(u, length=6.283185307179586, order=1)[source]#

Spectral derivative of periodic grid values via the FFT.

Multiplies each Fourier coefficient by \((i k)^{\text{order}}\) (numpy.fft.rfft() / irfft()). For an even number of points the unpaired Nyquist mode is zeroed for odd orders, so the derivative of real data stays real (Trefethen 2000, Ch. 3).

Parameters:

u (array-like, shape (n,)) – Samples at x_j = j * length / n, j = 0, ..., n-1.
length (float) – Period.
order (int) – Derivative order.

Return type:

Returns:

ndarray, shape (n,)

Examples

>>> import numpy as np
>>> x = np.linspace(0, 2 * np.pi, 16, endpoint=False)
>>> bool(np.allclose(fourier_derivative(np.sin(x)), np.cos(x)))
True

mathematicskit.pde.fourier_differentiation_matrix(n, length=6.283185307179586, order=1)[source]#

Dense Fourier spectral differentiation matrix on n periodic points (n even).

With \(h = 2\pi/n\) and \(d = i - j\),

\[D^{(1)}_{ij} = \tfrac12 (-1)^{d} \cot(d h / 2)\ (d \ne 0), \qquad D^{(2)}_{ij} = -\frac{(-1)^{d}}{2 \sin^2(d h/2)}\ (d \ne 0), \quad D^{(2)}_{ii} = -\frac{\pi^2}{3h^2} - \frac16,\]

rescaled by \((2\pi/\text{length})^{\text{order}}\) (Trefethen 2000, Ch. 3, eqs. (3.10)-(3.12)).

Parameters:

n (int) – Even number of grid points.
length (float) – Period.
order (int)

Return type:

Returns:

ndarray, shape (n, n)

Examples

>>> import numpy as np
>>> x = np.linspace(0, 2 * np.pi, 16, endpoint=False)
>>> D2 = fourier_differentiation_matrix(16, order=2)
>>> bool(np.allclose(D2 @ np.sin(3 * x), -9 * np.sin(3 * x)))
True

mathematicskit.pde.fourier_sine_coefficients(f, n_terms, length=1.0)[source]#

Fourier sine-series coefficients of f on \([0, L]\).

\[b_k = \frac{2}{L} \int_0^L f(x) \sin\!\left(\frac{k\pi x}{L}\right) dx, \qquad k = 1, \dots, n_{\text{terms}}\]

computed with scipy.integrate.quad(). See Fourier, Theorie analytique de la chaleur (1822), and Strauss, Ch. 5.1.

Parameters:

f (Callable[[float], float])
n_terms (int)
length (float)

Return type:

Returns:

ndarray, shape (n_terms,) – b[k - 1] is the coefficient of \(\sin(k\pi x/L)\).

Examples

>>> import numpy as np
>>> b = fourier_sine_coefficients(lambda x: np.sin(2 * np.pi * x), 3)
>>> np.round(b, 10).tolist()
[0.0, 1.0, 0.0]

mathematicskit.pde.heat_series_solution(coefficients, x, t, length=1.0, alpha=1.0)[source]#

Fourier’s series solution of the heat equation on a rod with zero end temperatures.

\[u(x, t) = \sum_{k \ge 1} b_k\, e^{-\alpha (k\pi/L)^2 t} \sin\!\left(\frac{k\pi x}{L}\right)\]

Parameters:

coefficients (array-like, shape (n_terms,)) – Sine coefficients of the initial condition, e.g. from fourier_sine_coefficients().
x (array-like)
t (float)
length (float)
alpha (float)

Return type:

Returns:

ndarray, same shape as x

Examples

>>> import numpy as np
>>> u = heat_series_solution([1.0], np.array([0.5]), t=0.1)
>>> bool(np.isclose(u[0], np.exp(-np.pi**2 * 0.1)))
True

mathematicskit.pde.laplacian_1d(n, dx, bc='dirichlet')[source]#

Second-difference matrix approximating \(d^2/dx^2\), \(O(dx^2)\).

\[(L u)_i = \frac{u_{i-1} - 2 u_i + u_{i+1}}{dx^2}\]

Parameters:

n (int) – Number of unknowns.
dx (float) – Grid spacing.
bc (str) – "dirichlet": the n unknowns are interior points and the (known) boundary values are dropped – add their contribution to the right-hand side separately. "periodic": u_{-1} = u_{n-1} and u_n = u_0. "neumann": the n unknowns include both endpoints, with zero slope imposed by reflecting ghost points (u_{-1} = u_1).

Return type:

csr_matrix

Returns:

scipy.sparse.csr_matrix, shape (n, n)

Examples

>>> L = laplacian_1d(4, 1.0)
>>> L.toarray().astype(int).tolist()
[[-2, 1, 0, 0], [1, -2, 1, 0], [0, 1, -2, 1], [0, 0, 1, -2]]

mathematicskit.pde.laplacian_2d(nx, ny, dx, dy)[source]#

Five-point Laplacian on an nx by ny grid of interior unknowns (Dirichlet).

Unknown u[i, j] is stored at flat index i * ny + j (row-major, matching u.ravel() for an "ij"-indexed array), so the operator is the Kronecker sum \(L_x \otimes I_{n_y} + I_{n_x} \otimes L_y\).

Parameters:

nx (int) – Number of interior unknowns along x and y.
ny (int) – Number of interior unknowns along x and y.
dx (float) – Grid spacings.
dy (float) – Grid spacings.

Return type:

csr_matrix

Returns:

scipy.sparse.csr_matrix, shape (nx * ny, nx * ny)

Examples

>>> import numpy as np
>>> L = laplacian_2d(3, 3, 1.0, 1.0)
>>> L.shape, float(L[4, 4])
((9, 9), -4.0)

mathematicskit.pde.max_amplification(scheme, number, n_xi=721)[source]#

Largest \(|G(\xi)|\) over \(\xi \in [0, \pi]\); the scheme is stable iff this is \(\le 1\).

Parameters:

scheme (str)
number (float)
n_xi (int) – Number of sampled phase angles.

Return type:

float

Returns:

float

Examples

>>> max_amplification("lax_wendroff", 0.8) <= 1.0
True
>>> max_amplification("ftcs_advection", 0.1) > 1.0
True

mathematicskit.pde.multigrid_poisson_2d(f, n=65, boundary=0.0, tol=1e-08, max_cycles=50, pre_smooth=2, post_smooth=2, omega=0.8)[source]#

Solve \(u_{xx} + u_{yy} = f\) on the unit square by multigrid V-cycles.

Discretization and boundary handling match solve_poisson_2d() (five-point Laplacian, Dirichlet data), so both converge to the same discrete solution.

Parameters:

f (float | Callable | ndarray) – Source term; a callable is evaluated as f(X, Y) on "ij" meshgrids.
n (int) – Grid points per side, boundaries included; must be 2**k + 1.
boundary (float | Callable | ndarray) – Dirichlet data.
tol (float) – Stop when the residual 2-norm falls below tol times its initial value.
max_cycles (int)
pre_smooth (int) – Weighted-Jacobi sweeps before and after each coarse-grid correction.
post_smooth (int) – Weighted-Jacobi sweeps before and after each coarse-grid correction.
omega (float) – Jacobi damping (4/5 is optimal for smoothing the 2D five-point Laplacian).

Return type:

Returns:

EllipticSolution – solver_result is an IterativeSolveResult whose residual_history holds the relative residual after each V-cycle.

Examples

>>> import numpy as np
>>> f = lambda X, Y: -2 * np.pi**2 * np.sin(np.pi * X) * np.sin(np.pi * Y)
>>> sol = multigrid_poisson_2d(f, n=33)
>>> sol.solver_result.converged, sol.solver_result.iterations < 15
(True, True)

mathematicskit.pde.solve_laplace_2d(boundary, n=(33, 33), x_range=(0.0, 1.0), y_range=(0.0, 1.0), method='direct', **kwargs)[source]#

Solve Laplace’s equation \(u_{xx} + u_{yy} = 0\) with Dirichlet data boundary.

Equivalent to solve_poisson_2d() with f = 0; see it for the parameters. By the discrete maximum principle, the solution’s extreme values lie on the boundary.

Return type:

Returns:

EllipticSolution

Parameters:

boundary (float | Callable | ndarray)
method (str)

Examples

>>> import numpy as np
>>> sol = solve_laplace_2d(lambda X, Y: X + 2 * Y, n=(9, 9))  # harmonic data is reproduced exactly
>>> X, Y = np.meshgrid(sol.x, sol.y, indexing="ij")
>>> bool(np.allclose(sol.u, X + 2 * Y))
True

mathematicskit.pde.solve_poisson_1d(f, a=0.0, b=1.0, n=101, ua=0.0, ub=0.0)[source]#

Solve \(u'' = f\) on \([a, b]\) with \(u(a) = u_a\), \(u(b) = u_b\).

Second-order central differences, \(O(\Delta x^2)\); the tridiagonal system is solved with scipy.sparse.linalg.spsolve().

Parameters:

f (float | Callable | ndarray) – Source term.
a (float)
b (float)
n (int) – Grid points, both ends included.
ua (float) – Boundary values.
ub (float) – Boundary values.

Return type:

Returns:

EllipticSolution

Examples

>>> import numpy as np
>>> sol = solve_poisson_1d(2.0, n=11)  # u'' = 2, u(0) = u(1) = 0  =>  u = x^2 - x
>>> bool(np.allclose(sol.u, sol.x**2 - sol.x))
True

mathematicskit.pde.solve_poisson_2d(f, n=(33, 33), x_range=(0.0, 1.0), y_range=(0.0, 1.0), boundary=0.0, method='direct', tol=1e-08, max_iter=20000, omega=None)[source]#

Solve \(u_{xx} + u_{yy} = f\) on a rectangle with Dirichlet boundary data.

Parameters:

f (float | Callable | ndarray) – Source term; a callable is evaluated as f(X, Y) on "ij" meshgrids.
n (tuple of int) – Grid points (nx, ny), boundaries included.
x_range (tuple of float)
y_range (tuple of float)
boundary (float | Callable | ndarray) – Dirichlet data; a callable is evaluated as g(X, Y), and for an array only the boundary entries are used.
method (str) – "direct" uses scipy.sparse.linalg.spsolve(); the others use ConjugateGradient, JacobiIteration, GaussSeidel, and SOR on the dense matrix.
tol (float) – Relative residual tolerance for the iterative methods.
max_iter (int) – Iteration cap for the iterative methods.
omega (float, optional) – SOR relaxation factor; defaults to optimal_sor_omega() (Young, 1950).

Return type: