Parabolic finite-difference methods

Let $E_x(v,w)$ and $E_y(v,w)$ denote the horizontal and vertical summands of $E_h^1(v,w)$ from [F1]. The horizontal telescoping calculation and its vertical counterpart in [F2] give the corresponding unweighted directional identities. Multiplying them by ${\kappa }_x$ and ${\kappa }_y$ yields $⟨A_{h,x}v,w{⟩}_h={\kappa }_x{E}_x(v,w)$ and $⟨A_{h,y}v,w{⟩}_h={\kappa }_y{E}_y(v,w)$, since all boundary values vanish.

By [F3] and step 1.1, $⟨A_{h}v,w{⟩}_h={\kappa }_x{E}_x(v,w)+{\kappa }_y{E}_y(v,w)$. Each edge sum is symmetric in $v,w$, so $⟨A_{h}v,w{⟩}_h=⟨v,A_{h}w{⟩}_h$.

Taking $w=v$ in step 2.1 gives $⟨A_{h}v,v{⟩}_h={\kappa }_x{E}_x(v,v)+{\kappa }_y{E}_y(v,v)$. Every summand is a positive mesh weight times a square, and ${\kappa }_x,{\kappa }_y>0$, hence this quantity is nonnegative.

Changing a representative $i_j$ by $N_j$ changes the exponent by the integer multiple $2\pi \mathrm{i}{\mathrm{\ell }}_j$, so [F2] makes ${\varphi }_{\mathrm{\ell }}$ well-defined on the periodic grid. It is nonzero because ${\varphi }_{\mathrm{\ell }}(0)=1$. Put ${\beta }_j=2\pi {\mathrm{\ell }}_j\mathrm{/}{N}_j$. The addition law gives ${\varphi }_{\mathrm{\ell }}(i\pm e_j)=e^{\pm \mathrm{i}{\beta }_j}{\varphi }_{\mathrm{\ell }}(i)$, including wrapped indices. By [F2], $e^{\mathrm{i}{\beta }_j}$ has modulus one and $e^{\mathrm{i}{\beta }_j}{e}^{-\mathrm{i}{\beta }_j}=1$, hence $e^{-\mathrm{i}{\beta }_j}=\overline{e^{\mathrm{i}{\beta }_j}}$. The Cartesian Euler formula in [F2] and [F3] therefore give $2-e^{\mathrm{i}{\beta }_j}-e^{-\mathrm{i}{\beta }_j}=2-2\cos {\beta }_j=4{\sin }^2({\beta }_j\mathrm{/}2)$. Substitution into [F1] yields the asserted eigenvalue on summing over $j$; it is nonnegative because each ${\kappa }_j\mathrm{/}{h}_j^2$ is positive and each sine square is nonnegative.

Put $P={\prod }_{j=1}^d{N}_j>0$. Repeated complex distributivity $a(b+c)=ab+ac$ moves a scalar through a finite sum. Applying this algebra law coordinate by coordinate, together with [F2] and the finite Fubini rule [F5], gives ${\sum }_{i\in G_h^{\mathrm{p}\mathrm{e}\mathrm{r}}}{\varphi }_{\mathrm{\ell }}(i)\overline{{\varphi }_m(i)}={\prod }_{j=1}^d{\sum }_{r=0}^{N_j-1}\exp (2\pi \mathrm{i}({\mathrm{\ell }}_j-m_j)r\mathrm{/}{N}_j)$. By [F4] this is $P$ when $\mathrm{\ell }=m$ and zero otherwise: distinct index tuples have a coordinate whose two standard representatives are not congruent modulo $N_j$. Thus the nonzero modes are orthogonal, also after multiplying the pairing by the positive grid weight ${\prod }_j{h}_j$. Apply the published one-dimensional inversion formula [F6] successively to each of the $d$ coordinates of any grid vector. Finite Fubini [F5] and the same distributivity combine these expansions into a linear combination of the displayed product modes ${\varphi }_{\mathrm{\ell }}$; hence they span the complexified grid space. If a linear combination of them is zero, multiply it by $\overline{{\varphi }_m}$ and sum over the grid. The displayed pairing identity leaves $P$ times its $m$th coefficient, so every coefficient is zero. The modes therefore form a basis.

Substituting $U^n=a_n{\varphi }_{\mathrm{\ell }}$ into the zero-source update gives $(a_{n+1}-a_n)+\mathrm{\Delta }t{\lambda }_{\mathrm{\ell }}((1-\theta )a_n+\theta a_{n+1})=0$. Since $1+\theta \mathrm{\Delta }t{\lambda }_{\mathrm{\ell }}\ge 1$, solving this scalar equation gives $a_{n+1}=g_{\theta }(\mathrm{\ell })a_n$.

