 |
 |
The following are animated GIFs showing the solution to problem 3, parts (b) and (c). These will show only in the HTML version.
Assuming that \(f\left ( v\right ) =\left ( a-v\right ) \left ( v-1\right ) v\), the equations solved are the following\begin {align*} \frac {\partial v}{\partial t} & =D\Delta v+f\left ( v\right ) -w+I\\ \frac {\partial w}{\partial t} & =\epsilon \left ( v-\gamma w\right ) \end {align*}
Click on image to see the animation run, it will open in new window.
| |
| |
The diffusion PDE is given by\[ u_{t}-D\Delta u=0 \] Where \(D\) is the diffusion constant. The ADI scheme in 3D 6 is given by\begin {align} \left ( I-\frac {D\Delta t}{3}L_{x}\right ) u^{\ast } & =\left ( I+\frac {D\Delta t}{3}L_{y}+\frac {D\Delta t}{3}L_{z}\right ) u^{n}\tag {1}\\ \left ( I-\frac {D\Delta t}{3}L_{y}\right ) u^{\ast \ast } & =\left ( I+\frac {D\Delta t}{3}L_{x}+\frac {D\Delta t}{3}L_{z}\right ) u^{\ast }\tag {2}\\ \left ( I-\frac {D\Delta t}{3}L_{z}\right ) u^{n+1} & =\left ( I+\frac {D\Delta t}{3}L_{x}+\frac {D\Delta t}{3}L_{y}\right ) u^{\ast \ast } \tag {3} \end {align}
Where \(L_{x},L_{y},L_{z}\) are each the 1D Laplacian given by \(\frac {1}{h^{2}}\begin {bmatrix} -2 & 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ \cdots & \cdots & \cdots & \ddots & \cdots & \cdots \\ 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 & -2 \end {bmatrix} \)
Assuming that the spatial frequencies in each of the three Cartesian directions (\(x,y,z\)) are given by \(\xi _{x},\xi _{y},\xi _{z}\) where \(\frac {\pi }{h}\leq \xi _{i}\leq \frac {\pi }{h}\) and by setting \(r=\frac {D\Delta t}{3h^{2}}\), \(u^{\ast }=g^{\ast }e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and \(u=e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and substituting these into Eq. (1) and dividing throughout by \(e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) gives the following\begin {align} g^{\ast }\left ( 1-r\left ( e^{-i\xi _{1}h}-2+e^{i\xi _{1}h}\right ) \right ) & =1+r\left ( e^{-i\xi _{2}h}-2+e^{i\xi _{2}h}\right ) +r\left ( e^{-i\xi _{3}h}-2+e^{i\xi _{3}h}\right ) \nonumber \\ g^{\ast } & =\frac {1-4r+r\left ( e^{-i\xi _{2}h}+e^{i\xi _{2}h}\right ) +r\left ( e^{-i\xi _{3}h}+e^{i\xi _{3}h}\right ) }{\left ( 1+2r-r\left ( e^{i\xi _{1}h}+e^{-i\xi _{1}h}\right ) \right ) }\nonumber \\ & =\frac {1-4r+2r\cos \left ( \xi _{2}h\right ) +2r\cos \left ( \xi _{3}h\right ) }{1+2r-2r\cos \left ( \xi _{1}h\right ) }\nonumber \\ & =\frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) } \tag {4} \end {align}
The last step above was obtained by the use of the relation \(\cos A=1-2\sin ^{2}\left ( \frac {A}{2}\right ) \).
Applying the same method used above to Eq. (2), but now letting \(u^{\ast }=g^{\ast }e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and \(u^{\ast \ast }=g^{\ast \ast }e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and dividing throughout by \(e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) gives\begin {align} g^{\ast \ast }\left ( 1-r\left ( e^{-i\xi _{2}h}-2+e^{i\xi _{2}h}\right ) \right ) & =1+r\left ( e^{-i\xi _{1}h}-2+e^{i\xi _{1}h}\right ) +r\left ( e^{-i\xi _{3}h}-2+e^{i\xi _{3}h}\right ) \nonumber \\ g^{\ast \ast } & =\frac {1-4r+r\left ( e^{-i\xi _{1}h}+e^{i\xi _{1}h}\right ) +r\left ( e^{-i\xi _{3}h}+e^{i\xi _{3}h}\right ) }{1+2r-r\left ( e^{i\xi _{2}h}+e^{-i\xi _{2}h}\right ) }g^{\ast }\nonumber \\ & =\frac {1-4r+2r\cos \left ( \xi _{1}h\right ) +2r\cos \left ( \xi _{3}h\right ) }{1+2r-2r\cos \left ( \xi _{2}h\right ) }g^{\ast }\nonumber \\ & =\frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) }g^{\ast } \tag {5} \end {align}
Again, applying the same method to Eq. (3), but now letting \(u^{\ast \ast }=g^{\ast \ast }e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and \(u^{n+1}=ge^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) and dividing by \(e^{i\left ( \xi _{1}x+\xi _{2}y+\xi _{3}z\right ) }\) gives\begin {align} g\left ( 1-r\left ( e^{-i\xi _{3}h}-2+e^{i\xi _{3}h}\right ) \right ) & =1+r\left ( e^{-i\xi _{1}h}-2+e^{i\xi _{1}h}\right ) +r\left ( e^{-i\xi _{2}h}-2+e^{i\xi _{2}h}\right ) \nonumber \\ g & =\frac {1-4r+r\left ( e^{-i\xi _{1}h}+e^{i\xi _{1}h}\right ) +r\left ( e^{-i\xi _{2}h}+e^{i\xi _{2}h}\right ) }{1-r\left ( e^{-i\xi _{3}h}-2+e^{i\xi _{3}h}\right ) }g^{\ast \ast }\nonumber \\ & =\frac {1-4r+2r\cos \left ( \xi _{1}h\right ) +2r\cos \left ( \xi _{2}h\right ) }{1+2r-2r\cos \left ( \xi _{3}h\right ) }g^{\ast \ast }\nonumber \\ & =\frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) }g^{\ast \ast } \tag {6} \end {align}
Substituting (4) into (5) and substituting the resulting expression into (6) gives the overall magnification factor for the ADI scheme: \begin {equation} g=\left [ \frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) }\right ] \left [ \frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) }\right ] \left [ \frac {1-4r\left ( \sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) +\sin ^{2}\left ( \frac {\xi _{3}h}{2}\right ) \right ) }{1+4r\sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) }\right ] \tag {7} \end {equation} Letting \(A\equiv \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) ,B\equiv \sin ^{2}\left ( \frac {\xi _{2}h}{2}\right ) ,C\equiv \sin ^{2}\left ( \frac {\xi _{1}h}{2}\right ) =C\) in Eq. (7) results in\begin {equation} g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) =\left ( \frac {1-4r\left ( A+B\right ) }{1+4rC}\right ) \left ( \frac {1-4r\left ( A+C\right ) }{1+4rB}\right ) \left ( \frac {1-4r\left ( B+C\right ) }{1+4rA}\right ) \tag {8} \end {equation} The scheme is conditionally stable if \(\left \vert g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) \right \vert \leq 1\) for some value of \(r\) and \(\left \vert g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) \right \vert >1\) for some other value of \(r\) (this is the same as using different values of \(\Delta t\,\) in place of \(r\), since \(r=\frac {D\Delta t}{3h^{2}}\) and \(h\) and \(D\) would be kept constant).
Now the scheme can be shown to be conditionally stable by letting \(A=B=C=1\) in Eq. (8) and then by finding one value of \(r\) which makes the magnification factor to become less than one and then by looking for another value of \(r\) which makes the magnification factor to becomes larger than one.
Therefore, when \(A=B=C=1\), Eq. (8) becomes\begin {equation} \left \vert g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) \right \vert =\left ( \frac {1-8r}{1+4r}\right ) \left ( \frac {1-8r}{1+4r}\right ) \left ( \frac {1-8r}{1+4r}\right ) \tag {8A} \end {equation} Now, putting \(r=2\) in the above gives \(\left \vert g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) \right \vert =2.744>1\) implying that the scheme is unstable.
Putting \(r=0.5\) in Eq. (8A) gives \(\left \vert g\left ( \xi _{1},\xi _{2},\xi _{3}\right ) \right \vert =0.125<1\) implying that the scheme is stable.
Hence the scheme is conditionally stable, because by fixing \(h\) and \(D\), it was possible to find a time step \(\Delta t\) which made some mode become unstable. If one mode is unstable, the overall scheme is also unstable. This result shows that the above given ADI scheme for 3D is conditionally stable.
Expectation: Temporal accuracy is expected to be \(O\left ( \Delta t\right ) \) since at each \(1/3\) time step there is one implicit step compared to two explicit steps. Starting from the main equations shown in part (a)\begin {align} \overset {implicit\ \text {(backward\ Euler)}}{\overbrace {\left ( I-rL_{x}\right ) u^{\ast }}} & =\overset {\text {explicit (2 forward Euler)}}{\overbrace {\left ( I+rL_{y}+rL_{z}\right ) }}u^{n}\tag {1}\\ \left ( I-rL_{y}\right ) u^{\ast \ast } & =\left ( I+rL_{x}+rL_{z}\right ) u^{\ast }\tag {2}\\ \left ( I-rL_{z}\right ) u^{n+1} & =\left ( I+rL_{x}+rL_{y}\right ) u^{\ast \ast } \tag {3} \end {align}
There will be an \(O\left ( \Delta t\right ) \) error resulting from the application of Euler approximation to each of the terms in each equation above. One of the implicit errors will cancel exactly with one of the errors from the explicit part of the equation (due to sign difference), leaving an extra \(O\left ( \Delta t\right ) \) error after each third step. Hence at the completion of one full time step, the temporal error will be \(3O\left ( \Delta t\right ) \) or \(O\left ( \Delta t\right ) \).
Explanation: The derivation below follows the method explained in class for the 2D ADI case, but being applied to the 3D case. Starting by pre-multiplying Eq. (1) by the operator \(\left ( I+rL_{x}+rL_{z}\right ) \) gives\[ \left ( I+rL_{x}+rL_{z}\right ) \left ( I-rL_{x}\right ) u^{\ast }=\left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] But since \(\left ( I+rL_{x}+rL_{z}\right ) \) commutes7 with \(\left ( I-rL_{x}\right ) \), then the two terms in the LHS of the above equation can be interchanged giving\[ \left ( I-rL_{x}\right ) \overset {\text {now replace this from (2)}}{\overbrace {\left ( I+rL_{x}+rL_{z}\right ) u^{\ast }}}=\left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] Replacing the term marked above by its LHS value from Eq. (2) yields\[ \left ( I-rL_{x}\right ) \left ( I-rL_{y}\right ) u^{\ast \ast }=\left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] Pre-multiplying the above by the operator \(\left ( I+rL_{x}+rL_{y}\right ) \) gives\[ \left ( I+rL_{x}+rL_{y}\right ) \left ( I-rL_{x}\right ) \left ( I-rL_{y}\right ) u^{\ast \ast }=\left ( I+rL_{x}+rL_{y}\right ) \left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] But since \(\left ( I+rL_{x}+rL_{y}\right ) \) commutes with \(\left ( I-rL_{x}\right ) \left ( I-rL_{y}\right ) \) the above can be written as\[ \left ( I-rL_{x}\right ) \left ( I-rL_{y}\right ) \overset {\text {replace this from (3)}}{\overbrace {\left ( I+rL_{x}+rL_{y}\right ) u^{\ast \ast }}}=\left ( I+rL_{x}+rL_{y}\right ) \left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] Replacing the term marked above by its LHS value from Eq. (3) gives\[ \left ( I-rL_{x}\right ) \left ( I-rL_{y}\right ) \left ( I-rL_{z}\right ) u^{n+1}=\left ( I+rL_{x}+rL_{y}\right ) \left ( I+rL_{x}+rL_{z}\right ) \left ( I+rL_{y}+rL_{z}\right ) u^{n}\] Expanding all terms by multiplying all operators and simplifying the result and using \(L=L_{x}+L_{y}+L_{z}\) gives the following
 |
