Use standard 3-point discretization of the Laplacian on a regular mesh to find numerical solution to the PDE below and perform refinement study to compute the error that shows the rate of convergence for both the 1-norm and the max-norm.
The differential equation with its boundary conditions is \(u_{xx}=e^{x}\) with \(u\left ( 0\right ) =0,u\left ( 1\right ) =1.\)
Finding the analytical solution
Let \(D\) be a differential operator \(D\equiv \frac{d}{dx}.\) Since \(\left ( D-1\right ) e^{x}=0\), then applying \(\left ( D-1\right ) \) to both sides of the differential equation results in\[ \left [ \left ( D-1\right ) D^{2}\right ] (u)=0 \] Thus the characteristic equation is \(\left ( r-1\right ) r^{2}=0\) with roots \(r_{1}=1,r_{2}=0\,\) and \(r_{3}=0\). Therefore the complete solution is \[ u\left ( x\right ) =c_{1}e^{r_{1}x}+c_{2}e^{r_{2}x}+xc_{3}e^{r_{3}x}\] It helps to designate in the above which part is the particular solution and which is the homogeneous\[ u\left ( x\right ) =\overset{u_{p}}{\overbrace{c_{1}e^{x}}}+\overset{y_{h}}{\overbrace{c_{2}+c_{3}x}}\] \(c_{1}\) is found by substituting the particular solution in the original differential equation giving \(c_{1}=1\). The solution becomes1 \[ u\left ( x\right ) =e^{x}+c_{2}+xc_{3}\] The remaining constants \(c_{2}\) and \(c_{3}\) are found by satisfying the boundary conditions on the above solution. Applying \(u\left ( 0\right ) =0\,\ \)gives \(c_{2}=-1\) and applying \(u\left ( 1\right ) =1\) gives \(c_{3}=2-e\). The solution becomes \[ u\left ( x\right ) =2x-\exp \left ( 1\right ) x+\exp \left ( x\right ) -1 \]
Scheme to use for the numerical solution
The numbering used is the one described in the class. The first grid point has an index \(j=0\), and the last grid point has an index \(j=N+1\). \(U\) is the unknown to solve for. It is the numerical solution of the differential equation at the grid points. Lower case \(u\) is the exact analytical solution found earlier.
Because this problem has Dirichlet boundary conditions, \(U\) is known at \(j=0\) and at \(j=N+1\), therefore only internal grid points are used to solve for \(U\). These internal points are numbered \(1\cdots N\). The spacing between each grid point is \(h=\frac{len}{N+1}\) where the length of the domain \(len\) is always taken to have value \(1\).
The physical x-coordinate of each grid point is given by \(x_{j}=jh\). For example, when \(j=0\), the first point will have a physical x-coordinate of 0, and when \(j=N+1\), the last grid point will have a physical x-coordinate of \(1\).
The diagram below helps illustrate these relations and the notations used.
The standard 3-point central difference formula \(\frac{U_{j-1}-2U_{j}+U_{j+1}}{h^{2}}\) was used to approximate \(u_{xx}\) on each internal grid point. This approximation has local truncation error (LTE) of \(O\left ( h^{2}\right ) \). Therefore, at each internal grid point \(j\) the equation \(u_{xx}=f\left ( x\right ) \) is approximated by \(\frac{U_{j-1}-2U_{j}+U_{j+1}}{h^{2}}=f_{x=jh}\).
For the special case of \(j=1\) (the first internal grid point on left side), the above formula was modified by replacing \(U_{o}\) by its given value (\(U_{0}\) is on the boundary and its value is known). Therefore, for \(j=1\) the discrete equation becomes\begin{align} \frac{\alpha -2U_{0}+U_{1}}{h^{2}} & =F_{1}\nonumber \\ \frac{-2U_{0}+U_{1}}{h^{2}} & =F_{1}-\frac{\alpha }{h^{2}} \tag{2} \end{align}
The same was done for \(j=N\) (the last internal grid point on right side). The standard 3 points formula was modified by replacing \(U_{N+1}\) by its given value. Therefore, for \(j=N\) the discrete equation becomes\begin{align} \frac{U_{N-1}-2U_{N}+\beta }{h^{2}} & =F_{N}\nonumber \\ \frac{U_{N-1}-2U_{N}}{h^{2}} & =F_{N}-\frac{\beta }{h^{2}} \tag{3} \end{align}
The above equations are collected and put in the form \(Au=f\) resulting in \[ \frac{1}{h^{2}}\begin{pmatrix} -2 & 1 & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & \cdot & \cdot & \cdot & 0\\ 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix}\begin{pmatrix} U_{1}\\ U_{2}\\ U_{3}\\ \cdot \\ \cdot \\ \cdot \\ U_{N}\end{pmatrix} =\begin{pmatrix} f_{h}-\frac{\alpha }{h^{2}}\\ f_{2h}\\ f_{3h}\\ \cdot \\ \cdot \\ f_{\left ( N-1\right ) h}\\ F_{Nh}-\frac{\beta }{h^{2}}\end{pmatrix} \] The above system is solved for the unknown vector \(U\). These are the values of the numerical solution at the internal grid points\(.\) The \(A\) matrix above is symmetric and tridiagonal and of size \(N\times N\) where \(N\) is the number of internal grid points.
Error norm calculation
The total error at a grid point \(j\) is given by \[ \mathbf{e}_{j}\mathbf{=U}_{j}\mathbf{-u}\left ( x_{j}\right ) \] Where \(\mathbf{U}_{j}\) is the numerical solution at the \(j^{th}\) grid point and \(\mathbf{u}\left ( x_{j}\right ) \) is the exact solution sampled at the same grid point location. \(\mathbf{e}\) is a vector of length \(N\) and represents the difference between \(U\) and \(u\) at each point.
To measure the size of the error vector \(\mathbf{e}\), a grid norm is used in place of the standard vector norm. The following are the definitions of the norms used.
Description of method used in Refinement study
The goal of the refinement study is to check that the scheme selected is stable. Since the scheme selected is consistent (it has an LTE of order \(O\left ( h^{2}\right ) \) therefore \(\left \Vert \tau \right \Vert \rightarrow 0\) as \(h\rightarrow 0\)), this implies that the scheme will be stable if it can be shown that the numerical solution converges to the exact solution2 .
For a stable scheme the following relation must hold\[ \left \Vert e\right \Vert \leq \left \Vert A^{-1}\right \Vert \left \Vert \tau \right \Vert \] therefore stability can be established by verifying that error norm \(\left \Vert e\right \Vert \) is also of order \(O\left ( h^{2}\right ) \). This implies that \(\left \Vert A^{-1}\right \Vert \) is \(O\left ( 1\right ) \). In addition, showing that \(\left \Vert e\right \Vert \) is of order \(O\left ( h^{2}\right ) \) implies that \(\left \Vert e\right \Vert \rightarrow 0\) as \(h\rightarrow 0\), which means convergence.
These are the steps carried out in the refinement study
Results of the refinement study
This is a plot of \(\log (h)\) against the log of the various error norms. The source code is in the appendix.
The error table is the following
Conclusion
The study was carried out using 3 different norms (max norm,1-norm and 2-norm). It is found that for each norm the total error ratio converged to \(4\), implying \(p=2\), hence the total error norm was \(O\left ( h^{2}\right ) \) the same as LTE. This shows that the numerical scheme is stable and the numerical solution converges to the exact solution.
Solve \(u_{xx}=2\cos ^{2}\left ( \pi x\right ) \) with \(u_{x}\left ( 0\right ) =0,u_{x}\left ( 1\right ) =1\)
The existence of a solution is first verified. This is done in this case since for Neumann boundary conditions a solution might not even exist if the problem is not well posed. Integrating both sides of the PDE gives\begin{align*}{\displaystyle \int \limits _{0}^{1}} u_{xx}\left ( x\right ) dx & ={\displaystyle \int \limits _{0}^{1}} 2\cos ^{2}\left ( \pi x\right ) dx\\ & =u_{x}\left ( 1\right ) -u_{x}\left ( 0\right ) \\ & =1 \end{align*}
Substituting the values of \(u_{x}\)on the boundaries, the above becomes zero, hence the problem is well posed. Now the analytical solution is found.
Integrating the differential equation once\begin{equation} u_{x}=x+\frac{\sin \left ( 2\pi x\right ) }{2\pi }+c_{1} \tag{1} \end{equation} And Integrating again\begin{equation} u\left ( x\right ) =\frac{x^{2}}{2}-\frac{\cos \left ( 2\pi x\right ) }{4\pi ^{2}}+c_{1}x+c_{2} \tag{2} \end{equation} Substituting the boundary condition \(u_{x}\left ( 0\right ) =0\) in the above3 gives \(c_{1}=0\), therefore the solution is\begin{equation} u\left ( x\right ) =\frac{x^{2}}{2}-\frac{\cos \left ( 2\pi x\right ) }{4\pi ^{2}}+c_{2} \tag{3} \end{equation} This solution is not unique. The constant \(c_{2}\) is arbitrary and an infinite number of solutions exist. A solution exist which is up to an arbitrary additive constant. To select a constant for the purpose of the numerical analysis in the refinement study, the constant is selected to give the solution zero mean. Hence\begin{align*} \int _{0}^{1}u\left ( x\right ) dx & =0\\ \int _{0}^{1}\left ( \frac{x^{2}}{2}-\frac{\cos \left ( 2\pi x\right ) }{4\pi ^{2}}+c_{2}\right ) dx & =0\\ \frac{1}{6}+c_{2} & =0 \end{align*}
Therefore \(c_{2}=-\frac{1}{6}\). Therefore, the exact solution used in the refinement study to compare the numerical solution against is \[ u\left ( x\right ) =\frac{x^{2}}{2}-\frac{\cos \left ( 2\pi x\right ) }{4\pi ^{2}}-\frac{1}{6}\]
Scheme to use for the numerical solution
Two schemes formulated, both having a local truncation error \(\left \Vert \tau \right \Vert \) of order \(O\left ( h^{2}\right ) \). Hence both schemes are consistent.
First scheme The standard 3 point centered difference formula\(\frac{U_{j-1}-2U_{j}+U_{j+1}}{h^{2}}\) is used for approximating \(u_{xx}\) on internal grid points. However in this problem \(U_{o}\) and \(U_{N+1}\)are not known at the left most and the right most grid points (grid points of index \(0\) and \(N+1\)), therefore the grid points used includes these 2 points in addition to the standard internal grid points \((1\cdots N)\) resulting in the \(A\) matrix having size \(N+2.\) (\(A\) will have 2 additional rows as compared to part (a)).
Now, two imaginary points are introduced into the domain, one to the left of \(j=0\), and one to the right of \(j=N+1\)
For \(j=0\) the standard 3 points formulate is used to approximate \(u_{xx}\)\begin{equation} \frac{U_{-1}-2U_{0}+U_{1}}{h^{2}}=f_{0} \tag{1} \end{equation} And \(u_{x}\left ( 0\right ) =\alpha \) the 2 points central difference scheme is used\begin{equation} \frac{U_{1}-U_{-1}}{2h}=\alpha \tag{2} \end{equation} Now \(U_{-1}\) is eliminated using (1) and (2). From equation (2) \[ U_{-1}=U_{1}-2h\alpha \] substituting the above into (1)\begin{align} \frac{\left ( U_{1}-2h\alpha \right ) -2U_{0}+U_{1}}{h^{2}} & =f_{0}\nonumber \\ \frac{-2U_{0}+2U_{1}}{h^{2}} & =f_{0}+\frac{2\alpha }{h} \tag{3} \end{align}
The above equation (3) is the discrete equation for node \(j=0\) (the first row in the \(A\) matrix). Similarly for the right-most node \(j=N+1\), \(u_{xx}\) is approximated using\begin{equation} \frac{U_{N}-2U_{N+1}+U_{N+2}}{h^{2}}=f_{N+1} \tag{4} \end{equation} And \(u_{x}\left ( 1\right ) =\beta \) at node \(j=N+1\) is approximated using 2 points central difference\begin{equation} \frac{U_{N+2}-U_{N}}{2h}=\beta \tag{5} \end{equation} \(U_{N+2}\)is eliminated using (4) and (5). From equation (5) \[ U_{N+2}=U_{N}+2h\beta \] substituting the above into (4)\begin{align} \frac{U_{N}-2U_{N+1}+\left ( U_{N}+2h\beta \right ) }{h^{2}} & =f_{N+1}\nonumber \\ \frac{2U_{N}-2U_{N+1}}{h^{2}} & =f_{N+1}-\frac{2\beta }{h} \tag{6} \end{align}
The above equation (6) is the discrete equation for node \(j=N+1\) (the last row in the \(A\) matrix). \(Au=f\) is now setup resulting in \[ \frac{1}{h^{2}}\begin{pmatrix} -2 & 2 & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & \cdot & \cdot & \cdot & 0\\ 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 2 & -2 \end{pmatrix}\begin{pmatrix} U_{0}\\ U_{1}\\ U_{2}\\ \cdot \\ \cdot \\ N_{N-1}\\ U_{N}\\ U_{N+1}\end{pmatrix} =\begin{pmatrix} f_{0}+\frac{2\alpha }{h}\\ f_{1}\\ f_{2}\\ \cdot \\ \cdot \\ f_{N}\\ f_{N+1}-\frac{2\beta }{h}\end{pmatrix} \] The \(A\) matrix is tridiagonal, not symmetrical and singular since \(null(A)=\mathbf{1}^{T}\). In some books, this matrix is called the Laplacian matrix.
Second scheme In this scheme, the first order derivative for the Neumann boundary condition at the left point and at the right point is approximated using\begin{align*} \alpha & =\frac{1}{h}\left ( \frac{3}{2}U_{0}-2U_{1}+\frac{1}{2}U_{2}\right ) \\ \beta & =\frac{1}{h}\left ( \frac{3}{2}U_{N+1}-2U_{N}+\frac{1}{2}U_{N-1}\right ) \end{align*}
respectively.
The derivation of the above formulas is shown in example 1.2 in the textbook. The above approximation has an LTE of \(O\left ( h^{2}\right ) \). For the other internal grid points, the standard 3 points centred difference \(\frac{U_{j-1}-2U_{j}+U_{j+1}}{h^{2}}\) is used to approximate \(u_{xx}\). Using the above approximations, the system \(AU=f\) becomes\[ \frac{1}{h^{2}}\begin{pmatrix} \frac{3}{2} & -2 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & \cdot & \cdot & \cdot & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & \frac{1}{2} & -2 & \frac{3}{2}\end{pmatrix}\begin{pmatrix} U_{0}\\ U_{1}\\ U_{2}\\ \cdot \\ \cdot \\ U_{N-1}\\ U_{N}\\ U_{N+1}\end{pmatrix} =\begin{pmatrix} \frac{\alpha }{h}\\ F_{1}\\ F_{2}\\ \cdot \\ \cdot \\ \cdot \\ F_{N}\\ \frac{\beta }{h}\end{pmatrix} \] The \(A\) matrix is tridiagonal, not symmetrical and is also singular since \(null(A)=\mathbf{1}^{T}\).
Augmenting the systems before solving The \(A\) matrix for both schemes above is singular. Both have \(u=\left ( 1,1,1,\cdots ,1,1\right ) ^{T}\) as the null vector. Using Matlab to verify, the null command can be used as follows
For the second scheme
Since the matrix \(A\) is singular, before solving by Gaussian elimination it is augmented as follows\[ A\rightarrow \begin{pmatrix} A & v\\ u & 0 \end{pmatrix} \] where \(v\) is the the null of \(A^{\ast }\) (the adjoint of \(A,\) or for a real matrix, the same as \(A^{T}\)) and \(u\) is the null of \(A.\)
\(v\) is found for the above 2 schemes. For the first scheme:
And for the second scheme
Now that \(v\) and \(u\) are found, the augmented \(A\) can be formulated.
The force vector is also augmented as follows\begin{equation} f\rightarrow \begin{pmatrix} f\\ 0 \end{pmatrix} \tag{7} \end{equation} The augmented system is now solved for \(U\) using Gaussian elimination. In Matlab \[ U=A\backslash f \] The resulting solution \(U\) will have an extra element at the end, which is not used in the solution. This element, called \(\lambda \) was verified to be equal \(\frac{\mathbf{v\cdot f}}{\mathbf{v\cdot v}}.\)
Notice that last element of \(f\) was set to \(0\) in (7), this is because it corresponds to having selected an exact solution with a zero mean. When the last element of \(f\) is set to \(0\), this corresponds to having \({\displaystyle \sum \limits _{i=0}^{N+1}} U_{i}=0\), which implies the mean of the numerical solution is zero. This follows because \(null(A)\) was found to be \(\mathbf{1}^{T}\), meaning that, from looking at \(\begin{pmatrix} A & v\\ \mathbf{1}^{T} & 0 \end{pmatrix}\begin{pmatrix} U\\ \lambda \end{pmatrix} =\begin{pmatrix} f\\ 0 \end{pmatrix} ,\) the \(dot\left ( \mathbf{1}^{T},U\right ) =0\) , which is the same as saying that \({\displaystyle \sum \limits _{i=0}^{N+1}} U_{i}=0\) or the mean of the numerical solution is zero.
To summarize: A constant in the analytical solution was earlier selected which resulted in the mean of the exact solution to be zero. This constant was found to be \(\frac{-1}{6}\). The last entry in the augment \(f\) vector is set to be zero to cause the numerical solution to have zero mean as well.
Now that the numerical solution is found, the error vector \(e=U-u\) is found and its norms are calculated to complete the refinement study. Below is the result for both schemes described above.
Refinement study results
Conclusion
Two different schemes were used for solving part(b). Both schemes had an LTE of \(O\left ( h^{2}\right ) \), therefore they both were consistent numerical schemes. The result shows that \(\left \Vert e\right \Vert \) had \(O\left ( h^{2}\right ) \) implying \(\left \Vert A^{-1}\right \Vert \) had order \(O\left ( 1\right ) \) and that both scheme were stable.
It is observed that the second scheme required more iterations to coverage to the exact solution compared to the first scheme. Attempt was made to find if any error was made in the derivation of the scheme or in the code, but none found. Therefore, this showed that when comparing 2 schemes with the same theoretical order of accuracy, this does not necessarily imply they will have the same exact performance on the same numerical example. Therefore, before selecting a scheme, it would be useful to always run a number of numerical tests to compare schemes against each others on different numerical test problems.
Propose a discretization scheme for \(u_{xx}=f\) with \(u_{x}\left ( 0\right ) -\alpha u\left ( 0\right ) =g\) and \(u\left ( 1\right ) =b\)
Answer:
This problem has Robin boundary conditions on the left side and Dirichlet boundary conditions on the right side. The condition for existence of a solution implies\[ u_{x}\left ( 1\right ) -u_{x}\left ( 0\right ) ={\displaystyle \int \limits _{0}^{1}} f\left ( x\right ) dx \] Substituting the boundary conditions into the above results in\[ -\alpha u\left ( 0\right ) -g={\displaystyle \int \limits _{0}^{1}} f\left ( x\right ) dx \] Hence any \(f\left ( x\right ) \) which satisfies the above will make the problem a well posed one. In the following, 3 different schemes are presented using different methods of approximating the Robin boundary conditions. The last scheme resulted in an \(A\) matrix with the most desirable properties being tridiagonal and symmetric.
\(u_{x}\left ( 0\right ) \) is approximated by 2 points centered difference by introducing an imaginary point \(U_{-1}\) resulting in\[ u_{x}\left ( 0\right ) =\frac{U_{-1}-U_{1}}{2h}\] Substituting the above in the Robin boundary condition gives\begin{equation} \frac{U_{-1}-U_{1}}{2h}-\alpha U_{0}=g \tag{1} \end{equation} \(u_{xx}\left ( 0\right ) \) is approximated using the standard 3 point formula\begin{equation} \frac{U_{-1}-2U_{0}+U_{1}}{h^{2}}=f_{0} \tag{2} \end{equation} \(U_{-1}\) is now eliminated From (1) and (2). From equation (1)\[ U_{-1}=2h\left ( g+\alpha U_{0}\right ) +U_{1}\] Substituting into (2)\[ \frac{\left [ 2h\left ( g+\alpha U_{0}\right ) +U_{1}\right ] +2U_{0}+U_{1}}{h^{2}}=f_{0}\] And simplifying\[ \frac{2U_{0}\left ( 1+\alpha h\right ) +2U_{1}}{h^{2}}=f_{0}-\frac{2g}{h}\] This will be the equation to use at node \(j=0\) which is the first row in the \(A\) matrix.
For node \(j=N\), \(u_{xx}\left ( N\right ) \) is approximated using the standard 3 point formula\[ \frac{U_{N-1}-2U_{N}+U_{N+1}}{h^{2}}=f_{N}\] Since \(U_{N+1}=b\), the above becomes\[ \frac{U_{N-1}+2U_{N}}{h^{2}}=f_{N}-\frac{b}{h^{2}}\] Which is the equation for node \(j=N\). For all the remaining internal grid point, the standard 3 point formula is used to approximate \(u_{xx}\). Therefore, the system which now contains \(N+1\) equations is completed\[ \frac{1}{h^{2}}\begin{pmatrix} 2\left ( 1+\alpha h\right ) & 2 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & \cdot & \cdot & \cdot & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix}\begin{pmatrix} U_{0}\\ U_{1}\\ U_{2}\\ \cdot \\ \cdot \\ \cdot \\ U_{N-1}\\ U_{N}\end{pmatrix} =\begin{pmatrix} f_{0}-\frac{2g}{h}\\ f_{1}\\ f_{2}\\ \cdot \\ \cdot \\ \cdot \\ f_{N-1}\\ f_{N}-\frac{b}{h^{2}}\end{pmatrix} \]
Approximating \(u_{x}\left ( 0\right ) \) by 3 points forward difference\[ u_{x}\left ( 0\right ) =\frac{1}{h}\left ( \frac{3}{2}U_{0}-2U_{1}+\frac{1}{2}U_{2}\right ) \] Substituting the above in the Robin boundary condition gives\begin{align*} \frac{1}{h}\left ( \frac{3}{2}U_{0}-2U_{1}+\frac{1}{2}U_{2}\right ) -\alpha U_{0} & =g\\ U_{0}\left ( \frac{3}{2}-\alpha h\right ) -2U_{1}+\frac{1}{2}U_{2} & =hg \end{align*}
Dividing both sides by \(h^{2}\) to allow writing the matrix \(A\) with a common \(\frac{1}{h^{2}}\) factored out\[ \frac{U_{0}\left ( \frac{3}{2}-\alpha h\right ) -2U_{1}+\frac{1}{2}U_{2}}{h^{2}}=\frac{g}{h}\] For the right side, the same scheme is used as in the first scheme above\[ \frac{U_{N-1}+2U_{N}}{h^{2}}=f_{N}-\frac{b}{h^{2}}\] And the system becomes\[ \frac{1}{h^{2}}\begin{pmatrix} \left ( \frac{3}{2}-\alpha h\right ) & -2 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & \cdot & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \cdot & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & \cdot & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix}\begin{pmatrix} U_{0}\\ U_{1}\\ U_{2}\\ \cdot \\ \cdot \\ \cdot \\ U_{N-1}\\ U_{N}\end{pmatrix} =\begin{pmatrix} \frac{g}{h}\\ f_{1}\\ f_{2}\\ \cdot \\ \cdot \\ \cdot \\ f_{N-1}\\ f_{N}-\frac{b}{h^{2}}\end{pmatrix} \]
Approximating \(u_{x}\left ( 0\right ) \) by 2 points forward difference which has an LTE of \(O\left ( h\right ) \)\[ u_{x}\left ( 0\right ) =\frac{1}{h}\left ( U_{1}-U_{0}\right ) \] Substituting the above in the Robin boundary condition gives\begin{align*} \frac{1}{h}\left ( U_{1}-U_{0}\right ) -\alpha U_{0} & =g\\ U_{1}-U_{0}-\alpha hU_{0} & =hg\\ U_{0}\left ( -1-\alpha h\right ) +U_{1} & =hg \end{align*}
Dividing both sides by \(h^{2}\) to allow writing the matrix \(A\) with a common \(\frac{1}{h^{2}}\) factored out\[ \frac{\left [ U_{0}\left ( -1-\alpha h\right ) +U_{1}\right ] }{h^{2}}=\frac{g}{h}\] For the right side, the same scheme was used as in the first scheme above, resulting in\[ \frac{1}{h^{2}}\begin{pmatrix} \left ( -1-\alpha h\right ) & 1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & -2 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & \cdot & \cdot & \cdot & 0\\ 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1\\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{pmatrix}\begin{pmatrix} U_{0}\\ U_{1}\\ U_{2}\\ \cdot \\ \cdot \\ \cdot \\ U_{N-1}\\ U_{N}\end{pmatrix} =\begin{pmatrix} \frac{g}{h}\\ f_{1}\\ f_{2}\\ \cdot \\ \cdot \\ \cdot \\ f_{N-1}\\ f_{N}-\frac{b}{h^{2}}\end{pmatrix} \] This scheme has an LTE of \(O\left ( h\right ) \) since an \(O\left ( h\right ) \) approximation was used for \(u_{x}\left ( 0\right ) \)
Each of the above 3 schemes is a consistent scheme because for the first 2 schemes the LTE is \(O\left ( h^{2}\right ) \) and for the third scheme the LTE is \(O\left ( h\right ) \). A consistent scheme is one in which \(\left \Vert \tau \right \Vert \rightarrow 0\) as \(h\rightarrow 0\). It expected that the final error norm \(\left \Vert e\right \Vert \) will behave as \(O\left ( h^{2}\right ) \) when using the first 2 scheme and as \(O\left ( h\right ) \) when using the third scheme. This is provided \(A\) is stable. Meaning that \(\left \Vert A^{-1}\right \Vert \) is of order \(O\left ( 1\right ) \).
Considering the third scheme only since it is symmetric, and using the grid 2-norm, then given that \begin{align*} \left \Vert e\right \Vert _{2} & =\left \Vert A^{-1}\tau \right \Vert _{2}\\ & \leq \left \Vert A^{-1}\right \Vert _{2}\left \Vert \tau \right \Vert _{2} \end{align*}
And since for a symmetric matrix \(\left \Vert A^{-1}\right \Vert _{2}=\rho \left ( A^{-1}\right ) =\frac{1}{\left \vert \lambda _{\min }\right \vert }\) where \(\left \vert \lambda _{\min }\right \vert \) is the smallest eigenvalue of \(A\) in absolute value, then showing that \(\left \vert \lambda _{\min }\right \vert \geq 1\) would imply that \(\left \Vert e\right \Vert _{2}\) \(\leq \left \Vert \tau \right \Vert _{2}.\)
There two analytical methods for finding the \(\lambda _{\min }\) of \(A\). By solving the eigenvalue problem \(L\phi =\lambda \phi \) for the given boundary conditions or by using matrix algebra to solve for the eigenvalue of \(A\) directly. Attempt are made at both of these methods, but more time is needed to complete an analytical proof.
Therefore, to answer the question for this problem, a numerical method was used instead, where a refinement study was done using the third scheme found above in an attempt to show numerically that \(\left \Vert e\right \Vert \) is of the same order as \(\left \Vert \tau \right \Vert \) as expected. A well formed problem is constructed and used for the purpose of the numerical analysis part. The results and the conclusion follows.
The following numerical problem was select. \(u_{xx}=f\) with \(u_{x}\left ( 0\right ) -\alpha u\left ( 0\right ) =g\) and \(u\left ( 1\right ) =b\) with the following parameter values\begin{align*} \alpha & =1\\ g & =1\\ b & =1\\ f & =\cos \left ( x\right ) \end{align*}
This problem is well posed, and has the following analytical solution \[ u\left ( x\right ) =\frac{1}{2}(1+b-g-x+b\ x+g\ x+\cos (1)+x\ \cos (1)-2\cos (x)) \] And it satisfies the necessary condition \(-\alpha u\left ( 0\right ) -g={\displaystyle \int \limits _{0}^{1}} f\left ( x\right ) dx\)
The following are the results of the refinement study. The source code is in the appendix.
A scheme was implemented for Robin boundary conditions. LTE used was \(O\left ( h\right ) \) since this specific scheme generated an \(A\) matrix with desirable form (symmetrical and tridiagonal) . The expectation of total error using this scheme was to be \(O\left ( h\right ) \). This expectation was confirmed by doing a refinement study on the selected test problem. This result showed that the scheme is stable and therefore we conclude that \(\left \Vert A\right \Vert \) is of order \(O\left ( 1\right ) \) based on this test case.
From definition \begin{equation} \left \Vert \mathbf{e}\right \Vert _{\infty }=-\left \Vert B{\tau }\right \Vert _{\infty }\tag{1} \end{equation} Where \(B=A^{-1}\) was generated by the use of Green function as described on page 27 of the text book. Writing the product of the matrix \(B\) with the vector \(\mathbf{\tau }\) as \[{B\mathbf{\tau }}=b_{1}\tau _{1}+b_{1}\tau _{1}+\cdots +b_{N}\ \tau _{N}\] where \(\mathbf{b}_{j}\) is the \(j^{th}\) column of matrix \(B\), equation (1) becomes (all norms are to be taken as the max-norm, hence for clarify, the \(\infty \) is dropped in what follows) \begin{align} \left \Vert e\right \Vert & =-\left \Vert b_{1}{\tau }_{1}\mathbf{+}b_{1}\tau _{1}+\cdots +b_{N}{\tau }_{N}\right \Vert \nonumber \\ & \leq -\left ( \left \Vert b_{1}{\tau }_{1}\right \Vert +\left \Vert b_{2}{\tau }_{2}\right \Vert +\cdots +\left \Vert b_{N}{\tau }_{N}\right \Vert \right ) \nonumber \\ & \leq -\left ( \left \Vert b_{1}\right \Vert \left \vert{\tau }_{1}\right \vert +\left \Vert b_{2}\right \Vert \left \vert{\tau }_{2}\right \vert +\cdots +\left \Vert b_{N}\right \Vert \left \vert{\tau }_{N}\right \vert \right ) \tag{3} \end{align}
\(\left \Vert \mathbf{b}_{j}\right \Vert \) is the maximum element in the \(j^{th}\) column of the matrix \(B\). These elements are located along the diagonal of \(B\). Hence \(\left \Vert \mathbf{b}_{j}\right \Vert =\left \vert B_{jj}\right \vert \) and (3) becomes\begin{equation} \left \Vert \mathbf{e}\right \Vert \leq -\left ( \left \vert B_{11}\right \vert \ \left \vert{\tau }_{1}\right \vert +\left \vert B_{22}\right \vert \ \left \vert{\tau }_{2}\right \vert +\cdots +\left \vert B_{NN}\right \vert \ \left \vert{\tau }_{N}\right \vert \right ) \tag{4} \end{equation} From equation 2.46 in the textbook, \begin{equation} B_{ij}=\left \{ \begin{array} [c]{ccc}h(x_{j}-1)x_{i} & & i=1,2,\cdots j\\ h(x_{i}-1)x_{j} & & i=j,j+1,\cdots ,N \end{array} \right . \tag{5} \end{equation} To answer the question, assume that \(\left \vert{\tau }_{1}\right \vert =O\left ( h^{p}\right ) \), and assume for the moment, for generality, that the LTE at all the other points was different, hence for other points let \(\left \vert{\tau }_{j}\right \vert =O\left ( h^{q}\right ) \). Equation (4) becomes\begin{equation} \left \Vert \mathbf{e}\right \Vert \leq -O\left ( h^{p}\right ) \left \vert B_{11}\right \vert \ -O\left ( h^{q}\right ){\displaystyle \sum \limits _{k=2}^{N}} \left \vert B_{kk}\right \vert \ \ \tag{6} \end{equation} From (5), it is found that at point \(j=1\), \(\left \vert B_{11}\right \vert =h\left ( h-1\right ) h=h^{3}-h^{2}=O\left ( h^{2}\right ) \), while at the internal points, that is, at points near the middle, \(B_{jj}\) is approximated by \(O\left ( h\right ) \), hence (6) becomes\[ \left \Vert \mathbf{e}\right \Vert \leq -O\left ( h^{p}\right ) O\left ( h^{2}\right ) \ -O\left ( h^{q}\right ){\displaystyle \sum \limits _{k=2}^{N}} O\left ( h\right ) \ \ \] But \({\displaystyle \sum \limits _{k=2}^{N}} O\left ( h\right ) \sim \left ( N-1\right ) \times O\left ( h\right ) =\) \(O\left ( 1\right ) \) since \(N=\frac{1}{h}\), hence the above becomes\begin{align} \left \Vert \mathbf{e}\right \Vert & \leq -O\left ( h^{p}\right ) O\left ( h^{2}\right ) \ -O\left ( h^{q}\right ) \nonumber \\ & =-O\left ( h^{p+2}\right ) -O\left ( h^{q}\right ) \tag{7} \end{align}
Therefore, if \(q={p}\), i.e. if the LTE was the same at each grid point, including the first, the above simplifies to\begin{align*} \left \Vert \mathbf{e}\right \Vert & \leq -O\left ( h^{p+2}\right ) -O\left ( h^{p}\right ) \\ & \sim O\left ( h^{p}\right ) \end{align*}
Hence, having the LTE at the edge grid point be \(O\left ( h^{p}\right ) \) resulted in \(\left \Vert \mathbf{e}\right \Vert _{\infty }\) having an order of accuracy \(O\left ( h^{p}\right ) \) which is the expected. The smallest value of \(p\) that can be used for the edge point while still giving overall \(\left \Vert \mathbf{e}\right \Vert =O\left ( h^{2}\right ) \) can be found from (7)\[ \left \Vert \mathbf{e}\right \Vert =O\left ( h^{2}\right ) =-O\left ( h^{p+2}\right ) -O\left ( h^{q}\right ) \] Setting \(p=0\) and \(q=2\) results in\begin{align*} \left \Vert \mathbf{e}\right \Vert & =-O\left ( h^{2}\right ) -O\left ( h^{2}\right ) \\ & \sim O\left ( h^{2}\right ) \end{align*}
Hence \(p=0\) or \(O\left ( 1\right ) \), is possible for the LTE at the edge while arriving at an overall \(\left \Vert \mathbf{e}\right \Vert \sim O\left ( h^{2}\right ) \).
From equation 2.46 in the textbook, \begin{equation} B_{ij}=\left \{ \begin{array} [c]{ccc}h(x_{j}-1)x_{i} & & i=1,2,\cdots j\\ h(x_{i}-1)x_{j} & & i=j,j+1,\cdots ,N \end{array} \right . \tag{5} \end{equation} To answer the question, assume \(\left \vert{\tau }_{j=N/2}\right \vert =O\left ( h^{p}\right ) \), where \(j\) is the middle point, i.e. at \(x_{j}=\frac{1}{2}\)and assume for other points \(\left \vert{\tau }_{j}\right \vert =O\left ( h^{q}\right ) \), hence starting from \[ \left \Vert \mathbf{e}\right \Vert \leq -\left ( \left \vert B_{11}\right \vert \ \left \vert{\tau }_{1}\right \vert +\left \vert B_{22}\right \vert \ \left \vert{\tau }_{2}\right \vert +\cdots +\left \vert B_{NN}\right \vert \ \left \vert{\tau }_{N}\right \vert \right ) \] which was derived in part(a), then the above equation becomes\begin{equation} \left \Vert \mathbf{e}\right \Vert \leq -O\left ( h^{q}\right ){\displaystyle \sum \limits _{\substack{k=1\\k\neq N/2}}^{N}} \left \vert B_{kk}\right \vert -O\left ( h^{p}\right ) \left \vert B_{N/2,N/2}\right \vert \ \ \ \tag{8} \end{equation} At the middle point from (5), \(\left \vert B_{N/2,N/2}\right \vert \) can be found to be \(h\left ( \frac{1}{2}-1\right ) \frac{1}{2}=\frac{-1}{4}h\sim O\left ( h\right ) \), while for the points away from the middle, \(\left \vert B_{jj}\right \vert \) is approximated by \(O\left ( h^{2}\right ) \), hence (8) becomes\[ \left \Vert \mathbf{e}\right \Vert \leq -O\left ( h^{q}\right ) \left [ \left ( N-1\right ) \times O\left ( h^{2}\right ) \right ] \ -O\left ( h^{p}\right ) O\left ( h\right ) \ \ \] But \(\left ( N-1\right ) \times O\left ( h^{2}\right ) \) is \(O\left ( h\right ) \) since \(N=\frac{1}{h}\), hence the above becomes\begin{align} \left \Vert \mathbf{e}\right \Vert & \leq -O\left ( h^{q}\right ) O\left ( h\right ) \ -O\left ( h^{p+1}\right ) \nonumber \\ & =-O\left ( h^{q+1}\right ) -O\left ( h^{p+1}\right ) \tag{9} \end{align}
Therefore, if \(q={p}\), i.e. if the LTE was the same at each grid point, including the middle, then (9) becomes\begin{align*} \left \Vert \mathbf{e}\right \Vert & =-O\left ( h^{p+1}\right ) -O\left ( h^{p+1}\right ) \\ & =O\left ( h^{p+1}\right ) \end{align*}
Hence, having the LTE at the middle grid point be \(O\left ( h^{p}\right ) \) resulted in \(\left \Vert \mathbf{e}\right \Vert _{\infty }\) having an order of accuracy \(O\left ( h^{p+1}\right ) \) which is higher than the expected \(O\left ( h^{p}\right ) \). The smallest value of \(p\) that can be used while still giving \(\left \Vert \mathbf{e}\right \Vert =O\left ( h^{2}\right ) \) can be found from (9) \[ \left \Vert \mathbf{e}\right \Vert =O\left ( h^{2}\right ) =-O\left ( h^{q+1}\right ) -O\left ( h^{p+1}\right ) \] Setting \(p=1\) and \(q=1\) results in \(\left \Vert \mathbf{e}\right \Vert \sim O\left ( h^{2}\right ) \).
Hence \(p=1\) or \(O\left ( h\right ) \) is possible for the LTE at the middle point while still resulting in an overall \(\left \Vert \mathbf{e}\right \Vert =O\left ( h^{2}\right ) \)
To perform the numerical tests to verify the above conclusion, a numerical problem is used with a scheme which had an LTE of \(O\left ( h^{2}\right ) \). The LTE at the edge point only was then modified to test part(a) and the LTE at the middle point only was modified to test part(b).
Refinement study was done to verify that the order of accuracy remained \(O\left ( h^{2}\right ) \) after each modification.
For the numerical test, problem 1, part (a) was used for this test. That problem had an LTE of \(O\left ( h^{2}\right ) \) at each grid point, and the standard 3-point central difference was used for approximating \(u_{xx}.\)
To modify the LTE at a specific grid point, given that the LTE is defined as\[ \tau _{j}=\frac{u(j-1)-2u\left ( j\right ) +u\left ( j+1\right ) }{h^{2}}-f_{j}\] Where \(u\) above is the exact solution sampled at the given grid points, then to force the LTE to be \(O\left ( h\right ) \) in the middle, the above equation was modified for \(x_{j}=\frac{1}{2}\)(middle row of the system) to become\[ \tau _{j}=\frac{u(j-1)-2u\left ( j\right ) +u\left ( j+1\right ) }{h^{2}}-f_{j}+h \] And to force the LTE to be \(O\left ( 1\right ) \) at the first grid point, the LTE is modified for that point only (i.e. \(x_{j}=h\), the first row of the system) to become\[ \tau _{j}=\frac{u(j-1)-2u\left ( j\right ) +u\left ( j+1\right ) }{h^{2}}-f_{j}+1 \] The refinement study which was done for problem 1, part a, is now repeated with the above modifications. The following are the results generated.
Part (A) refinement
Before making the LTE modification for the first point, the error table generated was
One can see that \(\left \Vert e\right \Vert =O\left ( h^{2}\right ) \) as the ratio is \(4\). Now the \(f\) vector was modified as described above, and here is the small code segment how this was done. Complete code listing is at the end.
And now the program was run again
Examining the second table above shows that \(\left \Vert e\right \Vert _{\infty }=O\left ( h^{2}\right ) ,\) since the ratio of the error norm still converged to \(4\), implying \(p=2\) for \(\left \Vert e\right \Vert \). To compare the loglog plots, here are plots before and after the modification is made to the LTE.
Part(b) refinement study
Before making the LTE modification for the middle point, the error table generated was
One can see that \(\left \Vert e\right \Vert =O\left ( h^{2}\right ) \) because the ratio is \(4.\) Now the \(f\) vector was modified as described above, to force the middle point to have an LTE of \(O\left ( h\right ) \), here is the small code segment showing the modification.
And now the program was run again
Examining the second table above, shows that \(\left \Vert e\right \Vert _{\infty }=O\left ( h^{2}\right ) ,\) since the ratio of the error norm still converged to \(4\), implying \(p=2\) for \(\left \Vert e\right \Vert \). To compare the loglog plots, here are plots before and after the modification is made to the LTE.
This concludes the refinement study verifying results from part(a) and part(c)