By step 1.2 every grid vector is a finite linear combination of these modes, and the update acts on each component by the factor from step 2.1. Thus the stated factors describe the full update.

At each interior node set $r_{ij}^n=u(x_i,y_j,t_n)$ and define ${\tau }^n=(r^{n+1}-r^n)\mathrm{/}\mathrm{\Delta }t+A_h^{B_{\theta }}((1-\theta )r^n+\theta r^{n+1})-((1-\theta )f^n+\theta f^{n+1})$, where $B_{\theta }=(1-\theta )B^n+\theta B^{n+1}$ and $A_h^{B_{\theta }}$ applies the centered stencil using that exact boundary row. Applying the PDE at $t_n$ and $t_{n+1}$ rewrites this as the weighted average of the two spatial stencil defects plus $(r^{n+1}-r^n)\mathrm{/}\mathrm{\Delta }t-((1-\theta )u_t^n+\theta u_t^{n+1})$.

For a one-variable $C^4$ function $v$, Taylor expansion at an interior node gives $(v(x+h)-2v(x)+v(x-h))\mathrm{/}{h}^2-v^{\mathrm{\prime }\mathrm{\prime }}(x)=h^2(v^{(4)}({\xi }_{+})+v^{(4)}({\xi }_{-}))\mathrm{/}24$ for some ${\xi }_{+}\in (x,x+h)$ and ${\xi }_{-}\in (x-h,x)$. Thus the absolute error is at most $h^2\mathrm{\parallel }{v}^{(4)}{\mathrm{\parallel }}_{\mathrm{\infty }}\mathrm{/}12$. The stencil uses the exact boundary samples at the first and last interior rows, so the same Taylor estimate applies there. Applying it in $x$ and $y$ bounds each spatial defect in step 1.1 by $({\kappa }_x{h}_x^2{M}_{x4}+{\kappa }_y{h}_y^2{M}_{y4})\mathrm{/}12$.

Put $g(t)=u_t(x_i,y_j,t)$. The quotient $(u(t_{n+1})-u(t_n))\mathrm{/}\mathrm{\Delta }t$ is the average of $g$ on the time interval. For $\theta \mathrm{\ne }1\mathrm{/}2$, both this average and $(1-\theta )g(t_n)+\theta g(t_{n+1})$ lie between the minimum and maximum of $g$ on that interval; their difference is therefore at most its oscillation, which is at most $\mathrm{\Delta }t{M}_{tt}$. For $\theta =1\mathrm{/}2$, let $L$ be the linear interpolant of $g$ at the two endpoints. The interpolation remainder satisfies $\mathrm{\mid }g(t)-L(t)\mathrm{\mid }\le M_{ttt}(t-t_n)(t_{n+1}-t)\mathrm{/}2$. Integrating and dividing by $\mathrm{\Delta }t$ gives $\mathrm{\mid }(\mathrm{\Delta }t{)}^{-1}{\int }_{t_n}^{t_{n+1}}g-(g(t_n)+g(t_{n+1}))\mathrm{/}2\mathrm{\mid }\le (\mathrm{\Delta }t{)}^2{M}_{ttt}\mathrm{/}12$. These are the required bounds for the temporal term in step 1.1.

Combining the two spatial bounds from step 2.1 with the temporal bound from step 2.2 in the residual identity of step 1.1 gives the two asserted estimates at every interior node. Taking the maximum over nodes and time steps proves the stated complete-residual bounds.

If $\theta =0$, the implicit matrix in [F1] is the identity, so for every right-hand side there is exactly one solution, namely that right-hand side.

Suppose $\theta >0$ and $(I+\theta \mathrm{\Delta }t{A}_h)v=0$. Taking the weighted inner product with $v$ gives $\mathrm{\parallel }v{\mathrm{\parallel }}_h^2+\theta \mathrm{\Delta }t⟨A_{h}v,v{⟩}_h=0$. By [F2] the second term is nonnegative, so $\mathrm{\parallel }v{\mathrm{\parallel }}_h^2=0$ and $v=0$.

The interior vector space is finite dimensional and the matrix is square. Its zero kernel from step 1.2 means its rank equals the dimension by rank-nullity; hence its range is the whole space. Thus every right-hand side has exactly one solution.

