Reduction of Parabolic Equations to a System of ODE's

Consider $\left\{\begin{array}{ll} \displaystyle \frac{\partial U}{\partial t}=\frac{\pa... ...and known BV's at }x=0 & \mbox { and } x= X\quad \forall t>0 \end{array}\right.$

Semi-discretizing, using center differences in $x$ (as a particular example)

$$\frac{dU(t)}{dt}=\frac{1}{h^2}\Big\{U(x-h, t)-2U(x,t)+U(x+h,t)\Big\}+O(h^2)$$

Subdivide interval $0\le x\le X$ into $N$ equal subintervals with $x_i=ih$,      $i=0\cdots N$ ,     where $Nh=X$. $u_i(t)$ is an approximation to $u_i(t)$, where $U_i=u(ih,t)$ so the $i=0$ and $i=N$ are boundary lines. The equation for $u_i$ at some time $t$ are

$$N-1$$$$\ \mbox { ODE's} \left\{\begin{array}{l} \frac{du_1(t)}{dt}=\fr... ... \frac{du_{N-1}}{dt}=\frac{1}{h^2}(u_{N-2}-2u_{N-1}+u_N) \end{array}\right.$$

$u_0$ and $u_N$ are known values (boundary values) let ${\bf V}(t)=[u_1, u_2\ldots u_{N-1}]^T$ then the system can be written as

$$\frac{d{\bf V}(t)}{dt}=A {\bf V}(t) +{\bf b}$$ (144)

$${\bf b}$$$$\mbox { is a column of zero's and known values of } u$$

$$A=\displaystyle \frac{1}{h^2}\left[ \begin{array}{cccccc} -2... ... & 0 & 0 & 1 & -2 \\ \end{array}\right]\quad (N-1)\times (N-1)$$ matrix

Remark: for $y=y(t)$, the solution to

$$\frac{dy}{dt}= ay+b$$

$$y(0)=g$$

$$\mbox {is } \quad y(t)=-\frac{b}{a}+\left(g+\frac{b}{a}\right)e^{at}$$

$\Box$

Solution to

$${\bf V}(t)=-A^{-1}{\bf b}+e^{tA}(\bf {g+A^{-1} b})$$ (145)

$$\therefore \quad {\bf V}(t+k)=-A^{-1}{\bf b}+e^{kA}e^{tA}(\bf g+A^{-1}{\bf b})$$ (146)

substitutions (146) in (147)

$${\bf V}(t+k)=-A^{-1}{\bf b}+e^{kA}\left({\bf V}(t)+A^{-1}{\bf b}\right)$$

Note: if ${\bf b}=0\Rightarrow {\bf V}(t+k)=e^{kA}{\bf V}(t)$.

Stability: perturb $g$ to $g^{*}$ then

$${\bf V}^{*}(t)=A^{-1}{\bf b}+e^{tA}({\bf g}^{*}+A^{-1}{\bf b})$$

subtracting $\underbrace{{\bf V}^{*}(t)-{\bf V}(t)}_{\bf e(t)}=e^{TA}\underbrace{ ({\bf g}^{*}-{\bf g})} _{\bf e(0)}$

$${\bf e}(t)=e^{tA}{\bf e}(0)$$

$$\mbox { so we require that }\vert\vert A\vert\vert\le 1 \mbox { for stability}.$$

Note on $$\frac{dV}{dt}=A{\bf V}+{\bf b}$$

Take $P$ a constant coefficient real $n\times n$ matrix, then

$$e^{P}=I+P+\frac{P^2}{2}+\frac{P^3}{3!}+\cdots \sum^ {\infty}_{m=0}\frac{P^m}{m!}$$ (147)

$$\mbox {here } P^0\equiv I \mbox { is the $n \times n$ identity matrix}$$

If $Q$ is a real $n\times n$ matrix such that $PQ=QP$ (commute) then

$$e^Pe^Q=e^Qe^P=e^{P+Q}$$

Hence $e^Pe^{-P}=e^{-P}e^P=e^0=I$

premultiplication of $e^Pe^{-P}=I$ by $(e^P)^{-1}$ then shows that

$$e^{-P}=(e^P)^{-1}$$

On putting $P=At$ in (148) and differentiating with regards to $t$ we get that

$$\frac{d}{dt}(e^{At})=Ae^{At}=e^{At}A$$

Now consider ${\bf V}(t)=e^{At}{\bf g}$ where ${\bf g}$ is independent of $t$. This clearly satisfies the condition ${\bf V}(0)={\bf g}$. Differentiation with regards to $t$ gives

$$\frac{d{\bf V}}{dt}=Ae^{At}g=A{\bf V}$$

In other words the solution of

$$\frac{d{\bf V}}{dt}=A{\bf V} \quad \mbox { with } {\bf V}(0)={\bf g}\quad , \quad \mbox {is}$$

$${\bf V}(t)=e^{At}{\bf g}$$

Similarly, the vector function $\left\{\begin{array}{ll} {\bf V}(t)=-A^{-1}{\bf b} +e^{tA}({\bf g}+A^{-1}{\bf b})\\ {\bf V}(0)={\bf g} \end{array}\right.$

is the solution of $$\frac{d{\bf V}}{dt}=A{\bf V}+{\bf b}$$

provided $b$ and $A$ are independent of $t$

Finite Difference Schemes from Systems of ODE's

For simplicity assume $g$ given and that the boundary values associated with $$\frac{\partial u}{\partial t}=\frac{\partial ^2 u}{\partial x^2}$$ are 0.

$$\mbox { for }0\le x\le X$$

$$\mbox { The solution is approximately}$$

$${\bf V}(t+k)=e^{k A}{\bf V}(t)\quad t=0, k, 2k, \ldots$$ (148)

where $A$ is given before. The FD comes in approximating $e^{kA}$. First, notice that

$$e^{kA}=I+kA+\frac12k^2A^2+\cdots$$

Then, an obvious approximation is $e^{kA}\approx I+kA$, so (149) is approximately

$${\bf V}(t+k)=(I+kA)V(t)$$ (149)

if $t=nk$ and $\mu=k/h^2$ then (150) is

$\left[\begin{array}{l} u^{n+1}_1\\ u^{n+1}_2\\ :\\ :\\ u^{n+1}_{M-1} \end{... ...ts & \ddots \\ & \mu & (1-2\mu) &\mu\\ &0 & \mu &(1-2\mu) \end{array}\right]$ $\left[\begin{array}{l} u^n_1\\ u_2^n\\ :\\ :\\ u^{n}_{M-1} \end{array}\right]$

or $u^{n+1}_m=\mu u^n_{m-1}+(1-2\mu)u^n_m+\mu u^2_{m+1},\quad m=1,2\ldots M-1$

   with $$Mh=X$$

We can use Padé Approximants to get better approximations to $e^{kA}$.

Padé Approximants to $e^{\theta}\;,\;$ where $\theta$ is real:

Assume $e^{\theta}$ can be approximated as $$\frac{1+p_1\theta}{1+q_1\theta}\quad ; p_1, q_1$$ constants.

then we need 2 equations to determine $p_1, q_1 .$

$$e^{\theta}=\frac{1+p_1\theta}{1+q_1\theta}+c_3\theta^3+c_4\theta^4 +\cdots$$

multiplying both sides by denominator

$$\therefore (1+q_1\theta)(1+\theta+\f... ...ac16\theta^3\cdots)=1+p_1\theta+(1+q_1 \theta)(c_3\theta^3+c_4\theta^4+\cdots)$$

Hence $(1+q_1-p_1)\theta +(\frac12+q_1)\theta^2+\left(\frac16+\frac12q_1-c_3 \right)\theta^3$ + higher order terms $=0$

This is uniquely satisfied to ${\cal O}(\theta^3)$ by $p_1=\displaystyle \frac12$      $q_1=-\displaystyle \frac12\quad c_3=-\displaystyle \frac{1}{12}$

Hence $\hat r_{1/1}\equiv\displaystyle \frac{1+\frac12 \theta}{1-\frac{1}{2}\theta}$ is a $(1,1)$ Padé Approximation of $e^{\theta}$ of order $2$. It has leading-order error = $-\displaystyle \frac{1}{12} \theta^3$.

In general

$$e^{\theta}=\frac{1+p_1\theta +p_2\theta^2+\cdots p_T\theta^T}{1+q... ...{S+T+1}}_{\mbox {constant.}}\theta^{S+T+1}+{\cal O}\left(\theta^{S+t+2}\right)$$

Hence $\hat r_{1/S}=\displaystyle \frac{P_T(\theta)}{Q_S(\theta)}$ is the $(T,S)$ Padé Appointment of order $T+S$ to $e^{\theta}$

Exercise: Show that

$(T,S)$	$\hat r_{T/S}$	Principal
		Error Term
$(1,0)$	$1+\theta$	$$\frac12\theta^2$$
$(2,0)$	$1+\theta+\frac12\theta^2$	$$\frac16\theta^3$$
$(2,1)$	$\frac{1+\displaystyle \frac23\theta+\frac16\theta^2}{1-\displaystyle \frac13\theta}$	$-\frac{1}{72}\theta^4$

$\Box$

Example:

Approximate ${\bf V}(t+k)=e^{kA}{\bf V}(t)$ using $\hat v_{1/1}$

$$\Rightarrow {\bf V}(t+k)=(I-\frac12kA)^{-1}(I+\frac12kA){\bf V}(t)$$

or $(I-\frac{1}{2}kA){\bf V}(t+k)=(I+\frac12 kA){\bf V}(t)$

or $-\displaystyle \mu u^{n+1}_{m-1}+2(1+\mu)u^{n+1}_m-\mu u^{n+1}_{m+1}=\mu u^n_{m-1}+ 2(1-\mu)u^n_m+\mu u^{n}_{m+1} m=1\cdots M-1$

$$\mbox {Crank-Nicholson!}$$

The $\hat r_{1/1}$ is also called a ``Unitary'' approximation, which is important property of the Schroedinger equation which is important to preserve.

$A$ and $L$ Stability

Continuing our discussion of $$\frac{\partial U}{\partial t}=\frac{\partial ^2 U}{\partial ^2 x}\quad t>0<X$$

with $u(0,x)=g(x)$ and assume for simplicity that boundary values are zero.

Take ${\bf V}(0)=[g_1, g_2\cdots g_{M-1}]^T\;$$\mbox { assume that } Mh=X$

$${\bf V}(t) \mbox { is an approximating vector to } {\bf U}(t)$$

$${\bf V}(t_n+k)=\hat r_{T/S} (k{\bf A}){\bf V}(t_n)$$

but

$${\bf V}(t_n)=\hat r_{T/S} (k{\bf A}){\bf A}(t_{n-1}) ,$$$$\mbox {and so on, }$$

which leads recursively to

$${\bf V}(t_n)=[\hat r_{T/S}(kA)]^n{\bf V}(0)$$ (150)

The eigenvalues of $A \Big($$\mbox {for center difference approx to } \displaystyle \frac{\partial ^2}{\partial x^2}\Big)$ are

$$-\frac{4}{h^2}\sin^2\left(\frac{\ell\pi}{2M}\right)\quad s=1,2, \cdots M-1$$

and are all different. Hence the eigenvectors of $A$ are independent and a basis for the $(M-1)$-dimensional space of the vector $g$ of initial values $\therefore$

$${\bf g}=\sum^{M-1}_{\ell=1} c_{\ell}\phi_{\ell}$$

$\therefore$ (151) can be expressed as

$${\bf V}(t_n)=[\hat r_{T/S}(kA)]^n\sum^{M-1}_{l=1}c_l\phi_l=\sum^{M-1}_{l=1} c_l[\hat r_{T/S}(kA)]^n\phi_l$$ (151)

Since $A\phi_l=\lambda_l\phi_l$ and we know that $f(A)\phi_l=f(\lambda_l)\phi_l$ it follows that (152) can be expressed as

$${\bf V}(t_n)=\sum^{M-1}_{l=1}c_l[\hat r_{T/S}(k\lambda_l)]^n\phi_l.$$ (152)

(153) shows that ${\bf V} (t_n)$ will tend to the null vector as $n\rightarrow\infty$ if and only if

$$\vert\hat r_{T/S}(k\lambda_l)\vert<1\quad l=1, 2, \cdots M-1$$

If this condition is subject to $\mu=k/h^2$ value, the equations are ``CONDITIONALLY STABLE.''

   When $$\vert\hat r_{T/S}(\lambda_lk)\vert<1\quad \forall\quad \mu>0\Rightarrow \lq\lq A-$$$$\mbox {Stable'' or unconditionally stable}$$

Although $A$ stability implies that $-1<\hat r_{T/S} (k\gamma_s)<1$ for real $\hat r_{T/S}$, it is possible that some values of $\hat r_{T/S (k\lambda_l)}$ be close to $-1$ and hence for these $\hat r_{T/S}(k\lambda_l)$ will alternate in sign as $n$ increases and diminish in amplitude only very slowly. This phenomenon is particularly pronounced in the $x$-neighborhoods of points of the discontinuity either in the initial values or between boundary and intial values.

The real coefficients of $\hat r_{T/S}$ would clearly be free of unwanted oscillations if $0<\hat r_{T/S}<1$ and $\hat r_{T/S}(k\lambda_l)\rightarrow 0$ monotonically as $k\lambda_l$ increase in magnitude. The $(0,1)$ Padé Approximant

$\hat r_{0/1}=\displaystyle \frac{1}{1-k\lambda_l} \quad , \quad \lambda_l$ real negative would clearly have this property. This corresponds to implicit (backwards) Euler.

If $\vert\hat r_{T/S}(k\lambda_\ell)\vert<1$ for      $l=1,\cdots M-1$     we say the scheme is $L_0$ stable

For Crank-Nicholson $$\hat r_{1/1}(-z)=\frac{1-\frac12 z}{1+\frac12 z}=\frac{2/z-1}{2/z+1}$$

$\vert\hat r_{1/1}(-z)\vert<1\quad \forall z>0$ but $\hat r_{1/1}(-z)\rightarrow -1$ as $z\rightarrow \infty$

$\therefore$ CN is $A-$ stable.

In order to avoid unwanted oscillations one can show that it is sufficient for $c_{m-1}\phi_{m-1}$ to decay to zero faster than the lowest component $c_1\phi_1$. This implies that

see Lawson, Morris (1978) J Num Anal SIAM 15, pp 1212-25.