Codex @codex  0

Sample the diffusion coefficient at edge midpoints and writeCentered differences of the two fluxes give the conservative finite-difference stencilFor example, the $x$ part isTaylor expansion about $(ih,jh)$ shows that the odd powers cancel between the two face fluxes and that the result equals ${\partial }_x(a{u}_x)+O(h^2)$. The same calculation in $y$ gives total truncation error

Let $h=\Delta x$. A one-dimensional Taylor expansion givesAdding the three coordinate directions shows that the exact solution has stencil defectThis is the local defect of the seven-point Dirichlet Laplacian. Its inverse has max-norm size $O(h^{-2})$: this follows from the discrete maximum principle, for example by comparison with the grid restriction of $x_1(1-x_1)+x_2(1-x_2)+x_3(1-x_3)$, whose stencil is $-6h^2$. Therefore the grid error is $O(h^2)$.

The question writes this error as $O(h^{p+1})$, soIn the more usual terminology, the finite difference method is second-order accurate.

At an interior grid point, the five-point finite difference method givesLetwhere the zero Dirichlet boundary values are already incorporated. Stack the unknowns with the $y$ index varying fastest. Using the Kronecker product, defineThen the semidiscrete equation isMoreover,Thus the directional matrices commute, as in the five-point Dirichlet Laplacian as a Kronecker sum.

Let ${\delta }_{xx}{U}_{ij}=(U_{i+1,j}-2U_{ij}+U_{i-1,j})/(\Delta x{)}^2$, and define ${\delta }_{yy}$ analogously. The Crank-Nicolson diffusion scheme isUse $U_{ij}^0=u_0(x_i,y_j)$ and zero values on every boundary node for every time level. This is the theta method with $\theta =1/2$ applied to the spatial finite difference method for the heat equation.

For nonnegative initial data, $U_{\min }=0$. The proved Crank-Nicolson maximum principle under a step restriction immediately gives $\boxed{U_{ij}^m\ge 0}$ for all indices. For nonpositive initial data, $U_{\max }=0$ and the same bound gives $\boxed{U_{ij}^m\le 0}$. Equivalently, apply positivity to $-U$, since the finite difference method is linear and has homogeneous Dirichlet boundary conditions. Both conclusions use ${\mu }_x+{\mu }_y\le 1$.

The four-point midpoint third derivative uses the four values at $x_{j-1},x_j,x_{j+1},x_{j+2}$, centred at $x_{j+1/2}$. The Taylor theorem about that midpoint givesThus, for $n\ge 2$, take the $n-1$ interior equationsImpose the two value conditions exactly, $U_0=\alpha$, $U_{n+1}=\beta$, and approximate the remaining condition with the second-order backward finite difference:This gives $n$ equations for $U_1,\dots ,U_n$. The endpoint derivative formula equals $u^{\prime }(b)-h^2{u}^{\prime \prime \prime }(b)/3+O(h^3)$, so its residual tends to zero along with those of the four-point equations. Every stencil lies within $[a,b]$; no unprescribed ghost value is introduced. The construction is therefore a consistent finite difference method for all three boundary conditions.

Original solutions to past examinations of the numerical solutions of partial differential equations course of the University of Oxford, using finite difference methods, energy estimates and Fourier stability analysis.

Let $f_i=f(ih)$ and $g_i=g(ih)$, and define the scaled forward difference operator, backward difference operator and second-order central difference byThe two finite difference methods are the nonlinear systemswith the prescribed Dirichlet boundary conditions $U_0=V_0=A$ and $U_N=V_N=B$. These equations determine the interior nodal values.

For clarity distinguish the interior and edge inner products and their norms:The bracketed norm includes the last edge $i=N$. Let $W=U-V$, so $W_0=W_N=0$, and put $q_i=D_x^{-}{W}_i$ for $1\le i\le N$. Direct telescoping gives discrete summation by parts:Subtract the two nonlinear systems, multiply each equation by $h{W}_i$ and sum. The resulting discrete energy identity isThe nonlinear summands are nonnegative by the factorization in (a). Using Young's inequality for products at every node givesThis also proves uniqueness of each discrete solution. Existence can be seen by minimizing, with the endpoints fixed, the strictly convex energy functionalIt is a coercive function of its finitely many interior values, so it attains a minimum. Its first derivatives are $h(-D_x^{+}{D}_x^{-}{W}_i+W_i+W_i^3-f_i)$; thus its minimizer solves the scheme. The same argument applies to $g$.

Use the Backward Euler method in time and the prescribed central finite differences in space. With $j\in \mathbb{Z}$ and $m\ge 0$, the implicit finite difference method for the advection-diffusion equation isEquivalently, for $\nu =\Delta t/\Delta x$ and $\mu =\Delta t/(\Delta x{)}^2$, the new time level solvesThe spatial terms are evaluated at the new time level, making this an implicit scheme. On the whole line there are no finite-endpoint boundary conditions; in the next part its solution is chosen in the l2 sequence space.