Since \(r=\frac {D\Delta t}{3}\) the above becomes\[ u^{n+1}-u^{n}=\frac {D\Delta t}{3}Lu^{n+1}+2\frac {D\Delta t}{3}Lu^{n}+O\left ( \Delta t^{2}\right ) +O\left ( \Delta t^{3}\right ) +\cdots \] Dividing the above equation by \(\Delta t\) gives\[ \frac {u^{n+1}-u^{n}}{\Delta t}=\frac {D}{3}Lu^{n+1}+\frac {2D}{3}Lu^{n}+O\left ( \Delta t\right ) +O\left ( \Delta t^{2}\right ) +\cdots \] Now adding \(\frac {D}{6}Lu^{n+1}\) and subtracting \(\frac {D}{6}Lu^{n+1}\) and subtracting \(\frac {D}{6}Lu^{n}\) and adding \(\frac {D}{6}Lu^{n}\) from the RHS of the above equation gives\begin {align*} \frac {u^{n+1}-u^{n}}{\Delta t} & =\frac {D}{3}Lu^{n+1}+\frac {D}{6}Lu^{n+1}+\frac {2D}{3}Lu^{n}-\frac {D}{6}Lu^{n}+\frac {D}{6}Lu^{n}-\frac {D}{6}Lu^{n+1}+O\left ( \Delta t\right ) +O\left ( \Delta t^{2}\right ) +\cdots \\ & \overset {C-N}{\overbrace {\frac {u^{n+1}-u^{n}}{\Delta t}=\frac {1}{2}\left ( Lu^{n+1}+Lu^{n}\right ) }}+\frac {1}{6}\left ( Lu^{n}-Lu^{n+1}\right ) +O\left ( \Delta t\right ) +O\left ( \Delta t^{2}\right ) +\cdots \end {align*}
The C-N scheme is known to be \(O\left ( \Delta t^{2}+h^{2}\right ) \). Multiplying the term \(\frac {1}{6}\left ( Lu^{n}-Lu^{n+1}\right ) \) by \(\frac {\Delta t}{\Delta t}\) in the above yields\[ u_{t}=\frac {1}{2}\Delta u+\overset {\text {from C-N part}}{\overbrace {O\left ( \Delta t^{2}\right ) +O\left ( h^{2}\right ) }}+\frac {\Delta t}{6}L\left ( \frac {u^{n}-u^{n+1}}{\Delta t}\right ) +O\left ( \Delta t\right ) +O\left ( \Delta t^{2}\right ) +\cdots \] Taking the limits \(\Delta t\rightarrow 0\) results in\begin {align*} u_{t} & =\frac {1}{2}\Delta u+\overset {\text {from C-N part}}{\overbrace {O\left ( \Delta t^{2}\right ) +O\left ( h^{2}\right ) }}+\overset {still\ O\left ( \Delta t\right ) }{\overbrace {\frac {\Delta t}{6}\frac {\partial ^{7}}{\partial ^{2}x\partial ^{2}y\partial ^{2}z\partial t}+O\left ( \Delta t\right ) }}+O\left ( \Delta t^{2}\right ) +\cdots \\ & =\frac {1}{2}\Delta u+\overset {\text {from C-N part}}{\overbrace {O\left ( \Delta t^{2}\right ) +O\left ( h^{2}\right ) }}+O\left ( \Delta t\right ) +O\left ( \Delta t^{2}\right ) +\cdots \\ & =\frac {1}{2}\Delta u+\overbrace {O\left ( h^{2}\right ) +O\left ( \Delta t\right ) }+O\left ( \Delta t^{2}\right ) +\cdots \end {align*}
Since in the above, the dominant temporal error term is \(O\left ( \Delta t\right ) \) the scheme is a first order in time accurate. It is also a second order in space accurate.
The following diagram shows the discretization using cell-centered scheme for the case of \(N=4\). The center of the cells moves closer to the physical boundaries of the unit square as \(N\) becomes larger.
The physical domain is always the unit square \(x=0\cdots 1\), \(y=0\cdots 1\), but the discrete solution domain is the one at corners of the red grid above. A small example is used below to help determine the layout of the operator used in the direct solver. The 2D ADI scheme for the diffusion problem is\begin {align} \left ( I-\frac {D\Delta t}{2}L_{x}\right ) \mathbf {u}^{\ast } & =\left ( I+\frac {D\Delta t}{2}L_{y}\right ) \mathbf {u}^{n}\tag {1A}\\ \left ( I-\frac {D\Delta t}{2}L_{y}\right ) \mathbf {u}^{n+1} & =\left ( I+\frac {D\Delta t}{2}L_{x}\right ) \mathbf {u}^{\ast }\nonumber \end {align}
Where \(L_{x}=L_{y}=\frac {1}{h^{2}}\begin {bmatrix} -1 & 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ \cdots & \cdots & \cdots & \ddots & \cdots & \cdots \\ 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 & -1 \end {bmatrix} \) for the case of homogenous Neumann boundary conditions. The solution given below uses an overall \(L\) operator which is used by the direct solver. Another approach would have been to use the above \(L_{x},L_{y}\) operator, and iterate over each each row and column applying the direct solver each time.
The following derives the overall \(L\) operator used. Eq. (1A) can be written as\begin {align*} u_{ij}^{\ast }-\frac {D\Delta t}{2}\frac {u_{i-1,j}^{\ast }-2u_{ij}^{\ast }+u_{i+1,j}^{\ast }}{h^{2}} & =u_{ij}^{n}+\frac {D\Delta t}{2}\frac {u_{i,j-1}^{n}-2u_{ij}^{n}+u_{i,j+1}^{n}}{h^{2}}\\ u_{ij}^{n+1}-\frac {D\Delta t}{2}\frac {u_{i,j-1}^{n+1}-2u_{ij}^{n+1}+u_{i,j+1}^{n+1}}{h^{2}} & =u_{ij}^{n\ast }+\frac {D\Delta t}{2}\frac {u_{i-1,j}^{\ast }-2u_{ij}^{\ast }+u_{i+1,j}^{\ast }}{h^{2}} \end {align*}
letting \(r=\frac {D\Delta t}{2h^{2}}\) and simplifying the above gives \begin {align} u_{ij}^{\ast }\left ( 1+2r\right ) -ru_{i-1,j}^{\ast }-ru_{i+1,j}^{\ast } & =u_{ij}^{n}\left ( 1-2r\right ) +ru_{i,j-1}^{n}+ru_{i,j+1}^{n}\tag {1}\\ u_{ij}^{n+1}\left ( 1+2r\right ) -ru_{i,j-1}^{n+1}-ru_{i,j+1}^{n+1} & =u_{ij}^{\ast }\left ( 1-2r\right ) +ru_{i-1,j}^{\ast }+ru_{i+1,j}^{\ast } \tag {2} \end {align}
The above finite difference equations are applied at all the grid points, except for those for the rows and columns at the boundaries. In order to determine the equations to use for the boundary grid points, the approximation \(L_{x}\approx \frac {-u_{1}+u_{2}}{h^{2}}\) is used. Similar one is used for \(L_{y}.\) The result of using the above approximation is the following finite difference equations used for the boundary grid points\begin {align*} u_{ij}^{\ast }-\frac {D\Delta t}{2}\frac {-u_{ij}^{\ast }+u_{i+1,j}^{\ast }}{h^{2}} & =u_{ij}^{n}+\frac {D\Delta t}{2}\frac {-u_{ij}^{n}+u_{i,j+1}^{n}}{h^{2}}\\ u_{ij}^{n+1}-\frac {D\Delta t}{2}\frac {-u_{ij}^{n+1}+u_{i,j+1}^{n+1}}{h^{2}} & =u_{ij}^{n\ast }+\frac {D\Delta t}{2}\frac {-u_{ij}^{\ast }+u_{i+1,j}^{\ast }}{h^{2}} \end {align*}
Simplifying the above gives\begin {align} u_{ij}^{\ast }\left ( 1+r\right ) -ru_{i+1,j}^{\ast } & =u_{ij}^{n}\left ( 1-r\right ) +ru_{i,j+1}^{n}\tag {1A}\\ u_{ij}^{n+1}\left ( 1+r\right ) -ru_{i,j+1}^{n+1} & =u_{ij}^{\ast }\left ( 1-r\right ) +ru_{i+1,j}^{\ast } \tag {2A} \end {align}
To help obtain \(L\), a small example is used to help see the structure of the matrix. This small example exhibits all the needed information to generate the pattern for \(L\) from. Using \(n_{x}=n_{y}=4\), the grid becomes
Eq. (1) and (1A) are now written for all the nodes resulting in the following 16 equations \begin {align*} u_{11}^{\ast }\left ( 1+r\right ) -ru_{21}^{\ast } & =u_{11}^{n}\left ( 1-r\right ) +ru_{12}^{n}\\ u_{21}^{\ast }\left ( 1+2r\right ) -ru_{1,1}^{\ast }-ru_{3,1}^{\ast } & =u_{21}^{n}\left ( 1-r\right ) +ru_{22}^{n}\\ u_{31}^{\ast }\left ( 1+2r\right ) -ru_{2,1}^{\ast }-ru_{4,1}^{\ast } & =u_{31}^{n}\left ( 1-r\right ) +ru_{32}^{n}\\ u_{41}^{\ast }\left ( 1+r\right ) -ru_{31}^{\ast } & =u_{41}^{n}\left ( 1-r\right ) +ru_{42}^{n}\\ u_{12}^{\ast }\left ( 1+r\right ) -ru_{22}^{\ast } & =u_{12}^{n}\left ( 1-2r\right ) +ru_{1,1}^{n}+ru_{1,3}^{n}\\ u_{22}^{\ast }\left ( 1+2r\right ) -ru_{1,2}^{\ast }-ru_{3,2}^{\ast } & =u_{22}^{n}\left ( 1-2r\right ) +ru_{2,1}^{n}+ru_{2,3}^{n}\\ u_{32}^{\ast }\left ( 1+2r\right ) -ru_{2,2}^{\ast }-ru_{4,2}^{\ast } & =u_{32}^{n}\left ( 1-2r\right ) +ru_{3,1}^{n}+ru_{3,3}^{n}\\ u_{42}^{\ast }\left ( 1+r\right ) -ru_{32}^{\ast } & =u_{42}^{n}\left ( 1-2r\right ) +ru_{4,1}^{n}+ru_{4,3}^{n}\\ u_{13}^{\ast }\left ( 1+r\right ) -ru_{23}^{\ast } & =u_{13}^{n}\left ( 1-2r\right ) +ru_{1,2}^{n}+ru_{1,4}^{n}\\ u_{23}^{\ast }\left ( 1+2r\right ) -ru_{1,3}^{\ast }-ru_{3,3}^{\ast } & =u_{23}^{n}\left ( 1-2r\right ) +ru_{2,2}^{n}+ru_{2,4}^{n}\\ u_{33}^{\ast }\left ( 1+2r\right ) -ru_{2,3}^{\ast }-ru_{4,3}^{\ast } & =u_{33}^{n}\left ( 1-2r\right ) +ru_{3,2}^{n}+ru_{3,4}^{n}\\ u_{43}^{\ast }\left ( 1+r\right ) -ru_{33}^{\ast } & =u_{43}^{n}\left ( 1-2r\right ) +ru_{4,2}^{n}+ru_{4,4}^{n}\\ u_{14}^{\ast }\left ( 1+r\right ) -ru_{24}^{\ast } & =u_{14}^{n}\left ( 1-r\right ) +ru_{13}^{n}\\ u_{24}^{\ast }\left ( 1+2r\right ) -ru_{1,4}^{\ast }-ru_{3,4}^{\ast } & =u_{24}^{n}\left ( 1-r\right ) +ru_{23}^{n}\\ u_{34}^{\ast }\left ( 1+2r\right ) -ru_{2,4}^{\ast }-ru_{4,4}^{\ast } & =u_{34}^{n}\left ( 1-r\right ) +ru_{33}^{n}\\ u_{44}^{\ast }\left ( 1+r\right ) -ru_{34}^{\ast } & =u_{44}^{n}\left ( 1-r\right ) +ru_{43}^{n} \end {align*}
In matrix form, the above gives \(Au^{\ast }=b\) which is then used to solve for \(u^{\ast }.\) The matrix \(A\) is now written out. To save space and to allow the matrix to fit on the page, the following terms are used \begin {align*} r & =\frac {D\Delta t}{2h^{2}}\\ \alpha & \equiv 1+r\\ \beta & \equiv 1+2r\\ \gamma & \equiv 1-r\\ \theta & \equiv 1-2r \end {align*}
 |