Since $0\le {\sin }^2(\pi {\mathrm{\ell }}_j\mathrm{/}{N}_j)\le 1$, [F1] gives $0\le {\lambda }_{\mathrm{\ell }}\le 4S$ for every mode. If $\mathrm{\Delta }tS\le 1\mathrm{/}2$, then $0\le \mathrm{\Delta }t{\lambda }_{\mathrm{\ell }}\le 2$, so $\mathrm{\mid }{g}_0(\mathrm{\ell })\mathrm{\mid }=\mathrm{\mid }1-\mathrm{\Delta }t{\lambda }_{\mathrm{\ell }}\mathrm{\mid }\le 1$.

The Fourier modes in [F1] are an orthogonal basis. Each mode is multiplied by a scalar of magnitude at most one in step 1.1; summing the squared modal coefficients therefore shows that one update is nonexpansive in the weighted grid $L^2$ norm, and every finite power remains nonexpansive.

If every $N_j$ is even, choose ${\mathrm{\ell }}_j=N_j\mathrm{/}2$ for every coordinate. Then each sine square in [F1] equals one, so ${\lambda }_{\mathrm{\ell }}=4S$ and $g_0(\mathrm{\ell })=1-4\mathrm{\Delta }tS$. It equals $-1$ at $\mathrm{\Delta }tS=1\mathrm{/}2$ and is strictly less than $-1$ when $\mathrm{\Delta }tS>1\mathrm{/}2$. Repeated updates multiply this mode by $g_0(\mathrm{\ell }{)}^n$, whose magnitude grows geometrically in the latter case.

By [F3] the step has a unique solution. Write $\delta U=U^{n+1}-U^n$. The identity $U_{\theta }=U^n+\theta \delta U$ gives $⟨\delta U,U_{\theta }{⟩}_h=\frac{1}{2}(\mathrm{\parallel }{U}^{n+1}{\mathrm{\parallel }}_h^2-\mathrm{\parallel }{U}^n{\mathrm{\parallel }}_h^2)+(\theta -\frac{1}{2})\mathrm{\parallel }\delta U{\mathrm{\parallel }}_h^2$ by expanding the two squared norms.

Taking the inner product of [F1] with $U_{\theta }$ and using step 1.1 gives $\frac{1}{2}(\mathrm{\parallel }{U}^{n+1}{\mathrm{\parallel }}_h^2-\mathrm{\parallel }{U}^n{\mathrm{\parallel }}_h^2)+(\theta -\frac{1}{2})\mathrm{\parallel }\delta U{\mathrm{\parallel }}_h^2+\mathrm{\Delta }t⟨A_h{U}_{\theta },U_{\theta }{⟩}_h=\mathrm{\Delta }t⟨q^n,U_{\theta }{⟩}_h$. The two terms after the norm difference on the left are nonnegative by $\theta \ge 1\mathrm{/}2$ and [F2]. Cauchy–Schwarz and the triangle inequality give $\mathrm{\mid }⟨q^n,U_{\theta }{⟩}_h\mathrm{\mid }\le \mathrm{\parallel }{q}^n{\mathrm{\parallel }}_h\mathrm{\parallel }{U}_{\theta }{\mathrm{\parallel }}_h\le \mathrm{\parallel }{q}^n{\mathrm{\parallel }}_h\max (\mathrm{\parallel }{U}^n{\mathrm{\parallel }}_h,\mathrm{\parallel }{U}^{n+1}{\mathrm{\parallel }}_h)$.

Put $a=\mathrm{\parallel }{U}^n{\mathrm{\parallel }}_h$, $b=\mathrm{\parallel }{U}^{n+1}{\mathrm{\parallel }}_h$, and $r=\mathrm{\parallel }{q}^n{\mathrm{\parallel }}_h$. If $b\le a$, the claimed bound is immediate. If $b>a$, step 2.1 gives $(b^2-a^2)\mathrm{/}2\le \mathrm{\Delta }trb$; dividing by $(a+b)\mathrm{/}2>0$ yields $b-a\le 2\mathrm{\Delta }trb\mathrm{/}(a+b)\le 2\mathrm{\Delta }tr$. This proves the forced estimate in all cases. When $q^n=0$, it gives $b\le a$; iteration over any finite number of steps preserves the same mesh-independent constant $1$.

