Partial differential equations#
Partial differential equations expand ODEs to include multiple independant variables (coordinates).
Other spatial coordinates are typically extensions of the ODE case and are usually posed and treated as boundary value problems.
The time coordinate is almost always posed as an initial value problem and tends to introduces complexity in the form of numerical stability issues.
Classification of PDEs#
Introduction of time as a coordinate in PDEs raises the question of how information propogates, from the initial condition or other events that occur during solution.
Elliptic equations#
Equations of the form, $\( A\frac{\partial^2 u}{\partial x^2} + B\frac{\partial^2 u}{\partial x \partial y} + C\frac{\partial^2 u}{\partial y^2} \)$ are termed elliptic. The lack of a rate of change term implies that information is not being propogated, but rather is either instantaneously available everywere in space, or all propagations have finished and you are left with a steady state.
Note this does not imply the system isn’t changing with time (imagine the coefficients could be functions of time, \(A(t)\)), just that the impact of the change has reach an equilibrium.
Examples include steady-state, or quasistatic equations:
Laplace’s and and Poisson’s equation: \(\nabla^2 u = 0\), \(\nabla^2 u = f(x)\)
Electrostatics: Electric potential with a charge
Heat transport: Steady state heat distribution with a source
Fluid flow: Velocity potential of an incompressible, irrotational fluid with a sink
Helmholtz equation, \(\nabla^2 u +k^2 u = 0\)
Accoustics: vibrations of a membrane of a cavity
Quantum mechanics: wave functions in a potential well
Inclusion of Helmholtz’s equation may be surprising since it describes wave phenomena! Recall that wave propogation can be solved by the method of separation of variables to separate the time and space dependance, the latter of which is of the Helmholtz’s form (and defines the eigenfunctions!)
Solution approaches#
Elliptic equations generally amount to solving a system of (nonlinear) equations simultaneously, which we can do with root finders. Let’s recall the shortfalls of rootfinders in this context:
A good initial guess: This can be a substantial problem in common systems. A mitigation strategy might be to parametrically ramp-up a parameter from an easily solvable condition to your desired condition. E.g.: gradually ‘turn up the heat’ on a natural convection problem.
Convergence: You may need to discretize the problem with a large mesh which is a high dimensional problem. The root finder may wander easily in this space and fall into local minima or step completely out of the feasible region. Trust-region methods can help mitigate this but consume additional resources.
There is another subtle issue with parallelization on modern high performance platforms. A common approach is to perform domain decomposition which partitions the domain geometrically and passes pieces to different nodes. Elliptic systems imply each node instantaneously talks to all the others, which can be an overhead / node synchronization nightmare! For this reason multigrid methods / preconditioners are tremendously effective.
Parabolic equations#
Equations of the form, $\( A \frac{\partial u}{\partial t} - C \frac{\partial^2 u}{\partial x^2} \)$
propogate information in one direction (forward in time) and the process is diffusive.
Diffusive processes tend to smooth out discontinuities and any small oscillations, which is a favourable feature but doesn’t gaurantee stability as we will see.
This still implies an infinite rate of propogation of information, but now in the flux which instantaneously knows about a force everywhere at once.
##Solution schemes
Solving diffusion problems generally amounts to a time-marching scheme in which time is discretized with the spactial equation solved at each time step. We still have to be wary of stiffness!
It is common to want to solve the inverse problem, in which one uses a result to attempt to back-calculate a parameter of the model (e.g.: Use the rate of cooling to infer the thermal conductivity). It is here that the tendancy of diffusion to smooth the features (information only flows forward in time) works against you. Many paths can lead to the same result (consider trying to sharpen a blurry image!) this can be very problematic to arrive at anything more than a statistical answer.
Hyperbolic equations#
Equations of the form, $\( A \frac{\partial^2 u}{\partial t^2} - C \frac{\partial^2 u}{\partial x^2} \)\( propogate information forward *and* backward in time. This implies time-reversal symmetry, wherein you could replace \)t\( with \)-t\( and the physics would remain unchanged. This describes waves (second order equations) and advective terms (first order equations) e.g.: \)A \frac{\partial u}{\partial t} - C \frac{\partial u}{\partial x}$
Solution schemes#
Second time derivatives can generally be decomposed into a series of parabolic equations, which is a standard approach. Systems with periodic boundary conditions are treatable with fourier spectral analysis which can be much more efficient.
Time reversal symmetry implies the model may be reversible which is interesting for the inverse problem and reconstruction of earlier times before what is measured / simulated.
#Solving time dependent PDEs
The general approach to solving time dependent PDEs is to transform them into a BVP through a time-marching scheme:
Discretize the time derivative as an IVP into current and previous solution(s) with the initial value known
Discretize space as a BVP and add the term from the time descritization.
Step forward in time, updating the BVP accordingly.
Consider the parabolic case, since we can usually use reduction of order to transform a hyperbolic equation into parabolic equations.
where \(\alpha\) is the thermal diffusivity, which is necessarily postive.
The finite difference scheme for the spatial component is: $\(\alpha \frac{\partial^2 T}{\partial x^2} = \alpha \frac{T_{i-1} - 2 T_i + T_{i+1}}{\Delta x^2}\)$
Notation#
The spatial index will be noted as a subscript and indexed with \(i\), e.g.: \(T_i\).
The temporal index will be index with superscript \(t\).
The known timestep will always be \(T^t\) and we will be calculating \(T^{t+1}\).
Explicit Euler timestepping#
As in the IVP section, the explicit Euler scheme calculates the next solution based on the current one, $\( \frac{\partial T}{\partial t} = \frac{T^{t+1} -T^{t}}{\Delta t}\)$
and the full equation becomes,
with $\( r = \frac{\alpha \Delta t} {\Delta x^2}\)$ for compactness.
As a stencil for forward Euler time stepping is:
This implies that the value at a point at the next timestep depends only on the values at the previous time-step.
Example of Forward Euler time stepping#
# prompt: Give me an example of forward euler timestepping for a 1D heat balance equation with a slider for the time step size
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider
def forward_euler_heat(T, alpha, dx, dt, num_steps):
"""
Solves the 1D heat equation using forward Euler timestepping.
Args:
T: Initial temperature distribution (numpy array).
alpha: Thermal diffusivity.
dx: Spatial step size.
dt: Time step size.
num_steps: Number of time steps.
Returns:
A list of temperature distributions at each time step.
"""
T_history = [T.copy()]
for _ in range(num_steps):
T_new = T.copy()
for i in range(1, len(T) - 1):
T_new[i] = T[i] + alpha * dt / (dx ** 2) * (T[i - 1] - 2 * T[i] + T[i + 1])
T = T_new
T_history.append(T.copy())
return T_history
def plot_heat_equation(dt):
"""
Plots the solution of the heat equation for a given time step size.
Args:
dt: Time step size.
"""
# Parameters
alpha = 0.1 # Thermal diffusivity
dx = 0.1 # Spatial step size
num_steps = 50
x = np.arange(0, 1, 1/num_steps)
T_initial = np.zeros_like(x)
T_initial[int(len(x) / 2)] = 1 # Initial heat source in the middle
# Solve the heat equation
T_history = forward_euler_heat(T_initial, alpha, dx, dt, num_steps)
# Plotting
plt.figure(figsize=(8, 6))
for i in range(0, len(T_history), 5):
plt.plot(x, T_history[i], label=f"Time step {i}")
plt.xlabel("Position (x)")
plt.ylabel("Temperature (T)")
plt.title("Forward Euler Solution of 1D Heat Equation")
plt.legend()
plt.grid(True)
plt.show()
# Interactive widget
interact(plot_heat_equation, dt=FloatSlider(min=0.001, max=0.05, step=0.001, value=0.01));
Error analysis#
Analysing the error using the Explicit Euler timestepper is more sophisticated than for the ODE since we now need to consider the behaviouir of a (spatial) function at a given time.
We can do this through von Neumann stability analysis which checks the stability of a solution to oscillations. Mathematically, we express the soultion as a Fourier series and see which waves damp and which waves grow.
Consider Fourier modes of the solution: \(^kT_i^t = e^{j k i \Delta x}\) where \(k\) is the wave number and \(j\) is the imaginary number. For a particular wave with number \(k\), we check how the future solution, \(^kT_i^{t+1}\) depends on the current, \(^kT_i^n\), and define the growth factor \(G\), \(^kT_i^{t+1} = G ^kT_i^{t}.\) A method is deemed stable if the amplitude of the oscillations don’t grow over time; i.e.: \(|G| \le 1\).
Substitute this into the update equation:
Looking at the cosign, we can see \(G\) is in the range \(\bigg[1-4\frac{\alpha \Delta t}{\Delta x^2}, 1\bigg]\). The maximum already satisfies the stability, but the minimum must satisfy, $\( \frac{\alpha \Delta t}{\Delta x^2} \le \frac{1}{2}\)$
or more simply $\( \Delta t \le \frac{\Delta x^2}{2 \alpha}\)$
For higher dimensions, additional waves add to the instabilities. For a square grid with step size \(h\),
2D: $\( \Delta t \le \frac{h}{4 \alpha}\)\( 3D: \)\( \Delta t \le \frac{h}{6 \alpha}\)$
Once again we have come up against a stiffness issue wherein our timestep is limited by numerical stability rather than accuracy. Unsurprisingly and unfortunatley, this is true of all our explicit time stepping schemes.
Implicit Euler timestepping#
As with the ODE case, we can attempt an Implicit Euler timestepping scheme: $\( \frac{\partial T}{\partial t} = \frac{T^{t+1} -T^{t}}{\Delta t}\)$
and now discretize space at the new timestep:
The stencil is:
A quick glance at the error analysis will see that it is now \(\frac{1}{G}\) which is in the range \(\bigg[1-4\frac{\alpha \Delta t}{\Delta x^2}, 1\bigg]\) and which always satisfies the stability condition.
Once again the Forward Euler method, and implicit method more generally, are unconditional stability but at the cost of solving a system of equations.
Example of Backward Euler time stepping#
# prompt: Repeat solving hte heat transport equation but thistime with a backward Euler method
import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interact, FloatSlider
def backward_euler_heat(T, alpha, dx, dt, num_steps):
"""
Solves the 1D heat equation using backward Euler timestepping.
Args:
T: Initial temperature distribution (numpy array).
alpha: Thermal diffusivity.
dx: Spatial step size.
dt: Time step size.
num_steps: Number of time steps.
Returns:
A list of temperature distributions at each time step.
"""
T_history = [T.copy()]
for _ in range(num_steps):
A = np.diag(-alpha * dt / dx**2 * np.ones(len(T) - 1), -1) + \
np.diag(1 + 2 * alpha * dt / dx**2 * np.ones(len(T)), 0) + \
np.diag(-alpha * dt / dx**2 * np.ones(len(T) - 1), 1)
# Apply boundary conditions (Dirichlet in this example, T=0 at edges)
b = T.copy()
b[1:-1] = T[1:-1]
# plt.spy(A, markersize=1) # Creates a sparsity plot
# plt.title('Sparsity Pattern of Matrix A')
# plt.show()
# Solve the linear system for T at the next time step
T_new = np.linalg.solve(A, b)
# Update T and store the solution
T = T_new.copy()
T_history.append(T.copy())
return T_history
def plot_heat_equation(dt):
"""
Plots the solution of the heat equation for a given time step size.
Args:
dt: Time step size.
"""
# Parameters
alpha = 0.1 # Thermal diffusivity
dx = 0.1 # Spatial step size
num_steps = 50
x = np.arange(0, 1, 1/num_steps)
T_initial = np.zeros_like(x)
T_initial[int(len(x) / 2)] = 1 # Initial heat source in the middle
# Solve the heat equation
T_history = backward_euler_heat(T_initial, alpha, dx, dt, num_steps)
# Plotting
plt.figure(figsize=(8, 6))
for i in range(0, len(T_history), 5):
plt.plot(x, T_history[i], label=f"Time step {i}")
plt.xlabel("Position (x)")
plt.ylabel("Temperature (T)")
plt.title("Backward Euler Solution of 1D Heat Equation")
plt.legend()
plt.grid(True)
plt.show()
# Interactive widget
interact(plot_heat_equation, dt=FloatSlider(min=0.01, max=0.1, step=0.01, value=0.05));
Crank-Nicholson method#
One shortcoming of the backward Euler method is that it is only 1st order accurate in time whereas the central difference spatial discretization is second order. This means that one would need smaller time steps compared to spatial steps to ensure accuracy is optimized.
The Crank-Nicholson method is an alternative implicit method that is second order in both time and space. As with the central difference derivation (and many others!) the key is to average the explicit and implicit methods!
Take, $\( \frac{\partial T}{\partial t} = \frac{T^{t+1} -T^{t}}{\Delta t}\)$
and now average the spatial distretizations at \(t+1\) and \(t\):
Such that the increment scheme is: $\( (1 + 2r) T_i^{t+1} - r T_{i+1}^{t+1} - r T_{i-1}^{t+1} = (1 - 2r) T_i^{t} + r T_{i+1}^{t} + r T_{i-1}^{t}\)$
which seems ugly at first but is actually the exact same matrix sparsity as backward Euler, so achieves better accuracy for comparable computational cost.
The stencil is,
Example of Crank Nicholson accuracy vs. Euler#
# prompt: compare the backward euler and crank-nicholson method for 5 time steps using the previous example. Have the slider no change the time index instead of hte step.
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from ipywidgets import interact, IntSlider
# ... (Your existing code for forward_euler_heat, backward_euler_heat, etc.) ...
def crank_nicholson_heat(T, alpha, dx, dt, num_steps):
"""
Solves the 1D heat equation using the Crank-Nicholson method.
Args:
T: Initial temperature distribution (numpy array).
alpha: Thermal diffusivity.
dx: Spatial step size.
dt: Time step size.
num_steps: Number of time steps.
Returns:
A list of temperature distributions at each time step.
"""
T_history = [T.copy()]
r = alpha * dt / (2 * dx**2)
for _ in range(num_steps):
A = np.diag(1 + 2 * r * np.ones(len(T)), 0) + \
np.diag(-r * np.ones(len(T) - 1), -1) + \
np.diag(-r * np.ones(len(T) - 1), 1)
b = np.diag(1 - 2 * r * np.ones(len(T)), 0) @ T + \
np.diag(r * np.ones(len(T) - 1), -1) @ T + \
np.diag(r * np.ones(len(T) - 1), 1) @ T
plt.spy(A, markersize=1) # Creates a sparsity plot
plt.title('Sparsity Pattern of Matrix A')
plt.show()
# Apply boundary conditions (Dirichlet in this example, T=0 at edges)
b[0] = 0
b[-1] = 0
T_new = np.linalg.solve(A, b)
T = T_new.copy()
T_history.append(T.copy())
return T_history
def plot_comparison(time_index):
"""
Plots the solution of the heat equation using different methods for 5 timesteps.
Args:
time_index: Index of the time step to plot.
"""
# Parameters
alpha = 0.1
dx = 0.1
dt = 0.01
num_steps = 5
x = np.arange(0, 1, 1/10)
T_initial = np.zeros_like(x)
T_initial[int(len(x) / 2)] = 1
# Solve the heat equation using different methods
forward_euler_solution = forward_euler_heat(T_initial, alpha, dx, dt, num_steps)
backward_euler_solution = backward_euler_heat(T_initial, alpha, dx, dt, num_steps)
crank_nicholson_solution = crank_nicholson_heat(T_initial, alpha, dx, dt, num_steps)
# Plotting
plt.figure(figsize=(8, 6))
plt.plot(x, forward_euler_solution[time_index], label="Forward Euler")
plt.plot(x, backward_euler_solution[time_index], label="Backward Euler")
plt.plot(x, crank_nicholson_solution[time_index], label="Crank-Nicholson")
plt.xlabel("Position (x)")
plt.ylabel("Temperature (T)")
plt.ylim(-.1, 1.1)
plt.title(f"Comparison of Methods at Time Step {time_index}")
plt.legend()
plt.grid(True)
plt.show()
# Interactive widget
interact(plot_comparison, time_index=IntSlider(min=0, max=4, step=1, value=0));
#Summary
Once again, we have seen the value of implicit timesteppers for solving stiff problems. However, as discussed in the Initial Value Problem section, they require a root finder to converge on a system, which generally relies on linear solvers. For large systems, direct solution (LU) is intractable, and so iterative methods are required.
Explicit timesteppers are conditionally stable which limits the timestep, but they don’t require simultaneous solutions of equations which avoids the pitfalls of a root finder:
no question of initial guesses
always converges
parallelizes well on HPCs
much smaller memory usage
Does not need an iterative linear solver routine
At what point does the expense of interations in an implicit time stepper outweigh the moderate step size? -> It depends on the problem!
General multiphysics solvers rely on implicit timestepping since modern systems can use LU to root find. Specialized software (build for large systems) may prefer explicit methods which just work:
LS-Dyna - impact tests, crashes, drop tests, etc.
Several CFD packages (Flow 3D)