Starting the second ADI step to find \(u^{n+1}\), the process is similar to the one shown above, but the equations are written column-wise instead of row-wise as was the case earlier. For non-boundary grid points the following equation is used\[ u_{ij}^{n+1}\left ( 1+2r\right ) -ru_{i,j-1}^{n+1}-ru_{i,j+1}^{n+1}=u_{ij}^{\ast }\left ( 1-2r\right ) +ru_{i-1,j}^{\ast }+ru_{i+1,j}^{\ast }\] And for the boundary grid points the following equation is used\[ u_{ij}^{n+1}\left ( 1+r\right ) -ru_{i,j+1}^{n+1}=u_{ij}^{\ast }\left ( 1-r\right ) +ru_{i+1,j}^{\ast }\] Applying the above to each grid point results in the the following 16 equations \begin {align*} u_{11}^{n+1}\left ( 1+r\right ) -ru_{12}^{n+1} & =u_{11}^{\ast }\left ( 1-r\right ) +ru_{21}^{\ast }\\ u_{21}^{n+1}\left ( 1+r\right ) -ru_{22}^{n+1} & =u_{21}^{\ast }\left ( 1-2r\right ) +ru_{1,1}^{n}+ru_{3,1}^{n}\\ u_{31}^{n+1}\left ( 1+r\right ) -ru_{32}^{n+1} & =u_{31}^{\ast }\left ( 1-2r\right ) +ru_{2,1}^{n}+ru_{4,1}^{n}\\ u_{41}^{n+1}\left ( 1+r\right ) -ru_{42}^{n+1} & =u_{41}^{\ast }\left ( 1-r\right ) +ru_{31}^{\ast }\\ u_{12}^{n+1}\left ( 1+2r\right ) -ru_{1,1}^{n+1}-ru_{1,3}^{n+1} & =u_{12}^{\ast }\left ( 1-r\right ) +ru_{22}^{\ast }\\ u_{22}^{n+1}\left ( 1+2r\right ) -ru_{2,1}^{n+1}-ru_{2,3}^{n+1} & =u_{22}^{\ast }\left ( 1-2r\right ) +ru_{1,2}^{n}+ru_{3,2}^{n}\\ u_{32}^{\ast }\left ( 1+2r\right ) -ru_{3,1}^{\ast }-ru_{3,3}^{\ast } & =u_{32}^{n}\left ( 1-2r\right ) +ru_{2,2}^{n}+ru_{4,2}^{n}\\ u_{42}^{\ast }\left ( 1+2r\right ) -ru_{4,1}^{\ast }-ru_{4,3}^{\ast } & =u_{42}^{n}\left ( 1-r\right ) +ru_{32}^{n}\\ u_{13}^{\ast }\left ( 1+2r\right ) -ru_{1,2}^{\ast }-ru_{1,4}^{\ast } & =u_{13}^{n}\left ( 1-r\right ) +ru_{23}^{n}\\ u_{23}^{\ast }\left ( 1+2r\right ) -ru_{2,2}^{\ast }-ru_{2,4}^{\ast } & =u_{23}^{n}\left ( 1-2r\right ) +ru_{1,3}^{n}+ru_{3,3}^{n}\\ u_{33}^{\ast }\left ( 1+2r\right ) -ru_{3,2}^{\ast }-ru_{3,4}^{\ast } & =u_{33}^{n}\left ( 1-2r\right ) +ru_{2,3}^{n}+ru_{4,3}^{n}\\ u_{43}^{\ast }\left ( 1+2r\right ) -ru_{4,2}^{\ast }-ru_{4,4}^{\ast } & =u_{43}^{n}\left ( 1-r\right ) +ru_{33}^{n}\\ u_{14}^{\ast }\left ( 1+r\right ) -ru_{13}^{\ast } & =u_{14}^{n}\left ( 1-r\right ) +ru_{24}^{n}\\ u_{24}^{\ast }\left ( 1+r\right ) -ru_{23}^{\ast } & =u_{24}^{n}\left ( 1-2r\right ) +ru_{1,4}^{n}+ru_{3,4}^{n}\\ u_{34}^{\ast }\left ( 1+r\right ) -ru_{33}^{\ast } & =u_{34}^{n}\left ( 1-2r\right ) +ru_{2,4}^{n}+ru_{4,4}^{n}\\ u_{44}^{\ast }\left ( 1+r\right ) -ru_{43}^{\ast } & =u_{44}^{n}\left ( 1-r\right ) +ru_{34}^{n} \end {align*}
The above equations are now written as \(Au=b\) but the unknowns are listed column-wise in order to keep the tridiagonal form. The resulting matrix \(A\) is the following
 |
