dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation

Class CrankNicolsonHeatEquation1D

java.lang.Object

dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.CrankNicolsonHeatEquation1D

All Implemented Interfaces:

PDESolver

public class CrankNicolsonHeatEquation1D
extends Object
implements PDESolver

The Crank-Nicolson method is an algorithm for obtaining a numerical solution to parabolic PDE problems. It is based on central difference in space, and the trapezoidal rule in time, giving second-order convergence in time. That is, the solution is computed using this implicit equation: \[ \frac{u_{j}^{k + 1} - u_{j}^{k}}{\Delta t} = \frac{\beta}{2}\left( \frac{u_{j+1}^{k+1} - 2u_{j}^{k+1} + u_{j-1}^{k+1}}{\Delta x^2} + \frac{u_{j+1}^{k} - 2u_{j}^{k} + u_{j-1}^{k}}{\Delta x^2} \right) \] where \(u^k_j = u(t_k,x_j)\) is the estimate at \((k, j)\) in the grid.

This is a second-order method with a truncation error of order \(O(\Delta t^2 + \Delta x^2)\).

The solution is computed along a two-dimensional grid (e.g. time and space). Note that due to the way that the boundary values are given, the resulting grid will NOT contain the values at the boundaries and attempting to access \(u^k_0\) or \(u^k_{n+1}\) will result in an UnsupportedOperationException being thrown.

Grid sizes and interpretation

In solve(pde,M,N), M is the number of time intervals after time zero, and N is the number of interior spatial points. Thus \(\Delta t=T/M\) and \(\Delta x=a/(N+1)\). Supply M at least one and N at least two; the current stencil accesses neighboring interior values.

Read time levels 0 through M and interior spatial indices 1 through N. Coordinates x(0) and x(N+1) are available, but u(m,0) and u(m,N+1) throw UnsupportedOperationException. Query the supplied boundary functions separately when their values are needed. This implicit method does not make a coarse grid accurate: refine time and space and compare an analytic solution or a convergence study. Nonsmooth data can still cause oscillatory values.

Heat diffusion with fixed endpoint temperatures

The analytic sine-mode solution gives an independent accuracy check. The example also distinguishes boundary coordinates from unavailable boundary solution values.

 import dev.nm.analysis.function.rn2r1.univariate.AbstractUnivariateRealFunction;
 import dev.nm.analysis.differentialequation.pde.finitedifference.PDESolutionTimeSpaceGrid1D;
 import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.HeatEquation1D;
 import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.CrankNicolsonHeatEquation1D;
 AbstractUnivariateRealFunction initial = new AbstractUnivariateRealFunction() {
     public double evaluate(double x) { return Math.sin(Math.PI * x); }
 };
 AbstractUnivariateRealFunction zero = new AbstractUnivariateRealFunction() {
     public double evaluate(double t) { return 0; }
 };
 HeatEquation1D problem = new HeatEquation1D(0.1, 1, 0.1, initial, 0, zero, 0, zero);
 PDESolutionTimeSpaceGrid1D grid = new CrankNicolsonHeatEquation1D().solve(problem, 40, 39);
 assert grid.M() == 40 && grid.N() == 39;
 assert Math.abs(grid.t(40) - 0.1) < 1e-14;
 assert Math.abs(grid.x(20) - 0.5) < 1e-14;
 double exact = Math.exp(-0.1 * Math.PI * Math.PI * grid.t(40)) * Math.sin(Math.PI * grid.x(20));
 assert Math.abs(grid.u(40, 20) - exact) < 1e-4;
 assert grid.x(0) == 0.0;
 boolean boundaryUnavailable = false;
 try { grid.u(40, 0); } catch (UnsupportedOperationException expected) { boundaryUnavailable = true; }
 assert boundaryUnavailable;
 

Insulated endpoints

Zero flux at both ends preserves a constant profile. Neumann data specify derivatives through the signed boundary convention, not temperatures.

 import dev.nm.analysis.function.rn2r1.univariate.AbstractUnivariateRealFunction;
 import dev.nm.analysis.differentialequation.pde.finitedifference.PDESolutionTimeSpaceGrid1D;
 import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.HeatEquation1D;
 import dev.nm.analysis.differentialequation.pde.finitedifference.parabolic.dim1.heatequation.CrankNicolsonHeatEquation1D;
 AbstractUnivariateRealFunction initial = new AbstractUnivariateRealFunction() {
     public double evaluate(double x) { return 7; }
 };
 AbstractUnivariateRealFunction noFlux = new AbstractUnivariateRealFunction() {
     public double evaluate(double t) { return 0; }
 };
 HeatEquation1D problem = new HeatEquation1D(0.1, 1, 0.2, initial, 1, noFlux, 1, noFlux);
 PDESolutionTimeSpaceGrid1D grid = new CrankNicolsonHeatEquation1D().solve(problem, 4, 9);
 for (int n = 1; n <= grid.N(); ++n) {
     assert Math.abs(grid.u(grid.M(), n) - 7) < 1e-10;
 }
 

See Also:

Wikipedia: Crank-Nicolson method

Nested Class Summary

Nested Classes

Modifier and Type	Class and Description
static class	CrankNicolsonHeatEquation1D.Coefficients Gets the coefficients of a discretized 1D heat equation for each time step.

Constructor Summary

Constructors

Constructor and Description
CrankNicolsonHeatEquation1D()

Method Summary

Modifier and Type	Method and Description
PDESolutionTimeSpaceGrid1D	solve(HeatEquation1D pde, int M, int N) Solves the given one-dimensional heat equation.

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Constructor Detail

CrankNicolsonHeatEquation1D

public CrankNicolsonHeatEquation1D()

Method Detail

solve

public PDESolutionTimeSpaceGrid1D solve(HeatEquation1D pde,
                                        int M,
                                        int N)

Solves the given one-dimensional heat equation.

Parameters:

pde - the PDE problem

M - the number of grid points along the time-axis (excluding the initial condition)

N - the number of grid points along the space-axis (excluding the two boundaries)

Returns:

the solution \(u(t,x)\) at the grid points