The goal of a refinement study is to perform a numerical experiment to determine the order of accuracy of a given finite difference scheme. The appendix of this problem contain a review of the idea behind refinement study.
The problem asked us to determine the order of accuracy in time and in space. A program implementing the above scheme was run a number of times, each time with a different initial value for the space and time step. To verify order of accuracy for the C-N scheme, the space and the time step were divided by 2 simultaneously before the start of each run. To verify order of accuracy for the forward Euler scheme, the space step was divided by 2 but the time was divided by 4. For both schemes , the program generated ratios of successive errors between the numerical solutions at the end of each run (1 second long run).
Convergence of this ratio to the value 4 implied the results we are asked to demonstrate.
In the following, the C-N and the forward Euler finite difference schemes are derived, then the numerical results presented, followed by a conclusion.
The method of lines (MOL) was used to implement the C-N scheme to solve for the numerical solution \(u.\) The equations are solved using Matlab’s \(u=A\backslash b\) where \(A\) is a sparse matrix (the system update matrix) constructed based on the C-N discretization. An efficient algorithm to solve for \(u\) in this scheme is Thomas algorithm version of Gaussian elimination. It is understood that this will automatically be done by Matlab ”\(\backslash \)” operator when it recognizes that the \(A\) matrix is a tridiagonal giving an \(O\left ( n\right ) \) order for the solver where \(n\) is the number of unknowns.
Let the PDE be\begin {equation} d\ u_{t}-Du_{xx}+au=g\left ( x,t\right ) \tag {1} \end {equation} \(g\left ( x,t\right ) \) is an internal source with initial conditions as \(u\left ( x,t\right ) =u0\left ( x\right ) .\) The Dirichlet boundary conditions are\[ \left \{ \begin {array} [c]{c}u\left ( 0,t\right ) =\alpha \left ( t\right ) \\ u\left ( L,t\right ) =\beta \left ( t\right ) \end {array} \right . \] and Neumann boundary conditions are\[ \left \{ \begin {array} [c]{c}u_{t}\left ( 0,t\right ) =\alpha \left ( t\right ) \\ u_{t}\left ( L,t\right ) =\beta \left ( t\right ) \end {array} \right . \] The terms \(d,a\) above are constants, and \(D\) is the diffusion constant. For the C-N scheme (1) was discretized at point \(x_{j}\) with space step as \(h\) and with time step as \(k\) resulting in\[ d\ \frac {u_{j}^{n+1}-u_{j}^{n}}{k}=\frac {1}{2}\left ( f_{j}^{n}+f_{j}^{n+1}\right ) \] Where \(f^{n}\) is the RHS of the PDE at time \(t_{n}=nk\), so the above becomes\begin {align*} d\ \frac {u_{j}^{n+1}-u_{j}^{n}}{k} & =\frac {1}{2}\left [ \left ( Du_{xx}-au+g_{j}\right ) ^{n}+\left ( Du_{xx}-au+g_{j}\right ) ^{n+1}\right ] \\ & \\ & =\frac {1}{2}\left ( D\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}-au_{j}^{n}+g_{j}^{n}\right ) +\\ & \left ( D\frac {u_{j-1}^{n+1}-2u_{j}^{n+1}+u_{j+1}^{n+1}}{h^{2}}-au^{n+1}+g_{j}^{n+1}\right ) \\ & \\ & =\frac {D}{2h^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) +\frac {D}{2h^{2}}\left ( u_{j-1}^{n+1}-2u_{j}^{n+1}+u_{j+1}^{n+1}\right ) -\\ & \frac {a}{2}\left ( u_{j}^{n}+u_{j}^{n+1}\right ) +\frac {1}{2}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \end {align*}
collecting all terms at time \(n+1\) to the left gives\[ u_{j}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) -\frac {kD}{2dh^{2}}u_{j-1}^{n+1}-\frac {kD}{2dh^{2}}u_{j+1}^{n+1}=u_{j}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{j-1}^{n}\frac {kD}{2dh^{2}}+u_{j+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \] Let \begin {align*} r_{1} & =\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) \\ r_{2} & =\frac {kD}{2dh^{2}}\\ r_{3} & =\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) \\ r_{4} & =\frac {k}{2d} \end {align*}
Then the above becomes\begin {equation} r_{1}u_{j}^{n+1}-r_{2}u_{j-1}^{n+1}-r_{2}u_{j+1}^{n+1}=r_{3}u_{j}^{n}+r_{2}u_{j-1}^{n}+r_{2}u_{j+1}^{n}+r_{4}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \tag {2} \end {equation} The above algebraic equation (2) is the C-N finite difference scheme for (1) and is valid for \(x_{j}\) at the internal points. Considering the case of both ends having Dirichlet boundary conditions, and using the following grid numbering1
Then (2) above is valid at the internal nodes numbered \(j=2\cdots N-1\). Hence \(u_{1}^{n}\) will be the left boundary point and \(u_{N}^{n}\) will be the right boundary point. When the boundary conditions are Dirichlet, let \(u_{1}^{n}=\alpha \left ( t_{n}\right ) \) and \(u_{N}^{n}=\beta \left ( t_{n}\right ) \). Converting (2) to matrix form results in\[ \overset {\mathbf {A}}{\overbrace {\begin {bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{1} & -r_{2} & 0 & 0 & 0 & 0\\ 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0 & 0\\ 0 & 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \ddots & \vdots & 0\\ 0 & 0 & 0 & 0 & -r_{2} & r_{1} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end {bmatrix} }}\overset {\mathbf {x}}{\overbrace {\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-2}^{n+1}\\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} }}=\overset {\mathbf {b}}{\overbrace {\begin {bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{3} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{3} & r_{2} & 0 & 0 & 0\\ 0 & 0 & r_{2} & r_{3} & r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & 0\\ 0 & 0 & 0 & 0 & r_{2} & r_{3} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} \alpha ^{n+1}\\ r_{2}\tilde {\alpha }+r_{4}\tilde {g}_{2}\\ r_{4}\tilde {g}_{3}\\ r_{4}\tilde {g}_{4}\\ \vdots \\ r_{2}\tilde {\beta }+r_{4}\tilde {g}_{N-1}\\ \beta ^{n+1}\end {bmatrix} }}\] Where in (3) \(\alpha \left ( t_{n+1}\right ) +\alpha \left ( t_{n}\right ) \equiv \tilde {\alpha }\) and \(\beta \left ( t_{n+1}\right ) +\beta \left ( t_{n}\right ) \equiv \tilde {\beta }\) and \(g^{n}+g^{n+1}\equiv \tilde {g}\) .
Equation (3) is in the form \(Au^{n+1}=b\). And \(u^{n+1}\) is solved for using \(A\backslash b\). Notice that (3) is in the same form shown in class notes, which is \begin {equation} \left ( I-\frac {Dk}{2}L\right ) u^{n+1}=\left ( I+\frac {Dk}{2}L\right ) u^{n}+kf^{n+\frac {1}{2}} \tag {4} \end {equation} Where the update matrix \(B=\left ( I-\frac {Dk}{2}L\right ) ^{-1}\left ( I+\frac {Dk}{2}L\right ) \). \(L\) is the standard Laplace operator for \(1D\) problem given by\[\begin {bmatrix} -2 & 1 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 1 & -2 \end {bmatrix} \] Notice that (3) compared to (4), has additional terms included in the RHS in order to support the general form the parabolic PDE. Equation (4) represents the diffusion pde \(u_{t}-Du_{xx}=0.\)
\(u_{j}^{n}\) when \(n=0\) is obtained from initial conditions. The first step solves for \(u_{j}^{1}\) which is then used in the second second step to solve for \(u_{j}^{2}\) and so on, until the maximum time to solve for is reached. Since Dirichlet boundary conditions are used, \(u_{1}^{n}\) (the solution at the left edge) and \(u_{N}^{n}\), the solution at the right edge are always known. The above system is solved only for the internal nodes. Next section shows the numerical results.
Result for part(a) The above scheme was implemented with a GUI added to make it easier to use these algorithms. The following plot shows the numerical solution at \(t=1\).
The following is the ratio error table. This table shows that the ratio converged to 4.
# delt h ratio 2 0.2500000 0.2500000 1.0000e+000 3 0.1250000 0.1250000 2.9347e+000 4 0.0625000 0.0625000 3.6798e+000 5 0.0312500 0.0312500 4.2132e+000 6 0.0156250 0.0156250 4.1209e+000 7 0.0078125 0.0078125 4.0354e+000 8 0.0039063 0.0039063 4.0092e+000
The following is the loglog plot of the above result. The x-axis represents \(h\) and the y-axis the difference in errors (absolute). The slope of the line is seen to be 2 implying a second order accuracy.
Conclusion Since the ratio is \(4\) and since the time step and the space step were halved in each run, this implies C-N is second order in time and space.
Starting with the same PDE as in part (a)\begin {equation} d\ u_{t}-Du_{xx}+au=g\left ( x,t\right ) \tag {1} \end {equation} All terms and boundary conditions and the solution domain are as shown in part (a).
For the forward Euler scheme (1) was discretized at point \(x_{j}\) with space step as \(h\) and with time step as \(k\) as follows\begin {align*} d\ \frac {u_{j}^{n+1}-u_{j}^{n}}{k} & =f_{j}^{n}\\ & =D\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}-au_{j}^{n}+g_{j}^{n} \end {align*}
Hence\begin {align*} u_{j}^{n+1} & =u_{j}^{n}+\frac {kD}{dh^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) -\frac {ak}{d}u_{j}^{n}+\frac {k}{d}g_{j}^{n}\\ & =\frac {kD}{dh^{2}}u_{j-1}^{n}+u_{j}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{j+1}^{n}+\frac {k}{d}g_{j}^{n} \end {align*}
Let \(r_{1}=\left ( 1-2\frac {kD}{dh^{2}}-\frac {ak}{d}\right ) ,r_{2}=\frac {kD}{dh^{2}},\),\(r_{3}=\) \(\frac {k}{d}\), the above becomes
\begin {equation} u_{j}^{n+1}=r_{2}u_{j-1}^{n}+r_{1}u_{j}^{n}+r_{2}u_{j+1}^{n}+r_{3}g_{j}^{n} \tag {2} \end {equation}
The above algebraic equation (2) is the forward Euler finite difference scheme for (1) and is valid for \(x_{j}\) at the internal points.
Therefore, the stencil for the forward Euler scheme for the 1D parabolic PDE is
Considering the case of both ends having Dirichlet boundary conditions, and using the same numbering as in part (a) then (2) is valid at the internal nodes \(j=2\cdots N-1\).
\(u_{1}^{n}\) will be the left boundary point and \(u_{N}^{n}\) will be the right boundary point. Let \(u_{1}^{n}=\alpha \left ( t_{n}\right ) \) and \(u_{N}^{n}=\beta \left ( t_{n}\right ) \). Converting (2) to matrix form results in\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ \vdots \\ u_{N-2}^{n+1}\\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} =\begin {bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{1} & r_{2} & 0 & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{1} & r_{2} & 0 & 0 & 0 & 0\\ 0 & 0 & r_{2} & r_{1} & r_{2} & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \ddots & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & r_{2} & r_{1} & r_{2} & 0\\ 0 & 0 & 0 & 0 & 0 & r_{2} & r_{1} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N-2}^{n}\\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} \alpha ^{n+1}\\ r_{2}\alpha ^{n}+r_{3}g_{2}^{n}\\ r_{3}g_{3}^{n}\\ r_{3}g_{4}^{n}\\ \vdots \\ r_{3}g_{N-2}^{n}\\ r_{2}\beta ^{n}+r_{3}g_{N-1}^{n}\\ \beta ^{n+1}\end {bmatrix} \] And now \(u^{n+1}\) is found by direct matrix/vector multiplication as shown above. No matrix inversion is required in this case since this is an explicit method.
Looking at the stencil above, an idea is now suggested to determine stability directly from the stencil diagram. By imposing that the weight on each edge in the directed graph not exceed unity, and that the total algebraic sum of the weight of the edge also not exceed unity. This includes any combination of edges involved. If this is always the case, then \(u_{j}^{n+1}\) will always have an amplitude \(\leq \) \(u_{j}^{n}\) since the weights are never more than \(1\) no matter what combinations are used. This idea is applied to this problem with \(u_{t}+Du_{xx}=0\,\), hence \(a=0\) and \(d=1\). This gives that following conditions on the edges shown in the stencil diagram above\[ \left \{ \begin {array} [c]{lll}\left ( 1\right ) & \frac {kD}{h^{2}}\leq 1 & \text {Condition on }j-1\text { or }j+1\text { separately}\\ \left ( 2\right ) & 2\frac {kD}{h^{2}}\leq 1 & \text {Condition on }j-1\text { and }j+1\text { added together}\\ \left ( 3\right ) & \left \vert 1-2\frac {k\ D}{\ h^{2}}\right \vert \leq 1 & \text {Condition on the }j\text { edge}\\ \left ( 4\right ) & \left \vert \frac {kD}{h^{2}}+1-2\frac {k\ D}{\ h^{2}}\right \vert \leq 1 & \text {Condition on the }j\text { edge with either }j-1\text { or }j+1\\ \left ( 5\right ) & \left \vert 2\frac {kD}{h^{2}}+1-2\frac {k\ D}{h^{2}}\right \vert \leq 1 & \text {Condition that all edges sum to less than 1}\end {array} \right . \] Condition (1) is weaker than (2), hence not considered. Condition (3) results in \(\frac {kD}{h^{2}}\leq 1\) which is the same as (1). Condition (4) gives \(\left \vert 1-\frac {k\ D}{\ h^{2}}\right \vert \leq 1\) or \(\frac {k\ D}{\ h^{2}}\leq 2\) which is also weaker than (2). Condition (5) gives \(1\leq 1\) hence no information is obtained from it. Therefore, condition (2) remains, and that condition says that \(\frac {kD}{h^{2}}\leq \frac {1}{2}\), which is the strongest condition. Hence, this is the absolute condition for stability for forward Euler. This agrees with the method to determine this using Von Neumann analysis.
Result for part(b) The forward Euler results are below. The space step was divided by 2 and the time step was divided by 4.
# delt h ratio 2 0.0625000 0.2500000 1.0000e+000 3 0.0156250 0.1250000 3.0842e+000 4 0.0039063 0.0625000 3.7567e+000 5 0.0009766 0.0312500 4.1506e+000 6 0.0002441 0.0156250 4.0794e+000 7 0.0000610 0.0078125 4.0447e+000 8 0.0000153 0.0039063 3.9909e+000
The following is the corresponding loglog plot
Conclusion Since the ratio is \(4\) and since the time step was divided by 4 and the space step by 2, this implies forward Euler is first order in time and second order in space.
The appendix of this problem show the steady state analytical solution to the above PDE derived using Laplace transform method.
Review of refinement study process The idea behind refinement study is reviewed briefly. Assume the goal is to find the order of accuracy of a finite difference scheme with respect to the space step. The finite difference formula is first derived, and the exact solution is substituted into this formula. Terms that contain \(u\left ( x\pm h\right ) \), are replaced by Taylor series approximation. The result is simplified, and the error term is found. An small example is given to illustrate the idea.
To find the order of accuracy in space using forward Euler finite difference approximation to a derivative\[ u^{\prime }\left ( x\right ) =\frac {u\left ( x+h\right ) -u\left ( x\right ) }{h}\] Each term for \(u\) in the RHS above is replaced by its exact value, using Taylor series expansion where needed, resulting in\begin {align} u^{\prime }\left ( x\right ) & =\frac {\left [ u\left ( x\right ) +hu^{\prime }\left ( x\right ) +\frac {h^{2}}{2}u^{\prime \prime }\left ( x\right ) +\frac {h^{3}}{3!}u^{\prime \prime \prime }\left ( x\right ) +\cdots \right ] -u\left ( x\right ) }{h}\nonumber \\ & =u^{\prime }\left ( x\right ) +\overset {error}{\overbrace {\frac {h}{2}u^{\prime \prime }\left ( x\right ) +\frac {h^{2}}{3!}u^{\prime \prime \prime }\left ( x\right ) +\cdots }} \tag {1} \end {align}
The error term is hence. It is the amount that the RHS differs from the LHS. The leading error term in (1) (the dominant term) is \(\frac {h}{2}u^{\prime \prime }\left ( x\right ) \), but since \(x\) is a known value (the above is being evaluated at each grid point, hence \(x\) is known), then \(u^{\prime \prime }\left ( x\right ) \) is some constant, and the leading error term in (1) is of the form \(Ch\), where \(C\) is some constant. This is the same as saying that the error is of order \(h\).
The above method can be used to find the order of the error in approximation when the exact solution is know. In problem (1), the exact solution is not given and was difficult to obtain. Hence, instead of finding the order of accuracy using the above method, it was found using a numerical experiment (refinement study).
In the refinement study the error itself is determined, and from the error profile (as \(h\) is changed), the order is determined. But this error is the error between successive numerical solutions.
Once the numerical error is found (after running the refinement study), then one method to find \(p\) (order of the error ) is to take the logarithm resulting in \begin {align*} error & =Ch^{p}\\ \log \left ( error\right ) & =p\log \left ( h\right ) +\log C\\ & =p\log \left ( h\right ) +\text {constant} \end {align*}
and this represents an equation of the line \(Y=pX+k\), where \(p\) is the line slope which is the same as the order of accuracy. Hence, by generating different \(h\) values, and for each \(h\) determine the corresponding error, then \(p\) is found by measuring at the slope of line from the plot generated. If the slope is \(p=1\), then it is first order accuracy, and if the slope is \(p=2\), it is second order.
The above is a graphical method. Another method is as follows: Starting with some \(h\) value, the error \(e_{n-1}\) is found, then \(h\) is divided by half and the error, now called \(e_{n}\) is found again. The ratio \(\frac {e_{n-1}}{e_{n}}\) is found. If \(p\) happened to be \(2\), then the ratio will come out to be \(4\). This is because \(\frac {e_{n-1}}{e_{n}}=\frac {\left ( h_{n-1}\right ) ^{p}}{\left ( \frac {h_{n-1}}{2}\right ) ^{p}}=2^{p}\) and so if \(p=2\), then the ratio will be \(4\). If \(p=1\), then the ratio will be \(2\).
In the above description, errors are found using differences between successive solutions as follows \begin {align*} e_{n-1} & =\left \vert U_{n+1}-U_{n}\right \vert \\ e_{n} & =\left \vert U_{n}-U_{n-1}\right \vert \end {align*}
The norm used to measure \(U\), the approximate solution, is the Euclidean norm modified for the space grid\[ \left \Vert U\right \Vert =\sqrt {h}\left \Vert U\right \Vert _{2}\]
Steady state analytical solution to the PDE The following shows the steps used to determine the steady state solution for \begin {equation} u_{t}=au_{xx}+f\left ( t\right ) \tag {1} \end {equation} where \(a=\frac {1}{100}\) and \(f\left ( t\right ) =1-e^{-t}\) with initial conditions \(u\left ( 0,t\right ) =0\) and boundary conditions \(u\left ( 0,t\right ) =0\) and \(u\left ( 1,t\right ) =0.\)
The above is an inhomogeneous PDE (the source term \(1-e^{-t}\)). The boundary conditions are homogeneous, and with zero initial conditions.
Since this is an inhomogeneous PDE, separation of variables can not be used. But the steady state solution (the particular solution) can be found using an integral transform approach. Integral transformation is first applied to the PDE, resulting in an ODE which is then solved in the new transformed space, and the solution in time domain is found by inverse transforming back.
Since the spatial domain in this problem is a bounded interval (from 0 to 1), Fourier transformation will not be used because the spatial domain is bounded and does not match the Fourier transformation domain (from \(-\infty \) to \(+\infty \)), however, Laplace transformation (for \(t>0\)) can be used as it matches the time domain of the problem.
Therefore, taking the Laplace transform of (1) w.r.t time gives\[ -u\left ( x,0\right ) +sU\left ( x,s\right ) =a\frac {d^{2}U\left ( x,s\right ) }{dx^{2}}+\frac {1}{s}-\frac {1}{1+s}\] But \(u\left ( x,0\right ) =0\), hence the resulting ODE is\[ a\frac {d^{2}U\left ( x,s\right ) }{dx^{2}}-sU\left ( x,s\right ) =\frac {1}{1+s}-\frac {1}{s}\] With the boundary conditions \(U\left ( 0,s\right ) =0\) and \(U\left ( 1,s\right ) =0\) obtained from the spatial domain. The above ODE is a second order, linear ODE, a inhomogeneous ODE that can be solved for \(U\left ( x,s\right ) \), which results in the following (for the case \(a=\frac {1}{100})\)\[ U\left ( x,s\right ) =\frac {-e^{-10x\sqrt {s}}\left ( e^{10x\sqrt {s}-1}-1\right ) \left ( e^{10x\sqrt {s}}-e^{10\sqrt {s}}\right ) }{s^{2}\left ( 1+s\right ) \left ( 1+e^{10\sqrt {s}}\right ) }\] The steady state solution can now be found using the limit theorem for Laplace transform, giving\begin {align} u\left ( x,\infty \right ) & =\lim _{s->0}sU\left ( x,s\right ) \tag {1}\\ & =50x\left ( 1-x\right ) \end {align}
Here is a plot of the particular solution
x=0:0.01:1; plot(x,50*x.*(1-x))
Derivation of forward Euler for periodic boundary conditions From part(b) above,
\[ u_{j}^{n+1}=\frac {kD}{dh^{2}}u_{j-1}^{n}+u_{j}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{j+1}^{n}+\frac {k}{d}g_{j}^{n}\]
Periodic boundary conditions implies \(u\left ( 0,t\right ) =u\left ( 1,t\right ) \), Hence \(u_{j-1}\) when \(j\) is the first note on the left is the same as node \(N-1\). And \(u_{j+1}\) when \(j\) is the last node on the right is the same as node \(j=2\). As shown in the diagram below
Then \[ u_{1}^{n+1}=\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n}\]
And
\[ u_{N}^{n+1}=\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n}\]
Let \(r_{1}=\left ( 1-2\frac {kD}{dh^{2}}-\frac {ak}{d}\right ) ,r_{2}=\frac {kD}{dh^{2}},r_{3}=\) \(\frac {k}{d}\), Hence the system can be written as
\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} =\begin {bmatrix} r_{1} & r_{2} & 0 & 0 & 0 & r_{2} & 0\\ r_{2} & r_{1} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{1} & r_{2} & 0 & 0 & 0\\ 0 & 0 & r_{2} & r_{1} & r_{2} & 0 & 0\\ 0 & 0 & 0 & 0 & \ddots & 0 & 0\\ 0 & 0 & 0 & 0 & r_{2} & r_{1} & r_{2}\\ 0 & r_{2} & 0 & 0 & 0 & r_{2} & r_{1}\end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} 0\\ r_{3}g_{2}^{n}\\ r_{3}g_{3}^{n}\\ r_{3}g_{4}^{n}\\ \vdots \\ r_{3}g_{N-1}^{n}\\ 0 \end {bmatrix} \]
Derivation of forward Euler for Neumann boundary conditions both ends Using this numbering
Assume that \(u_{t}=\alpha \) at node \(1\) and \(u_{t}=\beta \) at node \(N\) (these are the Neumann boundary conditions).
Add a ghost node \(0\) to the left of node \(1\), and approximating \(\alpha \left ( t\right ) \) gives\[ \alpha =\frac {u_{0}-u_{2}}{2h}\] hence\begin {equation} u_{0}=2h\alpha +u_{2} \tag {1} \end {equation} But the PDE for node 1 is\begin {equation} u_{1}^{n+1}=\frac {kD}{dh^{2}}u_{0}^{n}+u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n} \tag {2} \end {equation}
Substitute (1) into (2) gives
\begin {align*} u_{1}^{n+1} & =\frac {kD}{dh^{2}}\left ( 2h\alpha +u_{2}\right ) +u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n}\\ & =u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +u_{2}^{n}\left ( 2\frac {kD}{dh^{2}}\right ) +2h\alpha \frac {kD}{dh^{2}} \end {align*}
Similarly for the right end. Add a ghost node \(N+1\) to the right of node \(N\), and approximating \(\beta \left ( t\right ) \) gives
\[ \beta =\frac {u_{N-1}-u_{N+1}}{2h}\]
hence
\begin {equation} u_{N+1}=2h\beta +u_{N-1} \tag {3} \end {equation}
But the PDE for node N is
\begin {equation} u_{N}^{n+1}=\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{N+1}^{n} \tag {4} \end {equation}
Substitute (3) into (4) gives
\begin {align*} u_{1}^{n+1} & =\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}\left ( 2h\beta +u_{N-1}\right ) \\ & =2\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +2h\beta \frac {kD}{dh^{2}} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}=\frac {kD}{dh^{2}}u_{j-1}^{n}+u_{j}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{j+1}^{n}+\frac {k}{d}g_{j}^{n}\]
Let \(r_{1}=\left ( 1-2\frac {kD}{dh^{2}}-\frac {ak}{d}\right ) ,r_{2}=\frac {kD}{dh^{2}},\),\(r_{3}=\) \(\frac {k}{d}\), Hence the system becomes (now nodes 1 and \(N\) are unknowns and added)
\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} =\begin {bmatrix} r_{1} & 2r_{2} & 0 & 0 & 0 & 0 & 0\\ r_{2} & r_{1} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r2 & r_{1} & r_{2} & 0 & 0 & 0\\ 0 & 0 & r_{2} & \ddots & r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & 0 & r_{2} & r_{1} & r_{2}\\ 0 & 0 & 0 & 0 & 0 & 2r_{2} & r_{1}\end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} 2h\alpha ^{n}\ r_{2}\\ r_{3}g_{2}^{n}\\ r_{3}g_{3}^{n}\\ r_{3}g_{4}^{n}\\ \vdots \\ r_{3}g_{N-1}^{n}\\ 2h\beta ^{n}\ r_{2}\end {bmatrix} \]
Derivation of forward Euler for Neumann on left and Dirichlet on right Using this numbering
Assume that \(u_{N}=\beta \) at node \(N\) and \(u_{t}=\alpha \) at node \(1\) (these are the Dirichlet and Neumann boundary conditions).
The unknowns in this case are nodes \(1\cdots N-1\), since \(u_{N}\) is known from Dirichlet boundary conditions. Add a ghost node \(0\) to the left of node \(1\), and approximating \(\alpha \) gives
\[ \alpha =\frac {u_{0}-u_{2}}{2h}\]
hence
\begin {equation} u_{0}=2h\alpha +u_{2} \tag {3} \end {equation}
But the PDE for node 1 is
\begin {equation} u_{1}^{n+1}=\frac {kD}{dh^{2}}u_{0}^{n}+u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n} \tag {4} \end {equation}
Substitute (3) into (4) gives
\begin {align*} u_{1}^{n+1} & =\frac {kD}{dh^{2}}\left ( 2h\alpha +u_{2}\right ) +u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{2}^{n}\\ & =u_{1}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\left ( 2\frac {kD}{dh^{2}}\right ) u_{2}^{n}+2\ h\ \alpha \frac {kD}{dh^{2}} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}=\frac {kD}{dh^{2}}u_{j-1}^{n}+u_{j}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{j+1}^{n}+\frac {k}{d}g_{j}^{n}\]
Let \(r_{1}=\left ( 1-2\frac {kD}{dh^{2}}-\frac {ak}{d}\right ) ,r_{2}=\frac {kD}{dh^{2}},r_{3}=\) \(\frac {k}{d}\), Hence the system becomes
\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} =\begin {bmatrix} r_{1} & 2r_{2} & 0 & 0 & 0 & 0 & 0\\ r_{2} & r_{1} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{1} & r_{2} & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \ddots & \vdots & 0\\ 0 & 0 & 0 & r_{2} & r_{1} & r_{2} & 0\\ 0 & 0 & 0 & 0 & r_{2} & r_{1} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} 2\ h\ \alpha ^{n}r_{2}\\ r_{3}g_{2}^{n}\\ r_{3}g_{3}^{n}\\ r_{3}g_{4}^{n}\\ \vdots \\ r_{3}g_{N-1}^{n}+r_{2}\beta ^{n}\\ \beta ^{n+1}\end {bmatrix} \]
Derivation of forward Euler for Neumann on right and Dirichlet on left Using this numbering
Assume that \(u_{0}=\alpha \) at node \(1\) and \(u_{t}=\beta \) at node \(N\) (these are the Dirichlet and Neumann boundary conditions).
The unknowns in this case are nodes \(2\cdots N\), since \(u_{1}\) is known from Dirichlet boundary conditions. Add a ghost node \(N+1\) to the right of node \(N\), and approximating \(\beta \left ( t\right ) \) gives
\[ \beta =\frac {u_{N-1}-u_{N+1}}{2h}\]
hence
\begin {equation} u_{N+1}=2h\beta +u_{N-1} \tag {3} \end {equation}
But the PDE for node N is
\begin {equation} u_{N}^{n+1}=\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{N+1}^{n} \tag {4} \end {equation}
Substitute (3) into (4) gives
\begin {align*} u_{N}^{n+1} & =\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}\left ( 2h\beta +u_{N-1}\right ) \\ & =2\frac {kD}{dh^{2}}u_{N-1}^{n}+u_{N}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +2h\beta \frac {kD}{dh^{2}} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}=\frac {kD}{dh^{2}}u_{j-1}^{n}+u_{j}^{n}\left ( 1-2\frac {k\ D}{d\ h^{2}}-\frac {ak}{d}\right ) +\frac {kD}{dh^{2}}u_{j+1}^{n}+\frac {k}{d}g_{j}^{n}\]
Let \(r_{1}=\left ( 1-2\frac {kD}{dh^{2}}-\frac {ak}{d}\right ) ,r_{2}=\frac {kD}{dh^{2}},\),\(r_{3}=\) \(\frac {k}{d}\), Hence the system becomes
\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} =\begin {bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{1} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{1} & r_{2} & 0 & 0 & 0\\ 0 & 0 & r_{2} & r_{1} & r_{2} & 0 & 0\\ 0 & \vdots & \vdots & \vdots & \ddots & \ddots & \vdots \\ 0 & 0 & 0 & 0 & r_{2} & r_{1} & r_{2}\\ 0 & 0 & 0 & 0 & 0 & 2r_{2} & r_{1}\end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} \alpha ^{n+1}\\ r_{2}\alpha ^{n}+r_{3}g_{2}^{n}\\ r_{3}g_{2}^{n}\\ r_{3}g_{3}^{n}\\ r_{3}g_{4}^{n}\\ \vdots \\ r_{3}g_{N-1}^{n}\\ 2h\beta ^{n}\ r_{2}\end {bmatrix} \]
Derivation of C-N for periodic boundary conditions From part(a) above,
\[ u_{j}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {a\Delta t}{2d}\right ) -\frac {kD}{2dh^{2}}u_{j-1}^{n+1}-\frac {kD}{2dh^{2}}u_{j+1}^{n+1}=u_{j}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{j-1}^{n}\frac {kD}{2dh^{2}}+u_{j+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \]
Periodic boundary conditions implies \(u\left ( 0,t\right ) =u\left ( 1,t\right ) \), Hence there is an extra one unknown (in addition to the internal nodes). Either \(u\left ( 0,t\right ) \) or \(u\left ( 1,t\right ) \) can be selected since they have the same value. When selecting the right end node, then \(u_{N+1}\) becomes an unknown to be added to the internal nodes. Using the following diagram
Then for node 2 \begin {align*} u_{2}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {a\Delta t}{2d}\right ) -\frac {kD}{2dh^{2}}u_{1}^{n+1}-\frac {kD}{2dh^{2}}u_{3}^{n+1} & =u_{2}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{1}^{n}\frac {kD}{2dh^{2}}+u_{3}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{2}^{n}+g_{2}^{n+1}\right ) \\ u_{2}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {a\Delta t}{2d}\right ) -\frac {kD}{2dh^{2}}u_{N}^{n+1}-\frac {kD}{2dh^{2}}u_{3}^{n+1} & =u_{2}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{N}^{n}\frac {kD}{2dh^{2}}+u_{3}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{2}^{n}+g_{2}^{n+1}\right ) \end {align*}
And for node N
\begin {align*} u_{N}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {a\Delta t}{2d}\right ) -\frac {kD}{2dh^{2}}u_{N-1}^{n+1}-\frac {kD}{2dh^{2}}u_{N+1}^{n+1} & =u_{N}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{N-1}^{n}\frac {kD}{2dh^{2}}+u_{N+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{N}^{n}+g_{N}^{n+1}\right ) \\ u_{N}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {a\Delta t}{2d}\right ) -\frac {kD}{2dh^{2}}u_{N-1}^{n+1}-\frac {kD}{2dh^{2}}u_{2}^{n+1} & =u_{N}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{N-1}^{n}\frac {kD}{2dh^{2}}+u_{2}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{N}^{n}+g_{N}^{n+1}\right ) \end {align*}
Let \begin {align*} r_{1} & =\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) \\ r_{2} & =\frac {kD}{2dh^{2}}\\ r_{3} & =\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) \\ r_{4} & =\frac {k}{2d}\\ g^{n}+g^{n+1} & \equiv \tilde {g} \end {align*}
Converting to matrix form gives\begin {align} \overset {\mathbf {A}}{\overbrace {\begin {bmatrix} r_{1} & -r_{2} & 0 & 0 & -r_{2}\\ -r_{2} & r_{1} & -r_{2} & 0 & 0\\ 0 & -r_{2} & r_{1} & 0 & 0\\ 0 & 0 & 0 & \ddots & \vdots \\ -r_{2} & 0 & 0 & -r_{2} & r_{1}\end {bmatrix} }}\overset {\mathbf {x}}{\overbrace {\begin {bmatrix} u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} }} & =\overset {\mathbf {b}}{\overbrace {\begin {bmatrix} r_{3} & r_{2} & 0 & 0 & r_{2}\\ r_{2} & r_{3} & r_{2} & 0 & 0\\ 0 & r_{2} & r_{3} & r_{2} & 0\\ 0 & 0 & 0 & \ddots & \vdots \\ r_{2} & 0 & 0 & r_{2} & r_{3}\end {bmatrix}\begin {bmatrix} u_{2}^{n}\\ u_{3}^{n}\\ u_{4}^{n}\\ \vdots \\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} r_{4}\tilde {g}_{2}\\ r_{4}\tilde {g}_{3}\\ r_{4}\tilde {g}_{4}\\ \vdots \\ r_{4}\tilde {g}_{N}\end {bmatrix} }}\tag {3}\\ u_{1}^{n+1} & =u_{N}^{n+1}\nonumber \end {align}
Derivation of C-N for Neumann boundary conditions both ends Using this numbering
Assuming that \(u_{t}=\alpha \) at node \(1\) and \(u_{t}=\beta \) at node \(N\) (these are the Neumann boundary conditions). Add a ghost node \(0\) to the left of node \(1\), and approximating \(\alpha \) gives
\[ \alpha =\frac {u_{0}-u_{2}}{2h}\]
hence
\begin {equation} u_{0}=2h\alpha +u_{2} \tag {1} \end {equation}
But the PDE for node 1 is
\begin {equation} r_{1}u_{1}^{n+1}-r_{2}u_{0}^{n+1}-r_{2}u_{2}^{n+1}=r_{3}u_{1}^{n}+r_{2}u_{0}^{n}+r_{2}u_{2}^{n} \tag {2} \end {equation}
substitute (1) into (2)
\begin {align*} r_{1}u_{1}^{n+1}-r_{2}\left ( 2h\alpha ^{n+1}+u_{2}^{n+1}\right ) -r_{2}u_{2}^{n+1} & =r_{3}u_{1}^{n}+r_{2}\left ( 2h\alpha ^{n}+u_{2}^{n}\right ) +r_{2}u_{2}^{n}\\ r_{1}u_{1}^{n+1}-2r_{2}u_{2}^{n+1} & =r_{3}u_{1}^{n}+2r_{2}u_{2}^{n}+2r_{2}h\alpha ^{n}+2r_{2}h\alpha ^{n+1} \end {align*}
Similarly for the right end. Add a ghost node \(N+1\) to the right of node \(N\), and approximating \(\beta \left ( t\right ) \) gives
\[ \beta =\frac {u_{N-1}-u_{N+1}}{2h}\]
hence
\begin {equation} u_{N+1}=2h\beta +u_{N-1} \tag {3} \end {equation}
But the PDE for node N is
\begin {equation} r_{1}u_{N}^{n+1}-r_{2}u_{N-1}^{n+1}-r_{2}u_{N+1}^{n+1}=r_{3}u_{N}^{n}+r_{2}u_{N-1}^{n}+r_{2}u_{N+1}^{n} \tag {4} \end {equation}
substitute (3) into (4)
\begin {align*} r_{1}u_{N}^{n+1}-r_{2}u_{N-1}^{n+1}-r_{2}\left ( 2h\beta ^{n+1}+u_{N-1}^{n+1}\right ) & =r_{3}u_{N}^{n}+r_{2}u_{N-1}^{n}+r_{2}\left ( 2h\beta +u_{N-1}^{n}\right ) \\ r_{1}u_{N}^{n+1}-2r_{2}u_{N-1}^{n+1} & =r_{3}u_{N}^{n}+2r_{2}u_{N-1}^{n}+2r_{2}h\beta +2r_{2}h\beta ^{n+1} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) -\frac {kD}{2dh^{2}}u_{j-1}^{n+1}-\frac {kD}{2dh^{2}}u_{j+1}^{n+1}=u_{j}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{j-1}^{n}\frac {kD}{2dh^{2}}+u_{j+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \] Let \begin {align*} r_{1} & =\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) \\ r_{2} & =\frac {kD}{2dh^{2}}\\ r_{3} & =\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) \\ r_{4} & =\frac {k}{2d} \end {align*}
gives
\[ r_{1}u_{j}^{n+1}-r_{2}u_{j-1}^{n+1}-r_{2}u_{j+1}^{n+1}=r_{3}u_{j}^{n}+r_{2}u_{j-1}^{n}+r_{2}u_{j+1}^{n}+r_{4}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \]
Then the above becomes
\[ \overset {\mathbf {A}}{\overbrace {\begin {bmatrix} r_{1} & -2r_{2} & 0 & 0 & 0 & 0\\ -r_{2} & r_{1} & -r_{2} & 0 & 0 & 0\\ 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & -r_{2} & r_{1} & -r_{2}\\ 0 & 0 & 0 & 0 & -2r_{2} & r_{1}\end {bmatrix} }}\overset {\mathbf {x}}{\overbrace {\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} }}=\overset {\mathbf {b}}{\overbrace {\begin {bmatrix} r_{3} & 2r_{2} & 0 & 0 & 0 & 0\\ r_{2} & r_{3} & r_{2} & 0 & 0 & 0\\ 0 & r_{2} & r_{3} & r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & r_{2} & r_{3} & r_{2}\\ 0 & 0 & 0 & 0 & 2r_{2} & r_{3}\end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} 2r_{2}h\alpha ^{n}+2r_{2}h\alpha ^{n+1}\\ r_{4}\left ( g_{2}^{n}+g_{2}^{n+1}\right ) \\ r_{4}\left ( g_{3}^{n}+g_{3}^{n+1}\right ) \\ \vdots \\ r_{4}\left ( g_{N-1}^{n}+g_{N-1}^{n+1}\right ) \\ 2r_{2}h\beta +2r_{2}h\beta ^{n+1}\end {bmatrix} }}\]
Derivation of C-N for Neumann on left and Dirichlet on right Using this numbering
Assume that \(u_{t}=\alpha \) at node \(1\) and \(u_{N}=\beta \) at node \(N\) (these are the Neumann and Dirichlet boundary conditions). Add a ghost node \(0\) to the left of node \(1\), and approximating \(\alpha \) gives
\[ \alpha =\frac {u_{0}-u_{2}}{2h}\]
hence
\begin {equation} u_{0}=2h\alpha +u_{2} \tag {1} \end {equation}
But the PDE for node 1 is
\begin {equation} r_{1}u_{1}^{n+1}-r_{2}u_{0}^{n+1}-r_{2}u_{2}^{n+1}=r_{3}u_{1}^{n}+r_{2}u_{0}^{n}+r_{2}u_{2}^{n} \tag {2} \end {equation}
substitute (1) into (2)
\begin {align*} r_{1}u_{1}^{n+1}-r_{2}\left ( 2h\alpha ^{n+1}+u_{2}^{n+1}\right ) -r_{2}u_{2}^{n+1} & =r_{3}u_{1}^{n}+r_{2}\left ( 2h\alpha ^{n}+u_{2}^{n}\right ) +r_{2}u_{2}^{n}\\ r_{1}u_{1}^{n+1}-2r_{2}u_{2}^{n+1} & =r_{3}u_{1}^{n}+2r_{2}u_{2}^{n}+2r_{2}h\alpha ^{n}+2r_{2}h\alpha ^{n+1} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) -\frac {kD}{2dh^{2}}u_{j-1}^{n+1}-\frac {kD}{2dh^{2}}u_{j+1}^{n+1}=u_{j}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{j-1}^{n}\frac {kD}{2dh^{2}}+u_{j+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \] Let \begin {align*} r_{1} & =\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) \\ r_{2} & =\frac {kD}{2dh^{2}}\\ r_{3} & =\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) \\ r_{4} & =\frac {k}{2d} \end {align*}
gives
\[ r_{1}u_{j}^{n+1}-r_{2}u_{j-1}^{n+1}-r_{2}u_{j+1}^{n+1}=r_{3}u_{j}^{n}+r_{2}u_{j-1}^{n}+r_{2}u_{j+1}^{n}+r_{4}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \]
Then the system becomes
\[ \overset {\mathbf {A}}{\overbrace {\begin {bmatrix} r_{1} & -2r_{2} & 0 & 0 & 0 & 0 & 0\\ -r_{2} & r_{1} & -r_{2} & 0 & 0 & 0 & 0\\ 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & 0\\ 0 & 0 & 0 & -r_{2} & r_{1} & -r_{2} & 0\\ 0 & 0 & 0 & 0 & -r_{2} & r_{1} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end {bmatrix} }}\overset {\mathbf {x}}{\overbrace {\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-2}^{n+1}\\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} }}=\overset {\mathbf {b}}{\overbrace {\begin {bmatrix} r_{3} & 2r_{2} & 0 & 0 & 0 & 0 & 0\\ r_{2} & r_{3} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{3} & r_{2} & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & 0\\ 0 & 0 & 0 & r_{2} & r_{3} & r_{2} & 0\\ 0 & 0 & 0 & 0 & 2r_{2} & r_{3} & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N-2}^{n}\\ u_{N-1}^{n}\\ u_{N}^{n}\end {bmatrix} +\begin {bmatrix} 2r_{2}h\alpha ^{n}+2r_{2}h\alpha ^{n+1}\\ r_{4}\left ( g_{2}^{n}+g_{2}^{n+1}\right ) \\ r_{4}\left ( g_{3}^{n}+g_{3}^{n+1}\right ) \\ \vdots \\ r_{4}\left ( g_{N-2}^{n}+g_{N-2}^{n+1}\right ) \\ r_{4}\left ( g_{N-1}^{n}+g_{N-1}^{n+1}\right ) +r_{2}\beta ^{n}\\ \beta ^{n+1}\end {bmatrix} }}\]
Derivation of C-N for Neumann on right and Dirichlet on left Using this numbering
Assume that \(u_{1}=\alpha \) at node \(1\) and \(u_{t}=\beta \) at node \(N\) (these are the Dirichlet and Neumann boundary conditions). Add a ghost node \(N+1\) to the right of node \(N\), and approximating \(\beta \) gives
\[ \beta =\frac {u_{N-1}-u_{N+1}}{2h}\]
hence
\begin {equation} u_{N+1}=2h\beta +u_{N-1} \tag {1} \end {equation}
But the PDE for node N is
\begin {equation} r_{1}u_{N}^{n+1}-r_{2}u_{N-1}^{n+1}-r_{2}u_{N+1}^{n+1}=r_{3}u_{N}^{n}+r_{2}u_{N-1}^{n}+r_{2}u_{N+1}^{n} \tag {2} \end {equation}
substitute (1) into (2)
\begin {align*} r_{1}u_{N}^{n+1}-r_{2}u_{N-1}^{n+1}-r_{2}\left ( 2h\beta ^{n+1}+u_{N-1}^{n+1}\right ) & =r_{3}u_{N}^{n}+r_{2}u_{N-1}^{n}+r_{2}\left ( 2h\beta ^{n}+u_{N-1}^{n}\right ) \\ r_{1}u_{N}^{n+1}-2r_{2}u_{N-1}^{n+1} & =r_{3}u_{N}^{n}+2r_{2}u_{N-1}^{n}+2r_{2}h\beta ^{n}+2r_{2}h\beta ^{n+1} \end {align*}
Nodes \(j=2\cdots N-1\) remain the same as before. In other words
\[ u_{j}^{n+1}\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) -\frac {kD}{2dh^{2}}u_{j-1}^{n+1}-\frac {kD}{2dh^{2}}u_{j+1}^{n+1}=u_{j}^{n}\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) +u_{j-1}^{n}\frac {kD}{2dh^{2}}+u_{j+1}^{n}\frac {kD}{2dh^{2}}+\frac {k}{2d}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \] Let \begin {align*} r_{1} & =\left ( 1+\frac {kD}{dh^{2}}+\frac {ak}{2d}\right ) \\ r_{2} & =\frac {kD}{2dh^{2}}\\ r_{3} & =\left ( 1-\frac {kD}{dh^{2}}-\frac {ak}{2d}\right ) \\ r_{4} & =\frac {k}{2d} \end {align*}
gives
\[ r_{1}u_{j}^{n+1}-r_{2}u_{j-1}^{n+1}-r_{2}u_{j+1}^{n+1}=r_{3}u_{j}^{n}+r_{2}u_{j-1}^{n}+r_{2}u_{j+1}^{n}+r_{4}\left ( g_{j}^{n}+g_{j}^{n+1}\right ) \]
Then the system becomes
\[ \overset {\mathbf {A}}{\overbrace {\begin {bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{1} & -r_{2} & 0 & 0 & 0 & 0\\ 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0 & 0\\ 0 & 0 & -r_{2} & r_{1} & -r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & 0\\ 0 & 0 & 0 & 0 & -r_{2} & r_{1} & -r_{2}\\ 0 & 0 & 0 & 0 & 0 & -2r_{2} & r_{1}\end {bmatrix} }}\overset {\mathbf {x}}{\overbrace {\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ u_{4}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} }}=\overset {\mathbf {b}}{\overbrace {\begin {bmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & r_{3} & r_{2} & 0 & 0 & 0 & 0\\ 0 & r_{2} & r_{3} & r_{2} & 0 & 0 & 0\\ 0 & 0 & r_{2} & r_{3} & r_{2} & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & 0\\ 0 & 0 & 0 & 0 & r_{2} & r_{3} & r_{2}\\ 0 & 0 & 0 & 0 & 0 & 2r_{2} & r_{3}\end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ \vdots \\ u_{N-2}^{n}\\ u_{N-1}^{n}\end {bmatrix} +\begin {bmatrix} \alpha ^{n+1}\\ r_{4}\left ( g_{2}^{n}+g_{2}^{n+1}\right ) +r_{2}\alpha ^{n}\\ r_{4}\left ( g_{3}^{n}+g_{3}^{n+1}\right ) \\ r_{4}\left ( g_{4}^{n}+g_{4}^{n+1}\right ) \\ \vdots \\ r_{4}\left ( g_{N-2}^{n}+g_{N-2}^{n+1}\right ) \\ 2r_{2}h\beta ^{n}+2r_{2}h\beta ^{n+1}\end {bmatrix} }}\]
The C-N scheme was programmed in Matlab and then run on the above problem. The following shows the result
In the above plot, the red line represents initial conditions (the step function shifted to the right by 0.5) and the blue line represents the final numerical solution at time \(t=1\) second. The following plot is a closer look at the grid near \(x=\frac {1}{2}\) showing the initial conditions
It is clear the numerical solution is not accurate as it does not match what is expected to occur physically which is for initial data to diffuse. The initial data contained high spatial frequency that should have been smoothed out rapidly. The final numerical solution is not smooth and contain high spatial frequency components which should have been attenuated by the time the run is completed.
The exact solution in the Fourier space is \(\hat {u}\left ( \xi ,t\right ) =\hat {u}\left ( \xi ,0\right ) e^{-D\xi ^{2}t}\), where \(\hat {u}\left ( \xi ,0\right ) \) is the spectrum or Fourier coefficients of the initial condition \(u\left ( x,0\right ) \). This shows that modes with large spatial frequency (large wave number \(\xi )\) will attenuate the fastest due to the negative exponential decay effect. But this was not observed in the above numerical solution.
C-N is a stable scheme (A-stable), but can be inaccurate if the time step used is large relative to the space step or if initial conditions contain large spatial frequency components. In C-N, the time step needs to be about the same order of value as the space step for the scheme to give accurate numerical results. (This is because in C-N, the order of accuracy of space and time are the same, as was found in problem 1).
Therefore, it appears that C-N scheme does not handle discontinuities in initial conditions well as this result shows.
In the next part, the amplification factor for C-N is determined, and a mathematical explanation for the above result is given.
The C-N scheme for \(u_{t}=Du_{xx}\) is given by
\begin {equation} -ru_{j-1}^{n+1}+u_{j}^{n+1}\left ( 1+2r\right ) -ru_{j+1}^{n+1}=ru_{j-1}^{n}+u_{j}^{n}\left ( 1-2r\right ) +ru_{j+1}^{n} \tag {1} \end {equation}
Where \(r=\frac {\Delta tD}{2h^{2}}\). Von Neumann analysis is used to determine the magnification factor2. Assume \(u_{j}^{n}=e^{i\xi x_{j}}\), and \(u_{j}^{n+1}=g\left ( \xi \right ) e^{i\xi x_{j}}\) then (1) becomes3\begin {align*} -r\left ( g\ e^{i\xi x_{j}}e^{-i\xi h}\right ) +g\ e^{i\xi x_{j}}\left ( 1+2r\right ) -r\left ( g\ e^{i\xi x_{j}}e^{i\xi h}\right ) & =r\left ( e^{i\xi x_{j}}e^{-i\xi h}\right ) +e^{i\xi x_{j}}\left ( 1-2r\right ) +re^{i\xi x_{j}}e^{i\xi h}\\ g\left ( -2r\cos \left ( \xi h\right ) +1+2r\right ) & =2r\cos \left ( \xi h\right ) +1-2r\\ g\left ( \xi \right ) & =\frac {1+2r\cos \left ( \xi h\right ) -2r}{1-2r\cos \left ( \xi h\right ) +2r} \end {align*}
Hence, the magnification factor is \[ g\left ( \xi \right ) =\frac {1+\frac {D\Delta t}{h^{2}}\left ( \cos \left ( \xi h\right ) -1\right ) }{1-\frac {D\Delta t}{h^{2}}\left ( \cos \left ( \xi h\right ) -1\right ) }\] Let \(\xi \) be written as \(\xi _{p}\) in the above in order to examine \(g\) in terms of specific wave number \(\xi _{p}\) (which has units of radians per unit length). The above becomes\begin {equation} g\left ( \xi _{p}\right ) =\frac {1+\frac {D\Delta t}{h^{2}}\left ( \cos \left ( \xi _{p}h\right ) -1\right ) }{1-\frac {D\Delta t}{h^{2}}\left ( \cos \left ( \xi _{p}h\right ) -1\right ) } \tag {2} \end {equation} But \(\xi _{p}=p\pi \) and \(p=1\cdots N\), with \(h=\frac {1}{N+1}\)where the line was discretized using the standard grid convention
Therefore \(\cos \left ( \xi _{p}h\right ) =\cos \left ( p\pi h\right ) =\cos \left ( \frac {p\pi }{N+1}\right ) .\) The largest frequency occurs when \(p=N\), because then \(h\) is smallest, and the smallest frequency occurs when \(p=1\). Hence, there are a total of \(N\) Fourier modes when representing initial data as Fourier series. Now the magnification factor in (2) becomes
\begin {equation} g\left ( \xi _{p}\right ) =\frac {1+\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }{1-\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) } \tag {3} \end {equation}
But since \(\left \vert \cos \left ( p\pi h\right ) \right \vert \leq 1\), then \(g\left ( \xi _{p}\right ) \) is less than \(1\) in magnitude for any \(p\), implying that C-N is stable. To determine the magnitude of \(g\left ( \xi _{p}\right ) \) when the mode has the largest frequency, let \(p\pi h=\frac {N}{N+1}\pi \approx \pi \) in (3), resulting in
\begin {equation} g\left ( \xi _N\right ) =\frac {1-2\frac {D\Delta t}{h^2}}{1+2\frac {D\Delta t}{h^2}} \tag {4} \end {equation}
When the time step \(\Delta t\) \(\gg h\), then \(\frac {D\Delta t}{h^{2}}\gg 1,\) and in the limit \(\left \vert g\left ( \xi _{N}\right ) \right \vert \rightarrow 1\). This shows that large frequency modes will decay very slowly because \(g\left ( \xi _{p}\right ) \) is now close to \(1\). No attenuation will occur between each application of the update matrix or between each time step.
The above explains the result seen in part (a). Large frequency components did not decay fast as was expected, because the time step used was much larger than the space step. The problem asked us to use \(\Delta t\) \(=0.1\) and \(h=0.01\), which gives \(\frac {D\Delta t}{h^{2}}=\frac {0.1}{0.01^{2}}=1000\) and hence \(g\left ( \xi _{N}\right ) \rightarrow \left \vert \frac {1-2000}{1+2000}\right \vert \rightarrow 0.9999\), and since this is almost one, then large frequency modes did not attenuate with each time step. The amplification factor needs to be small for attenuation to occur fast.
The following is a plot of \(g\left ( \xi _{p}\right ) \) showing how the amplification factor changes as function of \(p\) for the case of \(\Delta t=0.1\) and \(h=0.01\), and \(D=1\). It shows what was found above, that at large frequency where \(p\) is close to \(N\) will have a correspondingly large
To determine what value of \(\frac {D\Delta t}{h^{2}}\) is required to make the large frequency mode decay right away, let \(\frac {D\Delta t}{h^{2}}=\frac {1}{2}\) in (4), this gives \(g\left ( \xi _{N}\right ) =0\), which implies that large frequency mode will be knocked out right away. Here is a plot of \(g\left ( \xi _{N}\right ) =\frac {1-2\frac {D\Delta t}{h^{2}}}{1+2\frac {D\Delta t}{h^{2}}}\) as a function of \(\frac {D\Delta t}{h^{2}}\) showing that when \(\frac {D\Delta t}{h^{2}}=0\) then the magnification factor is minimum. This is only for mode \(p=N\).
Conclusion If initial data contained large difference in value over very short distances (in other words, large spatial frequencies) such as given in this problem, producing discontinuity in data and its space derivative, and when the time step is large compared to the space step, then the numerical solution produced by C-N will not be accurate since large frequency modes will not attenuate.
To compensate for large frequency present in initial data, the ratio\(\frac {D\Delta t}{h^{2}}\)needs to be made close to \(0.5\) as possible. It might be better not to use C-N at all in such case and look for a scheme which does not have this problem.
notice that condition that \(\frac {D\Delta t}{h^{2}}=\frac {1}{2}\) found above, is the same value for the upper limit for the absolute stability condition for the forward Euler discretization scheme for the 1D diffusion problem.
The time step was reduced and the program was run for \(1\) second. When the time step was reduced all the way to \(\Delta t=0.01\) then the final solution appeared smooth every where and in particular at \(x=0.5\). The following diagram illustrates this.
When using time step of \(0.01\,\ \sec ,\) \(100\) steps are used. From the last part it was found that \(g\left ( \xi _{N}\right ) =\frac {1-2\frac {D\Delta t}{h^{2}}}{1+2\frac {D\Delta t}{h^{2}}}\), hence \(\left \vert g\right \vert =\left \vert \frac {1-2\frac {0.01}{0.02^{2}}}{1+2\frac {0.01}{0.02^{2}}}\right \vert =\allowbreak 0.960\,78\), and therefore \(\allowbreak \left \vert g\right \vert ^{100}\rightarrow 0.960\,78^{100}\approx 0\) showing that by the 100\(^{th}\) iteration, the large mode frequency have completely smoothed out as verified by the above plots.
In the plot below, the magnification factor \(\left \vert g\left ( \xi \right ) \right \vert \) is shown for \(t=0.01\) and \(h=0.02\) showing that at \(p=50\) , \(g\left ( \xi \right ) =0.960\,78\). Compare this value with the one used in part(b) which was \(0.999.\)
To obtain smooth solution immediately after one time step, the required time step to accomplish this, can be determined from the condition for optimal amplification factor found in the last part which is given by \(\frac {D\Delta t}{h^{2}}=\frac {1}{2}.\) From this relation, and when \(h=0.02\) and \(D=1\), the time step will be \(\Delta t=0.0002\). This value of time step produces a smooth solution immediately (one step). This was confirmed, and here is the numerical solution, after only one time step, using \(\Delta t=0.0002,h=0.02\) and \(D=1\)
It is also possible to determine which \(\Delta t\) achieves a given specific attenuation of the high mode. Suppose it is required to attenuate the high mode to \(0.001\) of its initial amplitude at the end of \(1\) second run. Therefore, this means that\begin {align*} \left \vert g\left ( \xi _{N}\right ) \right \vert ^{\frac {1}{\Delta t}} & =0.001\\ \left \vert \frac {1-2\frac {D\Delta t}{h^{2}}}{1+2\frac {D\Delta t}{h^{2}}}\right \vert ^{\frac {1}{\Delta t}} & =0.001 \end {align*}
Taking logs, and using \(h=0.02\) and \(D=1\) results in \begin {align*} \frac {1}{\Delta t}\log \left ( \frac {1-5000\Delta t}{1+5000\Delta t}\right ) & =-3\\ \log \left ( 1-5000\Delta t\right ) -\log \left ( 1+5000\Delta t\right ) & =-3\Delta t \end {align*}
The above is not a linear equation, but can be numerically solved for the root \(\Delta t\). For the above example, \(\Delta t\) came out to be \(0.00761\) seconds. This means that when using \(\Delta t\) \(=0.00761\) sec, \(h=0.02\), and \(D=1\,,\) then the largest frequency harmonic will have its amplitude attenuated to \(0.01\%\) of its original value after \(1\) second run.
Making the time step smaller and smaller, while keeping the ratio \(\frac {\Delta t}{h}\) fixed, produces the following result (all runs are for one second). In this example, the ratio was kept at \(5\). The following sequence of ratios are used \(\left \{ \frac {0.1}{0.02},\frac {0.01}{0.002},\frac {0.001}{0.0002},\frac {0.0001}{0.00002}\right \} \) to generate the following solution after one second run
Now, recall from part (b) that the magnification factor for the largest Fourier mode was given by \[ g\left ( \xi _{N}\right ) =\frac {1-2\frac {D\Delta t}{h^{2}}}{1+2\frac {D\Delta t}{h^{2}}}\] When \(\frac {\Delta t}{h}\) is held constant, say \(C\), then the above becomes\[ g\left ( \xi _{N}\right ) =\frac {1-Ch^{-1}}{1+Ch^{-1}}\] In the limit, as \(h\) \(\rightarrow 0\) then \(g\left ( \xi _{N}\right ) \rightarrow 1\), which implies, as was found in part(b), that large Fourier modes (high \(p\) values) will not be attenuated. This is confirmed by the plots above.
Using backward Euler in place of C-N When using backward Euler. The finite difference scheme for \(u_{t}=Du_{xx}\) becomes \begin {align} \frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t} & =f\left ( u^{n+1}\right ) \tag {1}\\ & =D\frac {u_{j-1}^{n+1}-2u_{j}^{n+1}+u_{j+1}^{n+1}}{h^{2}}\nonumber \end {align}
Applying Von Neumann analysis, let \(u_{j}^{n}=e^{i\xi x_{j}}\), and \(u_{j}^{n+1}=g\left ( \xi \right ) e^{i\xi x_{j}}\), then (1) becomes\begin {align*} g\ e^{i\xi x_{j}}e^{i\xi h}-\ e^{i\xi x_{j}} & =\frac {D\Delta t}{h^{2}}g\left ( e^{i\xi x_{j}}e^{-i\xi h}-2e^{i\xi x_{j}}+e^{i\xi x_{j}}e^{i\xi h}\right ) \\ g\ e^{i\xi h}-1 & =\frac {D\Delta t}{h^{2}}g\left ( e^{-i\xi h}-2+e^{i\xi h}\right ) \\ g\left ( e^{i\xi h}-\frac {D\Delta t}{h^{2}}\left ( 2\cos \left ( \xi h\right ) -2\right ) \right ) & =1\\ g\left ( \xi \right ) & =\frac {1}{\left ( e^{i\xi h}-\frac {D\Delta t}{h^{2}}\left ( 2\cos \left ( \xi h\right ) -2\right ) \right ) } \end {align*}
Therefore\[ g\left ( \xi _{p}\right ) =\frac {1}{\left ( e^{i\xi _{p}h}-\frac {D\Delta t}{h^{2}}\left ( 2\cos \left ( \xi _{p}h\right ) -2\right ) \right ) }\] But \(\xi _{p}=p\pi \) and \(p=1\cdots N\), and to evaluate what happens to \(g\left ( \xi _{p}\right ) \) at the largest spatial frequencies, let \(\xi _{p}=N\pi \) and the above becomes\[ g\left ( \xi _{N}\right ) =\frac {1}{\left ( e^{iN\pi h}-\frac {D\Delta t}{h^{2}}\left ( 2\cos \left ( \frac {N}{N+1}\pi \right ) -2\right ) \right ) }\] But \(\frac {N}{N+1}\approx 1\) and \(e^{iN\pi h}\approx e^{i\pi }=\cos \left ( \pi \right ) +i\sin \pi =-1\), then the above becomes\[ g\left ( \xi _{N}\right ) =\frac {1}{\left ( -1-\frac {D\Delta t}{h^{2}}\left ( 2\cos \pi -2\right ) \right ) }\] Therefore\[ g\left ( \xi _N\right ) =\frac {1}{4\frac {D\Delta t}{h^2}-1} \] When \(\Delta t\gg h\) then \(\left \vert g\left ( \xi _{N}\right ) \right \vert \rightarrow \frac {1}{1-1}\rightarrow 0\) which implies that large frequency modes will be knocked out fast. This is opposite to the situation observed using C-N. Hence backward Euler does not have the same problem with large spatial frequencies in initial data. But notice that when \(\frac {D\Delta t}{h^{2}}=\frac {1}{2}\), now \(g\left ( \xi _{N}\right ) =\frac {1}{2-1}=1\), which means that large frequency mode will not decay or decay very slowly, This is also opposite of what was found for C-N.
Derivation of part (b) by direct application of DFT (the harder way) The C-N scheme for \(u_{t}=Du_{xx}\) given by\begin {equation} -ru_{j-1}^{n+1}+u_{j}^{n+1}\left ( 1+2r\right ) -ru_{j+1}^{n+1}=ru_{j-1}^{n}+u_{j}^{n}\left ( 1-2r\right ) +ru_{j+1}^{n} \tag {2} \end {equation} Assuming the problem is on the whole real line, then \begin {equation} u_{j}^{n}=\frac {1}{\sqrt {2\pi }}{\displaystyle \int \limits _{-\frac {\pi }{h}}^{\frac {\pi }{h}}} \hat {u}^{n}\left ( \xi \right ) e^{i\xi x_{j}}d\xi \tag {1A} \end {equation} and\begin {equation} u_{j}^{n+1}=\frac {1}{\sqrt {2\pi }}{\displaystyle \int \limits _{-\frac {\pi }{h}}^{\frac {\pi }{h}}} \hat {u}^{n+1}\left ( \xi \right ) e^{i\xi x_{j}}d\xi \tag {1B} \end {equation} Where \(\hat {u}^{n}\left ( \xi \right ) \) is the discrete Fourier transform (DFT) of \(u_{j}^{n}.\) In what follows, \(\hat {u}^{n}\) is written instead of \(\hat {u}^{n}\left ( \xi \right ) \) to make it easier to read the equations. The C-N finite difference scheme for \(u_{t}=Du_{xx}\) is given by\begin {equation} -ru_{j-1}^{n+1}+u_{j}^{n+1}\left ( 1+2r\right ) -ru_{j+1}^{n+1}=ru_{j-1}^{n}+u_{j}^{n}\left ( 1-2r\right ) +ru_{j+1}^{n} \tag {2} \end {equation} Where \(r=\frac {\Delta tD}{2h^{2}}\). Substitute (1A) and (1B) into (2), but leaving \(u_{j}^{n+1}\) as is gives
Hence \begin {align*} DFT\left ( u_{j}^{n+1}\left ( 1+2r\right ) \right ) & =\left ( \hat {u}^{n}\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) +2r\cos \left ( \xi h\right ) \hat {u}^{n+1}\right ) \\ DFT\left ( u_{j}^{n+1}\right ) & =\frac {1}{\left ( 1+2r\right ) }\left ( \hat {u}^{n}\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) +2r\cos \left ( \xi h\right ) \hat {u}^{n+1}\right ) \end {align*}
This implies that\[ \hat {u}^{n+1}=\frac {1}{\left ( 1+2r\right ) }\left ( \hat {u}^{n}\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) +2r\cos \left ( \xi h\right ) \hat {u}^{n+1}\right ) \] Solving for \(\hat {u}^{n+1}\) gives\begin {align*} \hat {u}^{n+1}-\frac {2r\cos \left ( \xi h\right ) \hat {u}^{n+1}}{\left ( 1+2r\right ) } & =\frac {1}{\left ( 1+2r\right ) }\hat {u}^{n}\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) \\ \hat {u}^{n+1}\left ( \frac {\left ( 1+2r\right ) -2r\cos \left ( \xi h\right ) }{\left ( 1+2r\right ) }\right ) & =\frac {\hat {u}^{n}\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) }{\left ( 1+2r\right ) } \end {align*}
Hence\begin {align*} \hat {u}^{n+1} & =\frac {\left ( 2r\left ( \cos \left ( \xi h\right ) -1\right ) +1\right ) }{\left ( 1+2r\right ) -2r\cos \left ( \xi h\right ) }\hat {u}^{n}\\ & =\frac {1+2r\cos \left ( \xi h\right ) -2r}{1-2r\cos \left ( \xi h\right ) +2r}\hat {u}^{n} \end {align*}
Therefore \(\hat {u}^{n+1}=g\left ( \xi \right ) \hat {u}^{n}\), where \(g\left ( \xi \right ) =\frac {1+2r\cos \left ( \xi h\right ) -2r}{1-2r\cos \left ( \xi h\right ) +2r}\), and since \(r=\frac {\Delta tD}{2h^{2}}\), then
\[ g\left ( \xi \right ) =\frac {1+\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }{1-\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }\] where \(\xi _{p}=p\pi \), is the wave number.
Another derivation for the magnification factor The magnification factor \(g\left ( \xi \right ) \) found above is the same as the eigenvalue of the update matrix of the C-N scheme. From \[ \left ( I-\frac {D\Delta t}{2}L\right ) u^{n+1}=\left ( I+\frac {D\Delta t}{2}L\right ) u^{n}\] or\[ u^{n+1}=\left ( I-\frac {D\Delta t}{2}L\right ) ^{-1}\left ( I+\frac {D\Delta t}{2}L\right ) u^{n}\] Where \(L\) is the 1D Laplacian grid operator which has eigenvalues \(\lambda _{p}=\frac {2}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) \), hence, let \(\mu _{p}\) be the eigenvalue of \(B\) above, then
\begin {align*} \mu _{p} & =\frac {1+\frac {D\Delta t}{2}\lambda _{p}}{1-\frac {D\Delta t}{2}\lambda _{p}}\\ & =\frac {1+\frac {D\Delta t}{2}\frac {2}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }{1-\frac {D\Delta t}{2}\frac {2}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }\\ & =\frac {1+\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) }{1-\frac {D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) } \end {align*}
Which is what was found in part(b)
The finite difference scheme for the diffusion problem is shown the appendix of this problem.
To obtain the absolute stability restriction, let\begin {align*} y^{n+1}-y^{n} & =\frac {\Delta t}{2}\left ( f\left ( y^{n}\right ) +f\left ( y^{\ast }\right ) \right ) \\ y^{n+1} & =y^{n}+\frac {\Delta t}{2}\left ( \lambda y^{n}+\lambda y^{\ast }\right ) \\ & =y^{n}+\frac {\Delta t}{2}\left ( \lambda y^{n}+\lambda \left ( y^{n}+\Delta t\lambda y^{n}\right ) \right ) \\ & =y^{n}+\frac {\Delta t}{2}\lambda y^{n}+\frac {\Delta t}{2}\lambda \left ( y^{n}+\Delta t\lambda y^{n}\right ) \\ & =y^{n}+\frac {\Delta t}{2}\lambda y^{n}+\frac {\Delta t}{2}\lambda y^{n}+\frac {\left ( \lambda \Delta t\right ) ^{2}}{2}y^{n}\\ & =\left ( 1+\Delta t\lambda +\frac {\left ( \lambda \Delta t\right ) ^{2}}{2}\right ) y^{n} \end {align*}
Assuming \(\Delta t\lambda =z\)\[ y^n+1=\left ( 1+z+\frac {1}{2}z^2\right ) y^n \] Hence \(R\left ( z\right ) =1+z+\frac {1}{2}z^{2}\) and for absolute stability it is required that \(\left \vert R\left ( z\right ) \right \vert \leq 1\) which leads to\begin {align*} -1 & \leq 1+z+\frac {1}{2}z^{2}\leq 1\\ -2 & \leq z+\frac {1}{2}z^{2}\leq 0 \end {align*}
A plot of \(z+\frac {1}{2}z^{2}\,\ \) shows that \(-2\leq z\leq 0\)
The above gives the interval of absolute stability for the eigenvalues. To obtain the region of absolute stability, assume \(\lambda \) can be complex in general (complex eigenvalue), which results in a disk of radius 1 centered at -1. This is the same region of stability as Forward Euler.
Notice that since \(R\left ( 0\right ) =1\) and \(\frac {dR\left ( z\right ) }{dz}=1\), then this method is also consistent and first order accurate in time.
To answer the question about any advantage of this method over forward Euler. Recalling that In forward Euler
\[ u^{n+1}=Bu^{n}\] Where \(B\) is the update matrix given by \(B=I+D\Delta tL_{x}\), where \(L_{x}\) is the 1-D Laplacian operator for \(u_{xx}\) with Dirichlet boundary conditions, with eigenvalue \(\lambda =\frac {2}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) \) where \(p=1\cdots N\) using the standard grid convention used before.
Let \(\mu \) be the eigenvalue of the above update matrix \(B\), hence\[ \mu =1+\frac {2D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) \] For stability, \(\left \vert u\right \vert \leq 1\) hence\begin {align*} \left \vert 1+\frac {2D\Delta t}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) \right \vert & \leq 1\\ \left \vert 1-4\frac {D\Delta t}{h^{2}}\right \vert & \leq 1 \end {align*}
Simplifying, this gives\[ 0\leq \frac {D\Delta t}{h^{2}}\leq \frac {1}{2}\] To compare the above with the R-K scheme in this problem, since\begin {align} \frac {u^{\ast }-u^{n}}{\Delta t} & =Lu^{n}\nonumber \\ \frac {u^{n+1}-u^{n}}{\Delta t} & =\frac {1}{2}\left ( Lu^{n}+Lu^{\ast }\right ) \tag {1} \end {align}
Expanding gives\begin {align*} u^{n+1} & =u^{n}+\frac {\Delta t}{2}\left ( Lu^{n}+L\left ( u^{n}+\Delta tLu^{n}\right ) \right ) \\ & =u^{n}+\frac {\Delta t}{2}\left ( Lu^{n}+Lu^{n}+\Delta tL^{2}u^{n}\right ) \\ & =\left ( I+\frac {\Delta t}{2}\left ( 2L+\Delta tL^{2}\right ) \right ) u^{n}\\ & =Bu^{n} \end {align*}
The eigenvalue of \(L\) is \(\frac {2}{h^{2}}\left ( \cos \left ( p\pi h\right ) -1\right ) \) but since \(\cos A-1=-2\sin ^{2}\frac {A}{2}\), hence \(\cos \left ( p\pi h\right ) -1=-2\sin ^{2}\frac {p\pi h}{2}\) and the eigenvalue is written as \(-\frac {4}{h^{2}}\sin ^{2}\frac {p\pi h}{2}\,.\)
let \(\sin ^{2}\frac {p\pi h}{2}\equiv \omega ,\) therefore the eigenvalue of \(L\) becomes \(\lambda =\) \(-\frac {4\omega }{h^{2}}\).
Using this notation, the above update matrix \(B\) for this scheme will have the following eigenvalue\begin {align*} \mu & =1+\frac {\Delta t}{2}\left ( -D\frac {8\omega }{h^{2}}+\Delta t\left ( -D\frac {4\omega }{h^{2}}\right ) ^{2}\right ) \\ & =1-\frac {4D\Delta t\omega }{h^{2}}+\Delta ^{2}t\frac {8D^{2}\omega ^{2}}{h^{4}} \end {align*}
For stability, \(\left \vert \mu \right \vert \leq 1.\) Maximum \(\mu \) will occur at \(\omega =1\) implying that \(\sin ^{2}\frac {p\pi h}{2}=1\) or \(p=N\), resulting in\[ \left \vert 1-\frac {4D\Delta t}{h^{2}}+\Delta ^{2}t\frac {8D^{2}}{h^{4}}\right \vert \leq 1 \] Therefore\begin {align*} -1 & \leq 1-\frac {4D\Delta t}{h^{2}}+\Delta ^{2}t\frac {8D^{2}}{h^{4}}\leq 1\\ -2 & \leq -\frac {4D\Delta t}{h^{2}}+\Delta ^{2}t\frac {8D^{2}}{h^{4}}\leq 0\\ 0 & \leq \frac {4D\Delta t}{h^{2}}-\Delta ^{2}t\frac {8D^{2}}{h^{4}}\leq 2\\ 0 & \leq \frac {D\Delta t}{h^{2}}-\Delta ^{2}t\frac {2D^{2}}{h^{4}}\leq \frac {1}{2}\\ 0 & \leq \frac {D\Delta t}{h^{2}}-2\left ( \frac {D\Delta t}{h^{2}}\right ) ^{2}\leq \frac {1}{2} \end {align*}
This shows that with the given RK scheme, stability implies the condition \(0\leq \frac {D\Delta t}{h^{2}}-\Delta ^{2}t\frac {2D^{2}}{h^{4}}\leq \frac {1}{2}.\) This is compared to \(0\leq \frac {D\Delta t}{h^{2}}\leq \frac {1}{2}\) for forward Euler.
What does this mean in terms of the time step? Will this allow the use of a larger time step than with FE while keep absolute stability?
Assume \(D=1\), This is a table showing the maximum value of \(\Delta t\) allowed for different \(h\) values
| scheme | \(h=1\) | \(h=0.1\) | \(h=0.01\) |
| FE \(0\leq \frac {\Delta t}{h^{2}}\leq \frac {1}{2}\) | \(0.5\) | \(0.05\) | \(0.005\) |
| RK \(0\leq \frac {\Delta t}{h^{2}}-\frac {2\Delta ^{2}t}{h^{4}}\leq \frac {1}{2}\) | \(0.5\) | \(0.05\) | \(0.005\) |
Hence, the largest time step does not change with this scheme when compared to forward Euler. It seems based on the above, that the explicit Runge-Kutta scheme for solving the diffusion PDE does not offer any advantage in handling the stiffness of the PDE since the time step remained constrained by \(h^{2}\) as with FE. \(\ \)It seems that explicit schemes are not suitable for stiff problems
Given\[ y_{t}=Dy_{xx}\] Then\begin {align*} \frac {y_{i}^{\ast }-y_{i}^{n}}{\Delta t} & =D\frac {y_{i-1}^{n}-2y_{i}^{n}+y_{i+1}^{n}}{h^{2}}\\ y_{i}^{\ast } & =y_{i}^{n}+\Delta t\overset {f\left ( y^{n}\right ) }{\overbrace {\frac {D}{h^{2}}\left ( y_{i-1}^{n}-2y_{i}^{n}+y_{i+1}^{n}\right ) }} \end {align*}
Hence\begin {align*} 2\frac {y_{i}^{n+1}-y_{i}^{n}}{\Delta t} & =\left ( f\left ( y^{n}\right ) +f\left ( y^{\ast }\right ) \right ) \\ & =\left [ \frac {D}{h^{2}}\left ( y_{i-1}^{n}-2y_{i}^{n}+y_{i+1}^{n}\right ) \right ] +\left [ \frac {D}{h^{2}}\left ( y_{i-1}^{\ast }-2y_{i}^{\ast }+y_{i+1}^{\ast }\right ) \right ] \\ & \\ & =\frac {D}{h^{2}}\left ( y_{i-1}^{n}-2y_{i}^{n}+y_{i+1}^{n}\right ) +\\ & \frac {D}{h^{2}}\left [ y_{i-1}^{n}+\Delta t\frac {D}{h^{2}}\left ( y_{i-2}^{n}-2y_{i-1}^{n}+y_{i}^{n}\right ) \right ] -\\ & 2\frac {D}{h^{2}}\left [ y_{i}^{n}+\Delta t\frac {D}{h^{2}}\left ( y_{i-1}^{n}-2y_{i}^{n}+y_{i+1}^{n}\right ) \right ] +\\ & \frac {D}{h^{2}}\left [ y_{i+1}^{n}+\Delta t\frac {D}{h^{2}}\left ( y_{i}^{n}-2y_{i+1}^{n}+y_{i+2}^{n}\right ) \right ] \end {align*}
Expand and simplify\begin {align*} 2\frac {y_{i}^{n+1}-y_{i}^{n}}{\Delta t} & =\frac {D}{h^{2}}y_{i-1}^{n}\left ( 1+1-2\Delta t\frac {D}{h^{2}}-2\Delta t\frac {D}{h^{2}}\right ) \\ & +\frac {D}{h^{2}}y_{i}^{n}\left ( -2+\Delta t\frac {D}{h^{2}}-2-4\Delta t\frac {D}{h^{2}}+\Delta t\frac {D}{h^{2}}\right ) \\ & +\frac {D}{h^{2}}y_{i+1}^{n}\left ( 1-2\Delta t\frac {D}{h^{2}}+1-2\Delta t\frac {D}{h^{2}}\right ) \\ & +\frac {D}{h^{2}}y_{i+2}^{n}\left ( \Delta t\frac {D}{h^{2}}\right ) \\ & \\ & =\frac {D}{h^{2}}\left [ y_{i-1}^{n}\left ( 2-4\Delta t\frac {D}{h^{2}}\right ) +y_{i}^{n}\left ( -4-2\Delta t\frac {D}{h^{2}}\right ) +y_{i+1}^{n}\left ( 2-4\Delta t\frac {D}{h^{2}}\right ) +y_{i+2}^{n}\left ( \Delta t\frac {D}{h^{2}}\right ) \right ] \end {align*}
Hence\begin {align*} y_{i}^{n+1} & =y_{i}^{n}+\frac {D\Delta t}{2h^{4}}\left [ y_{i-1}^{n}\left ( 2h^{2}-4\Delta tD\right ) +y_{i}^{n}\left ( -4h^{2}-2\Delta tD\right ) +y_{i+1}^{n}\left ( 2h^{2}-4\Delta tD\right ) +y_{i+2}^{n}\left ( \Delta tD\right ) \right ] \\ & =y_{i}^{n}+\frac {D^{2}\Delta ^{2}t}{2h^{4}}\left [ 2y_{i-1}^{n}\left ( \frac {h^{2}}{\Delta tD}-2\right ) -2y_{i}^{n}\left ( \frac {2h^{2}}{\Delta tD}+1\right ) +2y_{i+1}^{n}\left ( \frac {h^{2}}{\Delta tD}-2\right ) +y_{i+2}^{n}\right ] \\ & \\ & =y_{i-1}^{n}\left ( \frac {D^{2}\Delta ^{2}t}{2h^{4}}\left ( \frac {h^{2}}{\Delta tD}-2\right ) \right ) +y_{i}^{n}\left ( 1-\frac {D^{2}\Delta ^{2}t}{h^{4}}\left ( \frac {2h^{2}}{\Delta tD}+1\right ) \right ) +\\ & y_{i+1}^{n}\left ( \frac {D^{2}\Delta ^{2}t}{h^{4}}\left ( \frac {h^{2}}{\Delta tD}-2\right ) \right ) +y_{i+2}^{n}\frac {D^{2}\Delta ^{2}t}{2h^{4}} \end {align*}
Let \(\frac {D^{2}\Delta ^{2}t}{h^{4}}=r_{1}\) and let \(\frac {h^{2}}{\Delta tD}=r_{2}\) then4\begin {align*} y_{i}^{n+1} & =y_{i-1}^{n}\left ( \frac {r_{1}}{2}\left ( r_{2}-2\right ) \right ) +y_{i}^{n}\left ( 1-r_{1}\left ( 2r_{2}+1\right ) \right ) +y_{i+1}^{n}\left ( r_{1}\left ( r_{2}-2\right ) \right ) +y_{i+2}^{n}\frac {r_{1}}{2}\\ & =y_{i-1}^{n}\left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) +y_{i}^{n}\left ( 1-2r_{1}r_{2}-r_{1}\right ) +y_{i+1}^{n}\left ( r_{1}r_{2}-2r_{1}\right ) +y_{i+2}^{n}\frac {r_{1}}{2} \end {align*}
Using this 1-D grid
The matrix form of the scheme becomes\begin {align*} \begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N-1}^{n+1}\\ u_{N}^{n+1}\end {bmatrix} & =\begin {bmatrix} \left ( 1-2r_{1}r_{2}-r_{1}\right ) & \left ( r_{1}r_{2}-2r_{1}\right ) & \frac {r_{1}}{2} & 0 & 0 & 0\\ \left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) & \left ( 1-2r_{1}r_{2}-r_{1}\right ) & \left ( r_{1}r_{2}-2r_{1}\right ) & \frac {r_{1}}{2} & 0 & 0\\ 0 & \left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) & \left ( 1-2r_{1}r_{2}-r_{1}\right ) & \left ( r_{1}r_{2}-2r_{1}\right ) & \frac {r_{1}}{2} & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) & \left ( 1-2r_{1}r_{2}-r_{1}\right ) & \left ( r_{1}r_{2}-2r_{1}\right ) \\ 0 & 0 & 0 & 0 & \left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) & \left ( 1-2r_{1}r_{2}-r_{1}\right ) \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N}^{n}\end {bmatrix} +\\ & \begin {bmatrix} \left ( \frac {r_{1}r_{2}}{2}-r_{1}\right ) u_{0}^{n}\\ 0\\ 0\\ \vdots \\ y_{N+1}^{n}\frac {r_{1}}{2}\\ \left ( r_{1}r_{2}-2r_{1}\right ) u_{N+1}^{n}\end {bmatrix} \end {align*}
There is a problem above at node \(N\) as two additional nodes are needed to its right, but only one node exist. Need to look more into this later, as this part is not required for this HW.
The PDE is \(u_{t}+au_{x}=bu_{xx}\), \(\ \)with \(b>0\), the forward time, centered space discretization is\begin {equation} \frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+a\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h}=b\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}} \tag {1} \end {equation}
Applying Von Neumann analysis, let \(u_{j}^{n}=e^{i\xi x_{j}}\) and \(u_{j}^{n+1}=g\left ( \xi \right ) e^{i\xi x_{j}}\), then (1) becomes\begin {align*} \frac {ge^{i\xi x_{j}}-e^{i\xi x_{j}}}{\Delta t}+a\frac {e^{i\xi x_{j}}e^{i\xi h}-e^{i\xi x_{j}}e^{-i\xi h}}{2h} & =b\frac {e^{i\xi x_{j}}e^{-i\xi h}-2e^{i\xi x_{j}}+e^{i\xi x_{j}}e^{i\xi h}}{h^{2}}\\ ge^{i\xi x_{j}}-e^{i\xi x_{j}}+\frac {a\Delta t}{2h}\left ( e^{i\xi x_{j}}e^{i\xi h}-e^{i\xi x_{j}}e^{-i\xi h}\right ) & =\frac {b\Delta t}{h^{2}}\left ( e^{i\xi x_{j}}e^{-i\xi h}-2e^{i\xi x_{j}}+e^{i\xi x_{j}}e^{i\xi h}\right ) \\ g-1+\frac {\nu }{2}\left ( e^{i\xi h}-e^{-i\xi h}\right ) & =\mu \left ( e^{-i\xi h}-2+e^{i\xi h}\right ) \\ g\left ( \xi \right ) & =1-\frac {\nu }{2}\left ( e^{i\xi h}-e^{-i\xi h}\right ) +\mu \left ( e^{-i\xi h}-2+e^{i\xi h}\right ) \end {align*}
Hence5\[ g\left ( \xi \right ) =1+2\mu \left ( \cos \left ( \xi h\right ) -1\right ) -i\nu \sin \left ( \xi h\right ) \] But \(\cos A-1=-2\sin ^{2}\frac {A}{2}\), hence \(\cos \left ( \xi h\right ) -1=-2\sin ^{2}\frac {\xi h}{2}\) and the above becomes\[ g\left ( \xi \right ) =1-4\mu \sin ^{2}\left ( \frac {\xi h}{2}\right ) -i\nu \sin \left ( \xi h\right ) \] Therefore \begin {equation} \left \vert g\left ( \xi \right ) \right \vert ^{2}=\left [ 1-4\mu \sin ^{2}\left ( \frac {\xi h}{2}\right ) \right ] ^{2}+\nu ^{2}\sin ^{2}\left ( \xi h\right ) \tag {2} \end {equation} Using the trig identity \[ \sin ^{2}\left ( \xi h\right ) =4\sin ^{2}\left ( \frac {\xi h}{2}\right ) \left ( 1-\sin ^{2}\left ( \frac {\xi h}{2}\right ) \right ) \] then (2) becomes\[ \left \vert g\left ( \xi \right ) \right \vert ^{2}=\left [ 1-4\mu \sin ^{2}\left ( \frac {\xi h}{2}\right ) \right ] ^{2}+4\nu ^{2}\sin ^{2}\left ( \frac {\xi h}{2}\right ) \left ( 1-\sin ^{2}\left ( \frac {\xi h}{2}\right ) \right ) \] Let \(\sin ^{2}\left ( \frac {\xi h}{2}\right ) \equiv \omega \) in the above\[ \left \vert g\left ( \xi \right ) \right \vert ^{2}=\left ( 1-4\mu \omega \right ) ^{2}+4\nu ^{2}\omega \left ( 1-\omega \right ) \] The maximum of \(\left \vert g\left ( \xi _{p}\right ) \right \vert \) occurs when \(\xi h\approx \pi \), making \(\omega =1\), hence from above, the maximum of \(\left \vert g\left ( \xi \right ) \right \vert ^{2}\) is reduced to \(1-4\mu \), and then for stability\[ \left \vert 1-4\mu \right \vert \leq 1 \] therefore\begin {align*} -1 & \leq 1-4\mu \leq 1\\ -2 & \leq -4\mu \leq 0\\ 0 & \leq 4\mu \leq 2 \end {align*}
Hence\[ 0\leq \mu \leq \frac {1}{2} \] The above result can also be derived using the stencil diagram method. The stencil diagram for the above scheme (for internal nodes only) is
As was done in problem 1, By imposing that the weight on each edge in the above directed graph not exceed unity, and that the total algebraic sum of the weight of the edge also not exceed unity. This includes any combination of edges involved. For if this was the case, then \(u_{j}^{n+1}\) will always have an amplitude \(\leq \) \(u_{j}^{n}\). Applying this to the above diagram gives \[ \left \{ \begin {array} [c]{lll}\left ( 1\right ) & \frac {\upsilon }{2}+\mu \leq 1 & \text {Condition on }j-1\text { edge}\\ \left ( 2\right ) & \left \vert \mu -\frac {\upsilon }{2}\right \vert \leq 1 & \text {Condition on }j+1\text { edge}\\ (3) & 2\mu \leq 1 & \text {Condition on }j-1\text { added to }j+1\text { edge}\\ (4) & \left \vert 1-2\mu \right \vert \leq 1 & \text {Condition on the }j\text { edge}\\ \left ( 5\right ) & \left \vert 1+\frac {\upsilon }{2}-\mu \right \vert \leq 1 & \text {Condition on the }j\text { edge added to }j-1\\ \left ( 6\right ) & \left \vert 1-\frac {\upsilon }{2}-\mu \right \vert \leq 1 & \text {Condition on the }j\text { edge added to }j+1\\ \left ( 7\right ) & 1\leq 1 & \text {Condition that all edges sum to less than 1}\end {array} \right . \] Condition (7) gives no information. Condition (4) gives \(\mu \leq 1\) and hence weaker than (3), condition (5) is the same as (2) and gives \(\mu -\frac {\upsilon }{2}\leq 1\), condition (6) gives \(\frac {\upsilon }{2}+\mu \leq 2\) which is weaker than condition (1). Hence the following are remaining conditions (3),(1),(2).\[ \left \{ \begin {array} [c]{lll}\left ( 1\right ) & \frac {\upsilon }{2}+\mu \leq 1 & \text {Condition on }j-1\text { edge}\\ \left ( 2\right ) & \left \vert \mu -\frac {\upsilon }{2}\right \vert \leq 1 & \text {Condition on }j+1\text { edge}\\ (3) & 2\mu \leq 1 & \text {Condition on }j-1\text { added to }j+1\text { edge}\end {array} \right . \]
From the above, only condition (3) can provide useful information, which is that \(\mu \leq \frac {1}{2}\), which is what was found using Von Neumann analysis.
The scheme was implemented. The source code is in the appendix of this problem. The grid used is the standard grid
In the following, I will be use the following PDE \(cu_{t}+au_{x}=bu_{xx}\) (where\(c\) was added a parameter for the advection term).
Periodic boundary conditions implies \(u\left ( 0,t\right ) =u\left ( 1,t\right ) \), Hence there is an extra one unknown (in addition to the internal nodes). Either \(u\left ( 0,t\right ) \) or \(u\left ( 1,t\right ) \) can be selected since they have the same value. When selecting the right end node, then \(u_{N+1}\) becomes an unknown to be added to the internal nodes. Using the following diagram
The forward time, centered space discretization is\begin {equation} c\frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+a\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h}=b\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}\tag {1} \end {equation}
Therefore, for nodes \(2\cdots N\), the finite difference scheme is \begin {align*} c\frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+a\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h} & =b\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}\\ u_{j}^{n+1} & =u_{j}^{n}-a\Delta t\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2ch}+b\Delta t\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{ch^{2}}\\ u_{j}^{n+1} & =u_{j}^{n}-\frac {a\Delta t}{2ch}\left ( u_{j+1}^{n}-u_{j-1}^{n}\right ) +\frac {b\Delta t}{ch^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) \\ u_{j}^{n+1} & =u_{j-1}^{n}\left ( \frac {a\Delta t}{2ch}+\frac {b\Delta t}{ch^{2}}\right ) +u_{j}^{n}\left ( 1-2\frac {b\Delta t}{ch^{2}}\right ) +u_{j+1}^{n}\left ( -\frac {a\Delta t}{2ch}+\frac {b\Delta t}{ch^{2}}\right ) \end {align*}
Let \(\upsilon =\frac {a\Delta t}{ch}\), \(\mu =\frac {b\Delta t}{ch^{2}},\) the above becomes
\[ u_{j}^{n+1}=u_{j-1}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{j}^{n}\left ( 1-2\mu \right ) +u_{j+1}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \] Node \(j=1\) gives\begin {align*} u_{1}^{n+1} & =u_{0}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{1}^{n}\left ( 1-2\mu \right ) +u_{2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \\ & =u_{N+1}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{1}^{n}\left ( 1-2\mu \right ) +u_{2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \end {align*}
And for node \(j=N+1\)\begin {align*} u_{N+1}^{n+1} & =u_{N}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{N+1}^{n}\left ( 1-2\mu \right ) +u_{N+2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \\ & =u_{N}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{N+1}^{n}\left ( 1-2\mu \right ) +u_{1}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \end {align*}
The full system can now be written in matrix form\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N}^{n+1}\\ u_{N+1}^{n+1}\end {bmatrix} =\begin {bmatrix} \left ( 1-2\mu \right ) & \left ( \mu -\frac {\upsilon }{2}\right ) & 0 & 0 & 0 & \left ( \frac {\upsilon }{2}+\mu \right ) \\ \left ( \mu +\frac {\upsilon }{2}\right ) & \left ( 1-2\mu \right ) & \left ( \mu -\frac {\upsilon }{2}\right ) & 0 & 0 & 0\\ 0 & \left ( \mu +\frac {\upsilon }{2}\right ) & \left ( 1-2\mu \right ) & \left ( \mu -\frac {\upsilon }{2}\right ) & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 0 & 0 & 0 & \left ( \mu +\frac {\upsilon }{2}\right ) & \left ( 1-2\mu \right ) & \left ( \mu -\frac {\upsilon }{2}\right ) \\ \left ( \mu -\frac {\upsilon }{2}\right ) & 0 & 0 & 0 & \left ( \mu +\frac {\upsilon }{2}\right ) & \left ( 1-2\mu \right ) \end {bmatrix}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N}^{n}\\ u_{N+1}^{n}\end {bmatrix} \] The above can also be written as
\[\begin {bmatrix} u_{1}^{n+1}\\ u_{2}^{n+1}\\ u_{3}^{n+1}\\ \vdots \\ u_{N}^{n+1}\\ u_{N+1}^{n+1}\end {bmatrix} =\overset {update\ matrix}{\overbrace {\left ( I-\frac {a\Delta t}{2ch}\begin {bmatrix} 0 & 1 & 0 & 0 & 0 & -1\\ -1 & 0 & 1 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 & 0 & 0\\ \vdots & \vdots & -1 & \ddots & 1 & 0\\ 0 & 0 & 0 & -1 & 0 & 1\\ 1 & 0 & 0 & 0 & -1 & 0 \end {bmatrix} +\frac {b\Delta t}{ch^{2}}\begin {bmatrix} -2 & 1 & 0 & 0 & 0 & 1\\ 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & 0\\ 0 & 0 & 0 & 1 & -2 & 1\\ 1 & 0 & 0 & 0 & 1 & -2 \end {bmatrix} \right ) }}\begin {bmatrix} u_{1}^{n}\\ u_{2}^{n}\\ u_{3}^{n}\\ \vdots \\ u_{N}^{n}\\ u_{N+1}^{n}\end {bmatrix} \] \(u_{j}^{0}\) are taken from initial conditions. The above is in the form\[ u^{n+1}=Bu^{n}\] and is implemented directly as above in the code. Using the numerical values given in the problem\[ \Delta t=\frac {0.25h^{2}}{b}=0.25\frac {\left ( 0.05\right ) ^{2}}{1}=0.000\allowbreak 6\,25 \] and \[ \upsilon =\frac {a\Delta t}{h}=\frac {80\left ( 0.000\allowbreak 6\,25\right ) }{0.05}=1 \] and \[ \mu =\frac {b\Delta t}{h^{2}}=\frac {0.000\allowbreak 6\,25}{\left ( 0.05\right ) ^{2}}=0.25 \] since \(\mu \leq \frac {1}{2}\) the solution is expected to be stable since this is the condition derived in part (a).
The following is the result of running the program. The numerical solution grew with time. Here are few snap shots taken at increasing time steps showing the problem. After only about 10 time steps, the numerical solution can be seen to grow more than the initial conditions
The stability analysis that was done did not predict this based on the value of \(\mu \) which was \(\leq \frac {1}{2}\).
The eigenvalues \(u_{p}\) of the update matrix derived in part (b) are\[ u_{p}=1-i\frac {a\Delta t}{h}\sin \left ( \pi ph\right ) +\frac {2b\Delta t}{h^{2}}\left ( \cos \left ( \pi ph\right ) -1\right ) \] But \(\upsilon =\frac {a\Delta t}{h}\), \(\mu =\frac {b\Delta t}{h^{2}}\), then the above becomes\[ u_{p}=1-iv\sin \left ( \pi ph\right ) +2\mu \left ( \cos \left ( \pi ph\right ) -1\right ) \] Let \(\sin ^{2}\left ( \frac {\xi h}{2}\right ) \equiv \omega \), then the above can be written as\[ \left \vert u_{p}\right \vert ^{2}=\left ( 1-4\mu \omega \right ) ^{2}+4\nu ^{2}\omega \left ( 1-\omega \right ) \] As was done in part (a).
From part (a), it was found that when \(\omega =1\), this resulted in the condition of stability being \(0\) \(\leq \mu \leq \frac {1}{2}\), and when \(\omega =0\), the maximum eigenvalue is \(1.\) To find the condition of \(\left \vert u_{p}\right \vert \leq 1\) for the full range of \(\omega \), first expand \(\left \vert u_{p}\right \vert ^{2}\) into a quadratic in \(\omega \) and minimize \begin {align*} \left \vert u_{p}\right \vert ^{2} & =1+16\mu ^{2}\omega ^{2}-8\mu \omega +4\nu ^{2}\omega -4\nu ^{2}\omega ^{2}\\ & =1-4\omega \left ( 2\mu -\nu ^{2}\right ) +4\omega ^{2}\left ( 4\mu ^{2}-\nu ^{2}\right ) \end {align*}
Since \(\omega =\sin ^{2}\left ( \frac {\pi ph}{2}\right ) \), then \(\omega \) values are from \(0\cdots 1.\) Since the maximum eigenvalue occurs when \(\omega =0\), then \(\omega =0\) is the maximum point of the quadratic \(1-4\omega \left ( 2\mu -\nu ^{2}\right ) +4\omega ^{2}\left ( 4\mu ^{2}-\nu ^{2}\right ) \), hence the slope of this quadratic at \(\omega =0\) must be negative. But the slope is \begin {align*} \frac {d}{d\omega }\left \vert u_{p}\right \vert ^{2} & =\frac {d}{d\omega }\left ( 1-4\omega \left ( 2\mu -\nu ^{2}\right ) +4\omega ^{2}\left ( 4\mu ^{2}-\nu ^{2}\right ) \right ) \\ & =-4\left ( 2\mu -\nu ^{2}\right ) +8\omega \left ( 4\mu ^{2}-\nu ^{2}\right ) \end {align*}
For the above to be negative (so that the eigenvalue remain below 1) implies that \(2\mu -\nu ^{2}\) must be positive, i.e. \[ 2\mu -\nu ^{2}\geq 0 \] or\[ \nu ^{2}\leq 2\mu \] And since from part(a) it was found that \(\mu \leq \frac {1}{2}\), or \(2\mu \leq 1\) then the above become\[ \nu ^2\leq 2\mu \leq 1 \]
The solution for \(\Delta t=0.000\allowbreak 6\,25\), \(h=0.05\), \(a=80\), \(b=1\), with \(u\left ( x,0\right ) =\exp \left ( -20\left ( x-0.5\right ) ^{2}\right ) \) at \(t=0.01\) seconds is
The value of the numerical solution at that time is
-0.1727 0.1817 0.6603 1.1158 1.4092 1.4609 1.2868 0.9910 0.6912 0.4441 0.2527 0.1168 0.0451 0.0337 0.0532 0.0563 -0.0013 -0.1275 -0.2680 -0.3162 -0.1727
derivation of the convection-diffusion using general terms The PDE is \(cu_{t}+au_{x}=bu_{xx}\), \(\ \)with \(b>0\), the forward time, centered space discretization is\begin {equation} c\frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+a\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h}=b\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}\tag {1} \end {equation}
Therefore, for nodes \(2\cdots N\), the finite difference scheme is \begin {align*} c\frac {u_{j}^{n+1}-u_{j}^{n}}{\Delta t}+a\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2h} & =b\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{h^{2}}\\ u_{j}^{n+1} & =u_{j}^{n}-a\Delta t\frac {u_{j+1}^{n}-u_{j-1}^{n}}{2ch}+b\Delta t\frac {u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}}{ch^{2}}\\ u_{j}^{n+1} & =u_{j}^{n}-\frac {a\Delta t}{2ch}\left ( u_{j+1}^{n}-u_{j-1}^{n}\right ) +\frac {b\Delta t}{ch^{2}}\left ( u_{j-1}^{n}-2u_{j}^{n}+u_{j+1}^{n}\right ) \\ u_{j}^{n+1} & =u_{j-1}^{n}\left ( \frac {a\Delta t}{2ch}+\frac {b\Delta t}{ch^{2}}\right ) +u_{j}^{n}\left ( 1-2\frac {b\Delta t}{ch^{2}}\right ) +u_{j+1}^{n}\left ( -\frac {a\Delta t}{2ch}+\frac {b\Delta t}{ch^{2}}\right ) \end {align*}
Let \(\upsilon =\frac {a\Delta t}{ch}\), \(\mu =\frac {b\Delta t}{ch^{2}},\) the above becomes\[ u_{j}^{n+1}=u_{j-1}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{j}^{n}\left ( 1-2\mu \right ) +u_{j+1}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \] Node \(j=1\) gives\begin {align*} u_{1}^{n+1} & =u_{0}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{1}^{n}\left ( 1-2\mu \right ) +u_{2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \\ & =u_{N+1}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{1}^{n}\left ( 1-2\mu \right ) +u_{2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \end {align*}
And for node \(j=N+1\)\begin {align*} u_{N+1}^{n+1} & =u_{N}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{N+1}^{n}\left ( 1-2\mu \right ) +u_{N+2}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \\ & =u_{N}^{n}\left ( \frac {\upsilon }{2}+\mu \right ) +u_{N+1}^{n}\left ( 1-2\mu \right ) +u_{1}^{n}\left ( \mu -\frac {\upsilon }{2}\right ) \end {align*}
1This is slightly different from the standard numbering format we used before.
2The magnification factor is the term \(g\left ( \xi \right ) \) in the expression relating \(\hat {u}^{n+1}\) to \(\hat {u}^{n}\) in the expression \(\hat {u}^{n+1}=g\left ( \xi \right ) \) \(\hat {u}^{n}\)
3It is possible to derive the amplification factor using direct application of fourier transform, but the procedure is longer. The final result will be the same. The appendix of this problem contain this derivation.
4Notice that units of \(D\) are \(meter^{2}\) per second, hence \(\frac {h^{2}}{\Delta tD}\) is dimensionless, so we are ok.
5Notice that the first oder derivatives (or odd order in general) produces eigenvalues that are complex, and the even order ones produce real eigenvalues.