Since the edges of the domain are insulated, no concentration will diffuse to the outside. Therefore the result of diffusion will be that the concentration will diffuse internally and will spread out. Therefore, at steady state as \(t\rightarrow \infty \) the solution is known and given by \[ u\left ( x,y,\infty \right ) ={\displaystyle \int \limits _{h/2}^{1-\frac {h}{2}}} {\displaystyle \int \limits _{h/2}^{1-\frac {h}{2}}} u\left ( x,y,0\right ) dxdy \] The following plot shows the solution at \(t=1\) second with the steady state solution displayed as the blue horizontal flat surface superimposed on the same plot. The steady state solution is what would result if run time was made to be very long.
To verify that the numerical solution converges to the steady state solution, the plot below was generate which represents the solution of the above problem taken at \(t=4\) seconds. The gap in the diagram below is the difference between the steady state solution and the solution at \(t=4\) seconds. This gap became smaller the longer the time to run is made (keeping everything else fixed).
Refinement study was carried out to show that the 2D ADI scheme is a second order accurate in time and in space. The method used successive errors between numerical solutions. The algorithm of the refinement study is given below at the end of this part of the problem.
Recalling that In HW1, the spatial grid was divided by half each time. However, in this problem, since cell centered grid is used, \(h\) and \(\Delta t\) were divided by \(3\) each time. This was done so that the new grid will contain some grid locations that are still aligned in the same physical location as the previous step. The error between both solutions is obtained by taking the difference of only these points that are aligned. These points will be the grid point of the coarse grid. The following diagram illustrate this for the case of \(n=3\) and \(n=9\)
The result of the refinement study shows second order accuracy as the error ratio came out to be \(9\).
Below is the result obtained. In addition to the ratio table, it can be seen that the slope of the line in the log plot is \(2\), implying the scheme is second order accurate.
The following is the general outline of the algorithm used in the refinement study. The important part was to make sure when finding the error between the current and last solution, is to use the same physical locations that are aligned between both grids, and to use the coarse grid spacing when determining the grid norm of the error grid.
last_error = 0 h = 1/3 last_h = h delt = h last_u = Solve_2D_ADI(h,delt) LOOP h = h/3 delt = h current_u = Solve_2D_ADI(h,delt) -- now extract from current_u only locations that aligned with last_u grid current_u_mapped = extracted_u(last_u) error = last_h * norm(last_u - current_u_mapped,2) ratio = last_error/error last_h = h last_error = error loop_counter++ IF loop_counter > some_maximum THEN -- normally 5,6 iterations is enough EXIT LOOP END IF END LOOP
The spatial integral represents the total concentration in the domain. Since the boundary are insulated, matter will only diffuse internally and no loss will occur to the outside. Hence, from the conservation of mass principle, initial concentration will remain the same, but will spread out to the mean in space. Therefore, it is known physically total concentration will not change with time\[ \frac {d}{dt}{\displaystyle \int \limits _{\Omega }} u\left ( x,y,t\right ) dA=0 \] The problem is asking to show this mathematically.
Since \(u_{t}=D\Delta u\), then \begin {align*} I & =\frac {d}{dt}{\displaystyle \int \limits _{\Omega }} u\left ( x,y,t\right ) dA\\ & =D{\displaystyle \int \limits _{\Omega }} \Delta u\ dA \end {align*}
But \(\Delta u=\frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\), hence\[ \left ( \frac {1}{D}\right ) I={\displaystyle \int \limits _{y=0}^{1}} {\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\ dxdy \] To show that \(I=0\), the above is written as\begin {align} \left ( \frac {1}{D}\right ) I & ={\displaystyle \int \limits _{y=0}^{1}} {\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial x^{2}}dxdy+{\displaystyle \int \limits _{y=0}^{1}} {\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial y^{2}}\ dxdy\nonumber \\ & ={\displaystyle \int \limits _{y=0}^{1}} \left ( {\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial x^{2}}dx\right ) dy+{\displaystyle \int \limits _{x=0}^{1}} \left ( {\displaystyle \int \limits _{y=0}^{1}} \frac {\partial ^{2}u}{\partial y^{2}}\ dy\right ) dx \tag {1} \end {align}
By applying the fundamental theory of calculus (or using integration by parts) results in\[{\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial x^{2}}dx=\left . \frac {\partial u}{\partial x}\right \rfloor _{x=1}-\left . \frac {\partial u}{\partial x}\right \rfloor _{x=0}\] However \(\left . \frac {\partial u}{\partial x}\right \rfloor _{x=1}\) is the normal derivative at the right boundary, and \(\left . \frac {\partial u}{\partial x}\right \rfloor _{x=0}\) is the normal derivative at the left boundaries. These are both zero due to the homogenous Neumann boundary conditions given in the problem statement. therefore\begin {equation} {\displaystyle \int \limits _{x=0}^{1}} \frac {\partial ^{2}u}{\partial x^{2}}dx=0 \tag {2} \end {equation} Similar argument shows that\begin {equation} {\displaystyle \int \limits _{y=0}^{1}} \frac {\partial ^{2}u}{\partial y^{2}}\ dy=0 \tag {3} \end {equation} Substituting Eqs. (2) and (3) into (1) gives\[ \left ( \frac {1}{D}\right ) I=0 \] Therefore\[ \frac {d}{dt}{\displaystyle \int \limits _{\Omega }} u\left ( x,y,t\right ) dA=0 \]
The finite difference equations for the 2D ADI scheme is given by\begin {align} \left ( I-\frac {D\Delta t}{2}L_{x}\right ) \mathbf {u}^{\ast } & =\left ( I+\frac {D\Delta t}{2}L_{y}\right ) \mathbf {u}^{n}\tag {1}\\ \left ( I-\frac {D\Delta t}{2}L_{y}\right ) \mathbf {u}^{n+1} & =\left ( I+\frac {D\Delta t}{2}L_{x}\right ) \mathbf {u}^{\ast } \tag {2} \end {align}
Summing the equations over all entries in the 2D solution domain gives\begin {align} {\displaystyle \sum \limits _{i}} \sum _{j}\left ( I-\frac {D\Delta t}{2}L_{x}\right ) u^{\ast } & ={\displaystyle \sum \limits _{j}} \sum _{i}\left ( I+\frac {D\Delta t}{2}L_{y}\right ) u^{n}\tag {1A}\\{\displaystyle \sum \limits _{j}} \sum _{i}\left ( I-\frac {D\Delta t}{2}L_{y}\right ) u^{n+1} & ={\displaystyle \sum \limits _{i}} \sum _{j}\left ( I+\frac {D\Delta t}{2}L_{x}\right ) u^{\ast } \tag {2A} \end {align}
In the above \(i\) represents the row number and \(j\) represents the column number of the solution grid \(u\). The above two equations can be rewritten as\begin {align} {\displaystyle \sum \limits _{i}} \sum _{j}u_{ij}^{\ast }-\frac {D\Delta t}{2}{\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast } & ={\displaystyle \sum \limits _{j}} \sum _{i}u_{ij}^{n}+\frac {D\Delta t}{2}{\displaystyle \sum \limits _{j}} \sum _{i}L_{y}u_{ij}^{n}\tag {1B}\\{\displaystyle \sum \limits _{j}} \sum _{i}u_{ij}^{n+1}-\frac {D\Delta t}{2}{\displaystyle \sum \limits _{j}} \sum _{i}L_{y}u_{ij}^{n+1} & ={\displaystyle \sum \limits _{i}} \sum _{j}u_{ij}^{\ast }+\frac {D\Delta t}{2}{\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast } \tag {2B} \end {align}
Looking at the term \({\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast }\) from Eq. (1B) and rewriting this as follows \begin {align*} {\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast } & ={\displaystyle \sum \limits _{i}} \left ( \sum _{j}L_{x}u_{ij}^{\ast }\right ) \\ & ={\displaystyle \sum \limits _{i}} \overset {L_{x}\text { operator applied to }i^{th}\text {row}}{\overbrace {\left ( \sum _{j}L_{x}u_{ij}^{\ast }\right ) }} \end {align*}
In other words, \(\sum _{j}L_{x}u_{ij}^{\ast }\) is the result of applying \(L_{x}\) to each entry in the \(i^{th}\) row, then summing the result.
Therefore, \(L_{x}\) is applied to entry \(u^{\ast }\left ( i,1\right ) \) then to entry \(u^{\ast }\left ( i,2\right ) \) and so on, until the last entry in the row which is \(u^{\ast }\left ( i,n\right ) \).
How to find the result of applying \(L_{x}\) on each row? Given that \(L_{x}\) for 1D with homogenous Neumann boundary conditions is \[ L_{x}=\frac {1}{h^{2}}\begin {bmatrix} -1 & 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ \cdots & \cdots & \cdots & \ddots & \cdots & \cdots \\ 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 & -1 \end {bmatrix} \] Then applying the operator to each entry in the \(i^{th}\) row gives the following\[ \left ( \sum _{j}L_{x}u_{ij}^{\ast }\right ) =\begin {array} [c]{ccccccccccc}j=1 & j=2 & j=3 & 4 & 5 & 6 & 7 & \cdots & n-2 & n-1 & j=n\\ -u_{i1} & +u_{i2} & & & & & & & & & \\ +u_{i1} & -2u_{i2} & +u_{i3} & & & & & & & & \\ & +u_{i2} & -2u_{i3} & +u_{i4} & & & & & & & \\ & & +u_{i3} & -2u_{i4} & +u_{i5} & & & & & & \\ & & & +u_{i4} & -2u_{i5} & +u_{i6} & & & & & \\ & & & & +u_{i5} & -2u_{i6} & +u_{i7} & & & & \\ & & & & & & \ddots & & & & \\ & & & & & & & +u_{i,n-3} & -2u_{i,n-2} & +u_{i,n-1} & \\ & & & & & & & & +u_{i,n-2} & -2u_{i,n-1} & +u_{i,n}\\ & & & & & & & & & +u_{i,n-1} & -u_{i,n}\end {array} \] In the above, \(L_{x}\) was applied directly on the \(i^{th}\) row. The first line above shows the column index \(j\) which goes from \(1\cdots n\). The following diagram is made to help illustrate the above process, showing how \(L_{x}\) and \(L_{y}\) are applied to the solution in the \(u\) matrix.
One can see now that thesum is zero due to terms cancellation. The sum is zero in this case due to the homogenous Neumann boundary conditions which caused the first and last entries to cancel out.
Using the same procedure, then applying \(L_{y}\) to each column of \(u^{n}\) will also result in zero sum, since the north and south boundaries also have homogenous Neumann boundary conditions. Since boundary conditions do not change going from \(u^{\ast }\) to \(u^{n+1}\), the same result is obtained when applying \(L_{y}\) operator to each column of \(u^{n+1}\). From the above discussion, it is found that \begin {align*} {\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast } & =0\\{\displaystyle \sum \limits _{j}} \sum _{i}L_{y}u_{ij}^{n} & =0\\{\displaystyle \sum \limits _{j}} \sum _{i}L_{y}u_{ij}^{n+1} & =0\\{\displaystyle \sum \limits _{i}} \sum _{j}L_{x}u_{ij}^{\ast } & =0 \end {align*}
Substituting the above 4 equations back into Eqs. (1B),(2B) gives\begin {align} {\displaystyle \sum \limits _{i}} \sum _{j}u_{ij}^{\ast } & ={\displaystyle \sum \limits _{j}} \sum _{i}u_{ij}^{n}\tag {1C}\\{\displaystyle \sum \limits _{j}} \sum _{i}u_{ij}^{n+1} & ={\displaystyle \sum \limits _{i}} \sum _{j}u_{ij}^{\ast } \tag {2C} \end {align}
Substituting Eq. (1C) into (2C) gives\[ \sum _iju_ij^n+1=\sum _iju_ij^n \] Since the above is valid for any \(n\) (boundary conditions do not change with time) then setting \(n=0\) in the above results in\[ \sum _{ij}u_{ij}^{1}=\sum _{ij}u_{ij}^{0}\] Similarly, setting \(n=1\) results in\[ \sum _{ij}u_{ij}^{2}=\sum _{ij}u_{ij}^{1}\] and so on all the way any \(n\) value. Hence in general the following result is obtained\[ \sum _{ij}u_{ij}^{n}=\sum _{ij}u_{ij}^{n-1}\] By repeated back substitution on the RHS, the following is obtained\[ \sum _{ij}u_{ij}^{n}=\sum _{ij}u_{ij}^{0}\] Therefore, the discrete conservation property is satisfied.
Verification in code To verify part(d) in the code, a table was generated during one run, where \(\sum _{ij}u_{ij}^{n}\) was calculated at the end of each time step using the Matlab command sum(sum(u)), and this value was printed at each time step. The result shows that this value is constant implying the discrete conservation property is satisfied. Here is the result below
current_time sum(U(current_time) 0.00000 1897.85094 0.01235 1897.85094 0.02469 1897.85094 0.03704 1897.85094 0.04938 1897.85094 0.06173 1897.85094 0.07407 1897.85094 0.08642 1897.85094 0.09877 1897.85094 0.11111 1897.85094 .... 0.92593 1897.85094 0.93827 1897.85094 0.95062 1897.85094 0.96296 1897.85094 0.97531 1897.85094 0.98765 1897.85094
The equations to solve are the following\begin {align*} \frac {\partial v}{\partial t} & =D\Delta v+\left ( a-v\right ) \left ( v-1\right ) v-w+I\\ \frac {\partial w}{\partial t} & =\epsilon \left ( v-\gamma w\right ) \end {align*}
The first PDE \(\frac {\partial v}{\partial t}=D\Delta v+\left ( a-v\right ) \left ( v-1\right ) v-w+I\) was solved by the splitting method by solving the diffusion equation\(\frac {\partial v}{\partial t}=D\Delta v\) using ADI method separately and then by solving the reaction (non-linear) equation \(\frac {\partial v}{\partial t}=\left ( a-v\right ) \left ( v-1\right ) v-w+I\) along with \(\frac {\partial w}{\partial t}=\epsilon \left ( v-\gamma w\right ) \) separately. The following is a coupled first order non-linear differential equations system (the reaction ODE is nonlinear in votage \(v\))
\begin {align*} \frac {dv}{dt} & =\left ( a-v\right ) \left ( v-1\right ) v-w+I\\ \frac {dw}{dt} & =\epsilon \left ( v-\gamma w\right ) \end {align*}
The above system was solved using Runge-Kutta order 4. The following diagram illustrates the time line for one full splitting step.
ADI was described in problem 2, and the same function was reused for this problem for the diffusion solver. For solving the reaction system of equations, RK4 was implemented as follows. define\[ f\left ( v,w\right ) =\left ( a-v\right ) \left ( v-1\right ) v-w+I \] and also define \[ g\left ( v,w\right ) =\epsilon \left ( v-\gamma w\right ) \] Therefore, the RK4 solver for the above system becomes
\begin {align*} v^{n+1} & =v^{n}+\frac {1}{6}\left ( m_{1}+2m_{2}+2m_{3}+m_{4}\right ) \\ w^{n+1} & =w^{n}+\frac {1}{6}\left ( k_{1}+2k_{2}+2k_{3}+k_{4}\right ) \end {align*}
Where
\begin {align*} m_{1} & =\Delta tf\left ( v,w\right ) \\ m_{2} & =\Delta tf\left ( v+\frac {1}{2}m_{1},w+\frac {1}{2}k_{1}\right ) \\ m_{3} & =\Delta tf\left ( v+\frac {1}{2}m_{2},w+\frac {1}{2}k_{2}\right ) \\ m_{4} & =\Delta tf\left ( v+m_{3},w+k_{3}\right ) \end {align*}
And
\begin {align*} k_{1} & =\Delta tg\left ( v,w\right ) \\ k_{2} & =\Delta tg\left ( v+\frac {1}{2}m_{1},w+\frac {1}{2}k_{1}\right ) \\ k_{3} & =\Delta tg\left ( v+\frac {1}{2}m_{2},w+\frac {1}{2}k_{2}\right ) \\ k_{4} & =\Delta tg\left ( v+m_{3},w+k_{3}\right ) \end {align*}
Another point regarding the splitting method. It was required to decide which splitting method to use. Should a simple splitting, Strang splitting or the 2-step splitting method which was described in class be used? To make sure the second order accuracy of ADI 2D in time is preserved, simple splitting was not used (unless the operators commute, this would have caused the scheme to become first order accurate in time). Instead, the two step splitting method was used, as it was found to be simpler than Strang method to implement.
The program written in part(a) was run using the parameters given. The time step used was set to be the same as the space step. This time step is recommended for the ADI The diffusion solver as it is a second order accurate in time and space. This is the fast system (the stiff part of the system), hence making the time step larger than the space step would not give accurate results, even though it will remain a stable scheme. Keeping the time step the same as the space step seemed to be a good choice, as it kept the time resolution and the space resolution the same. The same time step was then used for the reaction solver, as was required by the splitting method to keep each step the same length.
The following shows the visualization of the voltage solution for up to 300 seconds as required using the surf() command.
The solution \(v\left ( t\right ) \) started from a peak value at one corner of the square. Shortly after, at about \(50\) seconds, a wave started to form, the wave front became large and it spread out and advanced with time. When the pulse reached the boundary on the other corner, it started to diffuse and by \(t=300\) seconds, the pulse has completely disappeared.
The following is the result of the simulation for this part.
In this simulation, the pulse that appeared after about \(50\) second, quickly became a spiral, and did not appear to diffuse as was the case in part (d). At the end of the simulation, the pulse was continuing to spiral in the same rotation direction it started with. The above phenomena seem to be termed an arrhythmia pulse.
One common theme between part(c) and part (b), is the formation of a wave like motion that traveles across the domain. The difference was in the shape of the pulse, the direction it moves to and the amount of diffusion that occured.
Derivation of ADI equations for 2D and 3D for the diffusion problem Given \(u_{t}=\Delta u\) (\(D\) is assumed 1), then in 2D forward Euler (explicit) gives\[ \frac {u_{ij}^{n+1}-u_{ij}^{n}}{k}=\left ( \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n}\] While C-N method gives\[ \frac {u_{ij}^{n+1}-u_{ij}^{n}}{k}=\frac {1}{2}\left [ \left ( \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n+1}+\left ( \frac {\partial ^{2}u}{\partial x^{2}}+\frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n}\right ] \] However, in ADI, the time step itself is divided by half, and in the first half step, one of the spatial second derivatives is implicit while and the other spatial second derivative is explicit. In the second half step, these are reversed. \begin {align*} \frac {u_{ij}^{n+\frac {1}{2}}-u_{ij}^{n}}{k/2} & =\overset {implicit}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial x^{2}}\right ) _{ij}^{n+\frac {1}{2}}}}+\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n}}}\\ \frac {u_{ij}^{n+1}-u_{ij}^{n+\frac {1}{2}}}{k/2} & =\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial x^{2}}\right ) _{ij}^{n+\frac {1}{2}}}}+\overset {implicit}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n+1}}} \end {align*}
Writing \(\frac {\partial ^{2}u}{\partial x^{2}}=L_{x}=\frac {u_{i-1,j}-2u_{ij}+u_{i+1,j}}{h^{2}}\) and \(\frac {\partial ^{2}u}{\partial y^{2}}=L_{y}=\frac {u_{i,j-1}-2u_{ij}+u_{i,j+1}}{h^{2}}\), the above 2 equations become\begin {align*} \frac {u_{ij}^{n+\frac {1}{2}}-u_{ij}^{n}}{k/2} & =\frac {u_{i-1,j}^{n+\frac {1}{2}}-2u_{ij}^{n+\frac {1}{2}}+u_{i+1,j}^{n+\frac {1}{2}}}{h^{2}}+\frac {u_{i,j-1}^{n}-2u_{ij}^{n}+u_{i,j+1}^{n}}{h^{2}}\\ \frac {u_{ij}^{n+1}-u_{ij}^{n+\frac {1}{2}}}{k/2} & =\frac {u_{i-1,j}^{n+\frac {1}{2}}-2u_{ij}^{n+\frac {1}{2}}+u_{i+1,j}^{n+\frac {1}{2}}}{h^{2}}+\frac {u_{i,j-1}^{n+1}-2u_{ij}^{n+1}+u_{i,j+1}^{n+1}}{h^{2}} \end {align*}
Moving all implicit terms to the LHS and rearranging results in\begin {align*} u_{ij}^{n+\frac {1}{2}}-\frac {k}{2}\frac {u_{i-1,j}^{n+\frac {1}{2}}-2u_{ij}^{n+\frac {1}{2}}+u_{i+1,j}^{n+\frac {1}{2}}}{h^{2}} & =u_{ij}^{n}+\frac {k}{2}\frac {u_{i,j-1}^{n}-2u_{ij}^{n}+u_{i,j+1}^{n}}{h^{2}}\\ u_{ij}^{n+1}-\frac {k}{2}\frac {u_{i,j-1}^{n+1}-2u_{ij}^{n+1}+u_{i,j+1}^{n+1}}{h^{2}} & =u_{ij}^{n+\frac {1}{2}}+\frac {k}{2}\frac {u_{i-1,j}^{n+\frac {1}{2}}-2u_{ij}^{n+\frac {1}{2}}+u_{i+1,j}^{n+\frac {1}{2}}}{h^{2}} \end {align*}
Hence, in operator form the above becomes\begin {align*} \left ( I-\frac {k}{2}L_{x}\right ) u_{ij}^{n+\frac {1}{2}} & =\left ( I+\frac {k}{2}L_{y}\right ) u_{ij}^{n}\\ \left ( I-\frac {k}{2}L_{y}\right ) u_{ij}^{n+1} & =\left ( I+\frac {k}{2}L_{x}\right ) u_{ij}^{n+\frac {1}{2}} \end {align*}
In the class notes, \(u_{ij}^{\ast }\) was used to represent \(u_{ij}^{n+\frac {1}{2}}\) but they are the same. The above is the ADI scheme for 2D. The 3D equations are now derived. Since three different directions exist now, the time step is divided into 3. This results in\begin {align*} \frac {u_{ij}^{n+\frac {1}{3}}-u_{ij}^{n}}{\Delta t/3} & =\overset {implicit}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial x^{2}}\right ) _{ij}^{n+\frac {1}{3}}}}+\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n}}}+\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial z^{2}}\right ) _{ij}^{n}}}\\ \frac {u_{ij}^{n+\frac {2}{3}}-u_{ij}^{n+\frac {1}{3}}}{\Delta t/3} & =\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial x^{2}}\right ) _{ij}^{n+\frac {1}{3}}}}+\overset {implicit}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n+\frac {2}{3}}}}+\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial z^{2}}\right ) _{ij}^{n+\frac {1}{3}}}}\\ \frac {u_{ij}^{n+1}-u_{ij}^{n+\frac {2}{3}}}{\Delta t/3} & =\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial x^{2}}\right ) _{ij}^{n+\frac {2}{3}}}}+\overset {\text {explicit}}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial y^{2}}\right ) _{ij}^{n+\frac {2}{3}}}}+\overset {implicit}{\overbrace {\left ( \frac {\partial ^{2}u}{\partial z^{2}}\right ) _{ij}^{n+1}}} \end {align*}
Similar to what done in the 2D case, the above are rearranged resulting in\begin {align} \left ( I-\frac {D\Delta t}{3}L_{x}\right ) \mathbf {u}^{n+\frac {1}{3}} & =\left ( I+\frac {D\Delta t}{3}L_{y}+\frac {D\Delta t}{3}L_{z}\right ) \mathbf {u}^{n}\tag {1}\\ \left ( I-\frac {D\Delta t}{3}L_{y}\right ) \mathbf {u}^{n+\frac {2}{3}} & =\left ( I+\frac {D\Delta t}{3}L_{x}+\frac {D\Delta t}{3}L_{z}\right ) \mathbf {u}^{n+\frac {1}{3}}\tag {2}\\ \left ( I-\frac {D\Delta t}{3}L_{z}\right ) \mathbf {u}^{n+1} & =\left ( I+\frac {D\Delta t}{3}L_{x}+\frac {D\Delta t}{3}L_{y}\right ) \mathbf {u}^{n+\frac {2}{3}} \tag {3} \end {align}
6Please see the appendix of this problem at the end of the HW report showing how these equations came about.
7To show these operators commute, similar argument can be made as was done for the 2D case in class, which is by saying that each operator \(L_{x},L_{y},L_{z}\) on its own commutes with the other 2, hence the result will follow.