The wave equation 1D acoustic and variable velocity $ V_j$ is defined as:
\begin{eqnarray}
\frac{\partial^2 U}{\partial t ^2} = V^2 \frac{\partial^2 U}{\partial x ^2}
\label{imp1}
\end{eqnarray}
And the implicit difference schema can be written as bellow for the unknown time $n+1$ (using the times $n$ and $n-1$ known):
\begin{eqnarray}
U_j^{n+1} = \left(\frac{\Delta t V_j^n}{\Delta s}\right)^2 \left( U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n - U_j^{n-1}
\label{imp2}
\end{eqnarray}
where $\Delta x = \Delta s $. \\
or as:
\begin{equation}
\left(U_j^{n-1} - 2 U_j^n \right) r_j^{2} = U_{j+1}^{n+1} - (2 + r_j^2) U_j^{n+1} + U_{j-1}^{n+1}
\end{equation}
With $r_j = \frac{\Delta s}{\Delta t V_j^n} $ \\
Applying the boundary condition $U_0^t=U_{N-1}^t=0$, the linear system of equations, in matrix form, for every unknown $ U_j^{n+1} $ for $j \in [0, N]$ being $NDs$ the 1D domain for the solution in $x$ is given by:
\begin{eqnarray}
\begin{pmatrix}
-(2+r_1^2) & 1 & & & \\
1 & -(2+r_2^2) & 1 & & \\
& 1 & -(2+r_3^2) & 1 & \\
& & 1 & -(2+r_4^2) & 1 \\
& & & 1 & -(2+r_5^2) \\
\end{pmatrix}
\begin{pmatrix}
\ U_{1}^{n+1} \ \\
\ U_{2}^{n+1} \ \\
\ U_{3}^{n+1} \ \\
\ U_{4}^{n+1} \ \\
\ U_{5}^{n+1} \ \\
\end{pmatrix}
=
\begin{pmatrix}
\left(U_j^{n-1} - 2 U_j^n \right) r_1^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_2^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_3^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_4^{2} \\
\left(U_j^{n-1} - 2 U_j^n \right) r_5^{2} \\
\end{pmatrix}
\end{eqnarray}
Simplified to the case where $N=7$ and $ U_{6} = U_{0} = 0 $
Von Neumann Stability Analysis
Since $ u(x,t) = e^{i(wt+kx)} $ is a solution for (\ref{imp1}) the discrete schema can be written as:
\begin{eqnarray}
u(x_j,t_n) &=& e^{i(wt_n+kx_j)} \nonumber \\
&=& e^{iwn\Delta t} e^{ikj\Delta s} \nonumber \\
&=& \epsilon^n e^{ikj\Delta s}
\label{imp3}
\end{eqnarray}
Where $\epsilon = e^{iw \Delta t} $ and $ i = \sqrt{-1} $ .
To maintain stability we should make sure that $ \epsilon \leq 1 $ not growing exponentially with increasing time steps. So applying (\ref{imp3}) in (\ref{imp2}) we can analyse the growth factor $\epsilon $ .
\begin{equation}
U_j^{n+1} = \left( \frac{\Delta t V_j^n}{\Delta s} \right) ^2 \left( U_{j+1}^{n+1} - 2 U_j^{n+1} + U_{j-1}^{n+1} \right) + 2 U_j^n - U_j^{n-1} \nonumber
\end{equation}
\begin{equation}
\begin{split}
\epsilon^{n+1} e^{ikj\Delta s} &= r^2 \left( \epsilon^{n+1} e^{ik(j+1)\Delta s} - 2 \epsilon^{n+1} e^{ikj\Delta s} + \epsilon^{n+1} e^{ik(j-1)\Delta s} \right) \\
&\quad + 2 \epsilon^n e^{ikj\Delta s} - \epsilon^{n-1} e^{ikj\Delta s} \nonumber
\end{split}
\end{equation}
\begin{eqnarray}
1 = r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 2 \epsilon^{-1} - \epsilon^{-2} \nonumber \\
\epsilon^{-2} - 2 \epsilon^{-1} - r^2 (e^{ik\Delta s} -2 + e^{-ik\Delta s}) + 1 = 0 \nonumber \\
\epsilon^{-2} - 2 \epsilon^{-1} + 4 r^2 \sin^2 \left( k \Delta s / 2 \right) + 1 = 0
\label{imp4}
\end{eqnarray}
Where $ r = \left( \frac{\Delta t V_j^n}{\Delta s} \right) $
At equation (\ref{imp4}) we can substitute $ \phi = \epsilon^{-1} $ turning it to a second degree bellow:
\begin{eqnarray}
\phi^2 - 2 \phi +4 r^2 \sin^2 \left(k\Delta s/2 \right) +1 = 0 \nonumber \\
\phi^2 - 2 \phi + c = 0 \nonumber
\end{eqnarray}
With:
$$ c = 1 +4 r^2 \sin^2 \left(k\Delta s/2 \right) $$
Has roots $ \phi^{'} $ and $ \phi^{''} $ as bellow. Replacing $c$ also.
\begin{eqnarray}
&=& \frac{2 \pm \sqrt{4 - 4c}}{2} = 1 \pm \sqrt{1-c} \nonumber \\
&=& 1 \pm \sqrt{-4 r^2 \sin^2 \left(k\Delta s/2 \right) } \nonumber \\
&=& 1 \pm i \ 2 r \sin \left(k\Delta s/2 \right) \nonumber \\
\phi^{'} &=& 1 + i \ 2 r \sin \left(k\Delta s/2 \right) \nonumber \\
\phi^{''} &=& 1 - i \ 2 r \sin \left(k\Delta s/2 \right) \nonumber
\end{eqnarray}
Analyzing the modulus, since $ \phi^{'} $ and $ \phi^{''} $ are conjugate pairs of the same complex number the modulus is the same for both.
We get:
$$ \| \phi \| = \| \phi^{'} \| = \| \phi^{''} \| = \sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}$$
Going back to $\epsilon$ and getting its modulus we get :
\begin{eqnarray}
\| \epsilon^{'} \| &=& \| \epsilon^{''} \| = \frac{1}{\| \phi \|} \nonumber \\
&=& \frac{1}{\sqrt{1+ 4 r^2 \sin ^2 \left(k\Delta s/2 \right)}}
\label{imp5}
\end{eqnarray}
Looking at (\ref{imp5}) we can see that always $ r^2 > 0 $ and $ sin^2(x) \in \ [0,1] $ so:
$$ \| \epsilon \| \leq 1 \ \ \forall \ r, \ k,\ \Delta s$$
It shows that FTCS in (\ref{imp2}) is unconditionally stable
Comments
Post a Comment