MATH 761 Lesson 02: ODE methods applied to semi-discretized system

Overview

• Semi-discretized PDE system behavior

  • Parabolic case revisited: Crank-Nicolson

  • Hyperbolic case: wave equation

  • Elliptic case: Poisson equation

Parabolic case: Crank-Nicolson

• For $q_t=\alpha q_{xx}$, forward Euler in time with centered in space finite differencing (second-order) led to time step condition $\Delta t⩽h^2/2$. This is computationally prohibitive, leading to small time steps as $h\to 0$

• Try an implicit-in-time scheme. Since spatial discretization is second-order in space try trapezoid which is second-order in time

  ${\displaystyle {𝒒}^{n+1}-{𝒒}^n=\frac{\alpha \Delta t}{2h^2}𝑩({𝒒}^{n+1}+{𝒒}^n),𝑩=\mathrm{diag}\left(\left[\begin{array}{lll}1& -2& 1\end{array}\right]\right)}$ ${\displaystyle {𝒒}^n={\left[\begin{array}{lll}{q}_1\left(t^n\right)& \dots & q_m\left(t^n\right)\end{array}\right]}^T,q_j\left(t\right)\cong q(t,x_j),x_j=jh,t^n=n\Delta t}$

• This is Crank-Nicolson, requiring solution of a linear system at each time step

  ${\displaystyle (𝑰-\nu 𝑩){𝒒}^{n+1}=(𝑰+\nu 𝑩){𝒒}^n}$

  Since linear system is tridiagonal, computational effort is small $𝒪\left(m\right)$ per time step

• Stability criterion: $𝑨𝒙=\lambda 𝒙$ with $\lambda \in ℂ$ within region of stability of numerical scheme

• Since ${𝒒}^{\text{'}}=𝑨𝒒$ is discretization of $q_t=q_{xx}$, eigenvalues of $𝑨$ related to those of ${\partial }_x^2$

• Example: Forward Euler in time with time step $\Delta t$, centered in space with step $h$

  ${\displaystyle q_j^{n+1}=q_j^n+\nu (q_{j+1}^n-2q_j^n+q_{j-1}^n),\nu =\frac{\alpha \Delta t}{h^2}\Rightarrow 𝑨=\mathrm{diag}\left(\left[\begin{array}{lll}\nu & 1-2\nu & \nu \end{array}\right]\right)}$

• Eigenvalues of $𝑩=\mathrm{diag}\left(\left[\begin{array}{lll}1& -2& 1\end{array}\right]\right)$, guess eigenvector is discretization of $e^{ik{x}_j}=e^{ijkh}$

  ${\displaystyle 𝑩{𝒚}_k={\mu }_k{𝒚}_k,y_j=e^{ijkh}\Rightarrow e^{i(j-1)kh}-2e^{ijkh}+e^{i(j+1)kh}={\mu }_k{e}^{ijkh}\Rightarrow }$ ${\displaystyle {\mu }_k=e^{-ikh}+e^{ikh}-2=2(\cos kh-1)=-4{\sin }^{2}kh}$

• Eigenvalues of $𝑨=𝑰+\nu 𝑩$ are

  ${\displaystyle {\lambda }_k=1+\nu {\mu }_k=1-\frac{4\Delta t}{h^2}{\sin }^{2}kh,\left|{\lambda }_k\right|>1\Rightarrow \Delta t⩽\frac{h^2}{2}}$