Hello mathematicians,
I focus on the numerical solution of partial differential equations (PDEs). Recently, I was studying the numerical stability of method of lines (MOL).
Method of lines transforms a given PDE into the system of ordinary differential equations (ODEs). The ODEs then can be solved using standard numerical methods (e.g. Runge-Kutta, etc.).
I found this publication: Numerical stability of MOL for PDEs by R. Vichnevetsky: https://scholarship.libraries.rutgers.edu/esploro/outputs/technicalDocumentation/Numerical-stability-of-methods-of-lines/991031549877904646
There are stability conditions for:
- a) parabolic equation 1D (Tab. 1, p. 17)
- b) hyperbolic equation 1D (Tab. 2, p. 21)
My problem is: I would like to derive the stability condition for the 2nd order hyperbolic equation, and I would like to ensure that it is correctly derived.
Firstly, I review the stability condition for 1st order hyperbolic PDE (from the article above, p. 18): $$ \frac{\partial u}{\partial t} + f\frac{\partial u}{\partial x}=0 $$
In equation (26) from the article, we have ($c$ are the characteristic values of $f$): $$ \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0 $$
We apply central differences: $$ \left[\frac{\partial u}{\partial x}\right]_n \approx \frac{u_{n+1} - u_{n-1}}{2\Delta x} $$
Thus: $$ \frac{\partial u_n}{\partial t} \approx -\frac{c}{2\Delta x} (u_{n+1} - u_{n-1}) $$
Spectral function of the right-hand side: $$ A(\omega) = -\frac{c}{2\Delta x} (e^{j\omega \Delta x} - e^{-j\omega \Delta x}) $$
Because: $e^{j\omega \Delta x} - e^{-j\omega \Delta x} = 2j \sin \phi$, we get: $$ A(\omega) = -\frac{jc}{\Delta x} \sin (\omega \Delta x) $$
Numerical stability condition is: $$ |\Delta t A(\omega)|_{max} = \left | \frac{c \Delta t}{\Delta x} \sin (\omega \Delta x) \right |_{max} = \frac{c \Delta t}{\Delta x} \leq |S_I| $$ and is purely-imaginary for all $\omega$.
Note, that $S_I$ is the intersection of the stability boundary of the integration method with the imaginary axis.
Now, I am trying to derive a similar stability condition for the wave equation -- 2nd order hyperbolic PDE:
$$ \frac{\partial^2 u}{\partial t^2} + \alpha^2\frac{\partial^2 u}{\partial x^2}=0 $$
We apply central differences: $$ \left[\frac{\partial^2 u}{\partial x^2}\right]_n \approx \alpha^2 \left(\frac{u_{n+1} - 2u_n + u_{n-1}}{\Delta x^2}\right) $$
Spectral function of the right-hand side: $$ A(\omega) = \frac{\alpha^2}{\Delta x^2} (-2+ e^{j\omega \Delta x} + e^{-j\omega \Delta x}) $$
Because $e^{j\omega \Delta x} + e^{-j\omega \Delta x} = 2 \cos \omega$, we get:
$$ A(\omega) = \frac{\alpha^2}{\Delta x^2} (-2+2 \cos (\omega \Delta x)) \\ A(\omega) = -\frac{2\alpha^2}{\Delta x^2} (1-\cos (\omega \Delta x)) $$
Because the maximum amplitude of $\cos$ is $\pm1$, then the bracket is $(1-\cos (\omega \Delta x)) = 2$.
Numerical stability condition is: $$ |\Delta t^2 A(\omega)|_{max} = \left | \frac{4 \alpha^2 \Delta t^2}{\Delta x^2}\right | \leq |S_R| $$ and is real for all values of $\omega$.
Note, that $S_R$ is the intersection of the stability boundary of the integration method with the real axis.
I am modelling the heat equation, so I am using the $\alpha$ parameter, which represents thermal diffusivity (e.g. copper = $11.3e^{-5} m^/s$.
I would like to kindly ask you: Is this stability condition correct, please?
Thank you very much in advance. Any help is highly appreciated :-).
Note that you also have a second order derivative in $t$. So to apply a standard numerical method you have to construct a first-order system. To get a small matrix norm or to see eigenvalues easily, distribute the factor equally. $$ \pmatrix{\hat u_t(ω,t)\\\hat v_t(ω,t)} =\frac{2α\sin(ωΔx/2)Δt}{Δx} \pmatrix{0&1\\1&0}\pmatrix{\hat u(ω,t)\\\hat v(ω,t)} $$ As this gives pairs of eigenvalues of opposite sign, the exact solution is not stable, and thus this can also not be preserved by the numerical method.
Note that what you started with is the Poisson equation, the wave equation has the opposite sign and then above calculations lead again to the imaginary trace of the stability region.