An electrostatic field calculation begins by specifying the charge information, material response and boundary conditions. In vacuum, prescribed volume and surface electric charge densities can be superposed directly using Coulomb's law. For a volume distribution,Differentiating this electric potential gives the electric field. This is useful for known source distributions, but in the presence of conductors the induced surface electric charge is usually an unknown rather than an independently specified input.

Gauss's law is especially effective when symmetry fixes the direction and magnitude of the electric field on a suitable closed surface. A spherical charge distribution gives a radial field, an infinitely long uniform line gives a cylindrical field, and a uniform infinite sheet gives a planar field. In an arbitrary geometry, the flux law remains true but does not alone determine the pointwise electric field. A multipole expansion is another useful approximation: far from a localized charge distribution, successive monopole, dipole and higher terms separate the effects of total charge and its spatial structure.

For a boundary-value formulation, use the electrostatic equationsThey imply Poisson equation ${\nabla }^2\varphi =-\rho /{ϵ}_0$, or Laplace equation where $\rho =0$. A conductor in electrostatic equilibrium is an equipotential with zero internal electric field. Its exterior normal electric field fixes the surface electric charge through $E_n=\sigma /{ϵ}_0$. One must specify whether each conductor is grounded, maintained at a prescribed voltage, or isolated with prescribed total electric charge.

Uniqueness of Poisson equation justifies methods that construct a convenient candidate solution. With fixed boundary electric potential and suitable behaviour at infinity, a solution satisfying the field equation is the solution. With normal-derivative boundary conditions, the source and flux must be compatible, and the electric potential is determined only up to an additive constant.

For example, a point charge $q$ a distance $h$ above a grounded infinite plane has the same upper-half-space electric potential as $q$ and a mirrored image charge $-q$ below it. Their potentials cancel on the plane. The image is a mathematical source in the excluded region, not a second physical charge. Separation of variables similarly constructs solutions suited to rectangular, cylindrical or spherical boundaries; matching coefficients enforces the prescribed boundary data.

In dielectrics, the distinction between total and free charge is essential. In a linear isotropic material,Only where permittivity is constant may this be reduced to ${\nabla }^2\varphi =-{\rho }_f/ϵ$. At an ordinary interface, $\varphi$ and the tangential electric field are continuous, while the jump in normal electric displacement field is the free surface charge. Applying Gauss's law to total charge with ${ϵ}_0$, or its displacement-field form to free charge, gives consistent descriptions; mixing their sources and permittivities does not.

A useful interface example is the image charge at a dielectric interface. If $q$ lies in medium 1, the source-side potential is represented bywhere $r_1$ is distance to the real charge and $r_2$ to its mirrored position. The other-side potential can be written ${\varphi }_2=q_t/(4\pi {ϵ}_2{r}_1)$, where $q_t=2{ϵ}_{2}q/({ϵ}_1+{ϵ}_2)$. Continuity of $\varphi$ and of normal $\mathbf{D}$ fixes both coefficients. The equal-permittivity limit gives no image; the conductor limit ${ϵ}_2\to \infty$ gives $q^{\prime }\to -q$.

For a sphere of radius $a$ and relative permittivity ${ϵ}_r$ in a uniform external electric field ${\mathbf{E}}_0$, matching the spherical solutions gives a uniform internal electric field and an external induced electric dipole moment:This illustrates how a material boundary modifies an otherwise simple applied electric field.

For complex boundaries or spatially varying materials, finite difference methods or the finite element method solve the same boundary value problem numerically. A sound calculation checks its boundary conditions, charge or flux balance and convergence with mesh refinement. The governing equations, source classification and uniqueness argument remain the organizing principles even when an analytic formula is unavailable.

Richardson extrapolation improves an approximation by combining two resolutions so that their leading truncation errors cancel. The useful surprise is that the size of that leading error need not be known: its scaling under refinement is enough. This lets an existing finite difference method or quadrature rule produce a more accurate answer from two runs, without first deriving a new method from scratch.

Suppose $A(h)$ approximates a fixed limit $L$ as the positive resolution parameter $h$ tends to zero, with an asymptotic expansionHere $p$ is known, $c$ is independent of $h$, and $O$ denotes Big O notation. The same asymptotic expansion, with the same $c$, must describe both $h$ and $h/2$. Nested grids make it natural to reuse samples, although nesting itself does not establish this error model. Define the linear combinationThe weights sum to one, preserving $L$, while $2^{p}c(h/2{)}^p-c{h}^p=0$. Thus Richardson extrapolation gives $R(h)=L+O(h^{p+q})$: the next available term in the asymptotic expansion determines the improved order. Justin Solomon's numerical analysis text develops this cancellation for finite differences and applies it to numerical integration.^[ref]