Put $s={\mu }_x+{\mu }_y$ and let $V^k$ denote the full grid values, equal to $U^k$ in the interior and $B^k$ on the boundary. At each interior node the zero-source update is $(1+2\theta s)U_{ij}^{n+1}-\theta {\mu }_x(V_{i+1,j}^{n+1}+V_{i-1,j}^{n+1})-\theta {\mu }_y(V_{i,j+1}^{n+1}+V_{i,j-1}^{n+1})=[1-2(1-\theta )s]U_{ij}^n+(1-\theta ){\mu }_x(V_{i+1,j}^n+V_{i-1,j}^n)+(1-\theta ){\mu }_y(V_{i,j+1}^n+V_{i,j-1}^n)$. The right side is known; its center coefficient is nonnegative by the stated condition, and its neighbor coefficients are nonnegative because $\theta \in [0,1]$ and ${\mu }_x,{\mu }_y>0$. Thus it is nonnegative whenever the old interior and boundary values are nonnegative.

Let $W$ be the unique interior solution of the implicit system with nonnegative right-hand side and nonnegative future boundary data, whose existence is [F2]. If an interior value attained a negative minimum $W_{ij}<0$, then every adjacent interior or boundary value would be at least $W_{ij}$, so $W_{ij}+\theta {\mu }_x(2W_{ij}-W_{i+1,j}-W_{i-1,j})+\theta {\mu }_y(2W_{ij}-W_{i,j+1}-W_{i,j-1})\le W_{ij}<0$. This is the corresponding component of the nonnegative right-hand side, a contradiction. Hence $W\ge 0$; for $\theta =0$ the same inequality reduces to $W_{ij}<0$ and gives the contradiction directly.

The zero-source scheme is linear. Let $M_0$ and $m_0$ be respectively the maximum and minimum of the finite collection of initial and prescribed boundary values at all time levels. The constant grid with interior and boundary values $M_0$ solves the zero-source scheme with constant boundary rows $M_0$, since every centered second difference of a constant is zero. Subtracting the numerical solution $U$ gives a solution with initial values $M_0-U^0\ge 0$ and boundary rows $M_0-B^n\ge 0$. Induction using steps 1.1–1.2 gives $U\le M_0$. Likewise, subtracting the constant grid with interior and boundary values $m_0$ from $U$ gives initial values $U^0-m_0\ge 0$ and boundary rows $B^n-m_0\ge 0$, so induction gives $U\ge m_0$. Thus $m_0\le U_{ij}^n\le M_0$ at every interior node and time level; the prescribed boundary nodes satisfy the same bounds by definition.

In the homogeneous-boundary case at $\theta =1\mathrm{/}2$, the stated sufficient maximum-principle restriction reduces to $\mathrm{\Delta }t({\kappa }_x\mathrm{/}{h}_x^2+{\kappa }_y\mathrm{/}{h}_y^2)\le 1$. By [F3] weighted $L^2$ stability holds for every $\mathrm{\Delta }t>0$ under the same zero-boundary setup. Thus the maximum-principle condition is a stricter sufficient condition; this comparison makes no claim that the maximum principle fails outside it.

Define the interior error $e^n=r^n-U^n$. Exact initial values give $e^0=0$, and exact boundary rows give $e^n=0$ on every boundary node. Subtracting the numerical equation from the exact sampled equation whose residual is [F3] shows that $e^n$ satisfies the same zero-Dirichlet θ update with source $q^n={\tau }^n$.

Applying [F2] at each time step gives $\mathrm{\parallel }{e}^{n+1}{\mathrm{\parallel }}_h\le \mathrm{\parallel }{e}^n{\mathrm{\parallel }}_h+2\mathrm{\Delta }t\mathrm{\parallel }{\tau }^n{\mathrm{\parallel }}_h$. Summing these inequalities from $n=0$ to $m-1$ and using $e^0=0$ yields $\mathrm{\parallel }{e}^m{\mathrm{\parallel }}_h\le 2\mathrm{\Delta }t{\sum }_{n=0}^{m-1}\mathrm{\parallel }{\tau }^n{\mathrm{\parallel }}_h$.

The weighted grid norm has total interior weight $h_x{h}_y(N_x-1)(N_y-1)<L_x{L}_y$, so [F3] gives $\mathrm{\parallel }{\tau }^n{\mathrm{\parallel }}_h\le \sqrt{L_x{L}_y}{R}_{\theta }$. Since $m\mathrm{\Delta }t\le T$, step 2.1 implies $\mathrm{\parallel }{e}^m{\mathrm{\parallel }}_h\le 2T\sqrt{L_x{L}_y}{R}_{\theta }$. Substituting the two explicit values of $R_{\theta }$ from [F3] proves the stated rates.