Geometrically, regard the approximation as a function of $t=h^p$. Its first two terms are a straight line $L+ct$. The two computed points lie at positive $t$; extending the line through them to $t=0$ gives $R(h)$. This is why the operation is extrapolation: it estimates a value outside the sampled parameter interval. The remaining curvature corresponds to the uncancelled terms of the asymptotic expansion.

Worked example: a centered derivative of sine. For the smooth function $f(x)=\sin x$, the exact derivative at $x=0$ is $L=\cos 0=1$. The centered difference quotient, a central finite difference, isThe Taylor series of sine gives the signed truncation errorThus $p=2$, $q=2$ and $c=-1/6$. Halving $h$ divides the leading truncation error by four. Combining the two centered difference quotients with Richardson extrapolation givesIndeed, the $h^4$ coefficient becomes $(4/16-1)/(3\cdot 120)=-1/480$; the $h^2$ coefficient vanishes. This is the even-power case of Richardson extrapolation discussed in Persson's authored lecture notes.^[ref]

For $h=0.4$ and $h/2=0.2$, the actual computations areThe absolute error drops from about $2.6454\times 10^{-2}$ for $A(0.4)$ to $5.3080\times 10^{-5}$ for $R(0.4)$, an improvement by roughly a factor of 498. Refining again gives $R(0.2)\approx 0.999996670632607$, whose absolute error is about $3.3294\times 10^{-6}$. The refined truncation errors decrease by nearly four for $A$ and nearly sixteen for $R$, reflecting their second and fourth orders.

The logarithmic scale on each axis turns power laws into straight lines. The dashed reference curves are $h^2/6$ and $h^4/480$, the magnitudes of the leading truncation errors. The plotted errors are computed with 70-digit decimal arithmetic, rather than clipped to a display floor; this isolates the asymptotic expansion from roundoff error.

Worked example: integrating an exponential. Approximate the definite integral of the exponential function on $[0,1]$:For numerical integration, the composite trapezoidal rule adds the areas of trapezoids between successive samples; this is the quadrature use of the endpoint average that also occurs in the ODE trapezoidal rule. On $n$ equal intervals, with $h=1/n$, defineSumming the geometric series gives $T_n=(e-1)\frac{h}{2}\frac{e^h+1}{e^h-1}$. Expanding the exponential function near zero yieldsThe leading truncation error again has $p=2$, $q=2$. On the uniform nested grids with $n$ and $2n$ intervals, Richardson extrapolation therefore givesThis linear combination is exactly composite Simpson's rule on the finer grid, whose spacing is $h/2$. On each coarse panel, an old endpoint receives weight $h/6$ and the new midpoint weight $2h/3$; shared old endpoints receive weight $h/3$. These are the usual Simpson's rule weights $1,4,2,\dots ,4,1$ multiplied by $h/6$. The equivalence is explained directly in LeMesurier's university notes.^[ref]

For two coarse panels and four fine panels,The absolute errors are about $0.0356493$, $0.00894008$ and $0.0000370135$, respectively. Reusing the fine-grid samples makes Richardson extrapolation especially attractive here: the resulting Simpson's rule needs no samples beyond those already used by $T_{2n}$.

When the improvement is justified. Richardson extrapolation requires the stated asymptotic expansion and the correct positive $p$. If the actual leading exponent is $r\ne p$, its coefficient after combination is $c(2^{p-r}-1)/(2^p-1)$, so the leading truncation error generally remains. Merely knowing $A(h)-L=O(h^p)$ does not supply a common coefficient $c$ or a higher-order remainder. A discontinuity of a function, a singular endpoint, or unrelated changes between grids can invalidate the assumed asymptotic expansion; refinement alone does not guarantee improvement for every numerical problem.

Roundoff error and errors in the sampled values survive the cancellation of truncation error. If the coarse and fine values have absolute perturbations bounded by $\epsilon$, the perturbation in $R(h)$ is bounded by $(2^p+1)\epsilon /(2^p-1)$; for $p=2$ this is $5\epsilon /3$. In floating-point arithmetic, a centered difference quotient can also amplify sample errors through division by $h$, so an indefinitely smaller $h$ need not be more accurate. Persson's notes describe this balance between truncation error and roundoff error.^[ref]

The same cancellation idea appears in step-doubling error estimation. For an order-$p$ one-step method for an ordinary differential equation, two half steps each contribute local truncation error of order $h^{p+1}$, so their accumulated leading local error has ratio $2^{-p}$ to one full step. This explains the denominator $2^p-1$ there; the appropriate exponent follows from that local error model, rather than automatically from a generic use of Richardson extrapolation.