Numerical Methods Qualification Exam Problems and Solutions (University of Maryland)/Jan08 667

Problem 4

Consider the system $x=1+h\frac{e^{-x^2}}{1+y^2}\text{,}y=.5+h\arctan (x^2+y^2)$.

Problem 4a

Show that if the parameter $h>0$ is chosen sufficiently small, then this system has a unique solution $(x^{*},y^{*})$ within some rectangular region.

Solution 4a

The system of equations may be expressed in matrix notation.

$f(x,y)=\left[\begin{array}{c}1+h\frac{e^{-x^2}}{1+y^2}\\ \\ .5+h\arctan (x^2+y^2)\end{array}\right]$

The Jacobian of $f(x,y)$, $J(f)$, is computed using partial derivatives

$J(f)=\left[\begin{array}{cc}\frac{h}{1+y^2}{e}^{-x^2}(-2x)& \frac{-h2y{e}^{-x^2}}{(1+y^2{)}^2}\\ h\cdot \frac{2x}{1+(x^2+y^2{)}^2}& h\cdot \frac{2y}{1+(x^2+y^2{)}^2}\end{array}\right]$

If $h$ is sufficiently small and $x,y$ are restricted to a bounded region $B$,

$\underset{C}{\underbrace{\Vert J(f)\Vert }}<1$

From the mean value theorem, for $\overrightarrow{x_1},\overrightarrow{x_2}\in B$, there exists $\overrightarrow{\xi }$ such that

$f(\overrightarrow{x_1})-f(\overrightarrow{x_2})=J(f(\overrightarrow{\xi }))(\overrightarrow{x_1}-\overrightarrow{x_2})$

Since $\Vert J(f)\Vert$ is bounded in the region $B$ give sufficiently small $h$

$\Vert f(\overrightarrow{x_1})-f(\overrightarrow{x_2})\Vert \le C\Vert \overrightarrow{x_1}-\overrightarrow{x_2}\Vert$

Therefore, $f$ is a contraction and from the contraction mapping theorem there exists a unique fixed point (solution) in a rectangular region.

Problem 4b

Derive a fixed point iteration scheme for solving the system and show that it converges. $$

Solution 4b

Use Newton Method

To solve this problem, we can use the Newton Method. In fact, we want to find the zeros of

$g(x,y)=\left[\begin{array}{c}x\\ \\ y\end{array}\right]-\left[\begin{array}{c}1+h\frac{e^{-x^2}}{1+y^2}\\ \\ .5+h\arctan (x^2+y^2)\end{array}\right]$

The Jacobian of $g(x,y)$, $J(g)$, is computed using partial derivatives

$J(g)=\left[\begin{array}{cc}1-\frac{h}{1+y^2}{e}^{-x^2}(-2x)& \frac{-h{e}^{-x^2}}{1+y^2}\\ h\cdot \frac{2x}{1+(x^2+y^2{)}^2}& 1-h\cdot \frac{2y}{1+(x^2+y^2{)}^2}\end{array}\right]$

Then, the Newton method for solving this linear system of equations is given by

$\left[\begin{array}{c}{x}_{k+1}\\ \\ y_{k+1}\end{array}\right]=\left[\begin{array}{c}{x}_k\\ \\ y_k\end{array}\right]-\underset{J(g{)}^{-1}}{\underbrace{\frac{1}{det(J(g))}\left[\begin{array}{cc}1-h\cdot \frac{2y}{1+(x^2+y^2{)}^2}& \frac{h{e}^{-x^2}}{1+y^2}\\ -h\cdot \frac{2x}{1+(x^2+y^2{)}^2}& 1-\frac{h}{1+y^2}{e}^{-x^2}(-2x)\end{array}\right]}}g(x_k,y_k)$

Show convergence by showing Newton hypothesis hold

Jacobian of g is Lipschitz

We want to show that $J(g)$ is a Lipschitz function. In fact,

$\begin{array}{cc}\Vert J(g)(\overrightarrow{x_1})-J(g)(\overrightarrow{x_2})\Vert & =\Vert (I-J(f))(\overrightarrow{x_1})-(I-J(f))(\overrightarrow{x_2})\Vert \\ & =\Vert \overrightarrow{x_1}-\overrightarrow{x_2}-J(f)\overrightarrow{x_1}+J(f)\overrightarrow{x_2}\Vert \\ & \le \Vert \overrightarrow{x_1}-\overrightarrow{x_2}\Vert +\Vert J(f)\overrightarrow{x_2}-J(f)\overrightarrow{x_1}\Vert \end{array}$

and now, using that $J(f)$ is Lipschitz, we have

$\Vert J(g)(\overrightarrow{x_1})-J(g)(\overrightarrow{x_2})\Vert \le (1+C)\Vert (\overrightarrow{x_1})-(\overrightarrow{x_2})\Vert$

Jacobian of g is invertible

Since $J(f)$ is a contraction, the spectral radius of the Jacobian of $f$ is less than 1 i.e.

$\rho (J(f))<1$.

On the other hand, we know that the eigenvalues of $J(g)$ are $\lambda (J(g))=1-\lambda (J(f))$.

Then, it follows that $det(I-(J(g)))\ne 0$ or equivalently $J(g)$ is invertible.

inverse of Jacobian of g is bounded

$\underset{J(g{)}^{-1}}{\underbrace{\frac{1}{det(J(g))}\left[\begin{array}{cc}1-h\cdot \frac{2y}{1+(x^2+y^2{)}^2}& \frac{h{e}^{-x^2}}{1+y^2}\\ -h\cdot \frac{2x}{1+(x^2+y^2{)}^2}& 1-\frac{h}{1+y^2}{e}^{-x^2}(-2x)\end{array}\right]}}$

Since $J(g{)}^{-1}$ exists, $\text{det}(J(g))\ne 0$. Given a bounded region (bounded $x,y$), each entry of the above matrix is bounded. Therefore the norm is bounded.

norm of (Jacobian of g)^-1 (x_0) * g(x_0) bounded

$\Vert J(g{)}^{-1}(x_0)*g(x_0)\Vert <\mathrm{\infty }$ since $J(g{)}^{-1}$ and $g(x)$ is bounded.

Then, using a good enough approximation $(x_0,y_0)$, we have that the Newton's method is at least quadratically convergent, i.e,

$\Vert (x_{k-1},y_{k+1})-(x^{*},y^{*})\Vert \le K\Vert (x_k,y_k)-(x^{*},y^{*}){\Vert }^2.$

Problem 5a

Outline the derivation of the Adams-Bashforth methods for the numerical solution of the initial value problem $y^{\prime }=f(x,y)\text{, }y(x_0)=y_0$.

Solution 5a

We want to solve the following initial value problem: $y^{\prime }=f(x,y)\text{, }y(x_0)=y_0$.

First, we integrate this expression over $[t_n,t_{n+1}]$, to obtain

$y(t_{n+1})=y(t_n)+\underset{t_n}{\overset{t_{n+1}}{\int }}f(t,y(t))dt$,

To approximate the integral on the right hand side, we approximate its integrand $f(t,y(t))$ using an appropriate interpolation polynomial of degree $p$ at $t_n,t_{n-1},...,t_{n-p}$.

This idea generates the Adams-Bashforth methods.

$y_{n+1}=y_n+\underset{t_n}{\overset{t_{n-1}}{\int }}\sum _{k=0}^{p}g(t_{n-k})l_k(t)dt$,

where, $y_n$ denotes the approximated solution, $g(t)\equiv f(t,y)$ and $l_k(t)$ denotes the associated Lagrange polynomial.

Problem 5b

Derive the Adams-Bashforth formula $y_{i+1}=y_i+h[-\frac{1}{2}f(x_{i-1},y_{i-1})+\frac{3}{2}f(x_i,y_i)](1)$

Solution 5b

From (a) we have

$y_{i+1}=y_i+\underset{x_i}{\overset{x_{i+1}}{\int }}fdx$ where $\underset{x_i}{\overset{x_{i+1}}{\int }}fdx\approx \underset{x_i}{\overset{x_{i+1}}{\int }}[\frac{x-x_i}{x_{i-1}-x_i}{f}_{i-1}+\frac{x-x_{i-1}}{x_i-x_{i-1}}{f}_i]dx$

Then if we let $x=x_i+sh$, where h is the fixed step size we get

$\begin{array}{cc}dx& =hds\\ x-x_i& =sh\\ x-x_{i-1}& =(1+s)h\\ x_i-x_{i-1}& =h\end{array}$

$h\underset{0}{\overset{1}{\int }}[\frac{sh}{-h}{f}_{i-1}+\frac{(1+s)h}{h}{f}_i]ds=h[-\frac{1}{2}{f}_{i-1}+\frac{3}{2}{f}_i]$

So we have the Adams-Bashforth method as desired

$y_{i+1}=y_i+h[-\frac{1}{2}f(x_{i-1},y_{i-1})+\frac{3}{2}f(x_i,y_i)]$

Problem 5c

Analyze the method (1). To be specific, find the local truncation error, prove convergence and find the order of convergence.

Solution 5c

Find Local Truncation Error Using Taylor Expansion

Note that $y_i\approx y(t_i)$. Also, denote the uniform step size $\mathrm{\Delta }t$ as h. Hence,

$t_{i-1}=t_i-h$

$t_{i+1}=t_i+h$

Therefore, the given equation may be written as

$y(t_i+h)-y(t_i)\approx h[-\frac{1}{2}{y}^{\prime }(t_i-h)+\frac{3}{2}{y}^{\prime }(t_i)]$

Expand Left Hand Side

Expanding about $t_i$, we get

$\begin{array}{cccc}\text{Order}& y(t_i+h)& y(t_i)& \mathrm{\Sigma }\\ & & & \\ 0& y(t_i)& y(t_i)& 0\\ & & & \\ 1& -y^{\prime }(t_i)h& 0& -y^{\prime }(t_i)h\\ & & & \\ 2& \frac{1}{2!}{y}^{″}(t_i)h^2& 0& \frac{1}{2}{y}^{″}(t_i)h^2\\ & & & \\ 3& -\frac{1}{3!}{y}^{‴}(t_i)h^3& 0& -\frac{1}{6}{y}^{‴}(t_i)h^3\\ & & & \\ 4& 𝒪(h^4)& 0& 𝒪(h^4)\\ & & & \end{array}$

Expand Right Hand side

Also expanding about $t_i$ gives

$\begin{array}{ccccc}\text{Order}& \frac{13}{12}h{x}^{\prime }(t_i+2h)& -\frac{5}{3}h{x}^{\prime }(t_i+h)& -\frac{5}{12}h{x}^{\prime }(t_i)& \mathrm{\Sigma }\\ & & & & \\ 0& 0& 0& 0& 0\\ & & & & \\ 1& \frac{13}{12}h{x}^{\prime }(t_i)& -\frac{5}{3}h{x}^{\prime }(t_i)& -\frac{5}{12}h{x}^{\prime }(t_i)& -1\cdot h{x}^{\prime }(t_i)\\ & & & & \\ 2& \frac{13}{12}(2h)h{x}^{″}(t_i)& -\frac{5}{3}hh{x}^{″}(t_i)& 0& \frac{1}{2}{h}^2{x}^{″}(t_i)\\ & & & & \\ 3& \frac{13}{12}(2h{)}^2\frac{h{x}^{‴}(t_i)}{2}& -\frac{5}{3}{h}^2\frac{h{x}^{‴}(t_i)}{2}& 0& \frac{4}{3}{h}^3{x}^{‴}(t_i)\\ & & & & \\ 4& 𝒪(h^4)& 𝒪(h^4)& 0& 𝒪(h^4)\end{array}$

Show Convergence by Showing Stability and Consistency

A method is convergent if and only if it is both stable and consistent.

Stability

It is easy to show that the method is zero stable as it satisfies the root condition. So stable.

Consistency

Truncation error is of order 2. Truncation error tends to zero as h tends to zeros. So the method is consistent.

Order of Convergence

Dahlquist principle: consistency + stability = convergent. Order of convergence will be 2.

Problem 6

Consider the problem $-u^{″}+u=f(x),0\le x\le 1,u(0)=u(1)=0.(2)$

Problem 6a

Give a variational formulation of (2), i.e. express (2) as $u\in H,B(x,y)=F(v),\mathrm{\forall }v\in H(3)$ Define the Space H, the bilinear form B and the linear functional F, and state the relation between (2) and (3).

Solution 6a

Multiplying (2) by a test function and using integration by parts we obtain:

$-\underset{0}{\overset{1}{\int }}{u}^{″}vdx+\underset{0}{\overset{1}{\int }}uvdx=\underset{0}{\overset{1}{\int }}fvdx$

$\underset{0}{\underbrace{{-u^{\prime }v|}_0^1}}+\underset{0}{\overset{1}{\int }}{u}^{\prime }{v}^{\prime }dx+\underset{0}{\overset{1}{\int }}uvdx=\underset{0}{\overset{1}{\int }}fvdx$

Thus, the weak form or variational form associated with the problem (2) reads as follows: Find $u\in H$ such that

$B(u,v)=\underset{0}{\overset{1}{\int }}{u}^{\prime }{v}^{\prime }dx+\underset{0}{\overset{1}{\int }}uvdx=\underset{0}{\overset{1}{\int }}fvdx=F(v)$ for all $v\in H,$

where $H:=H_0^1(0,1)=\{v:v\in H^1,v(0)=v(1)=0\}$.

Problem 6b

Let $0=x_0<x_1<\cdots <x_n=1$ be a mesh on $[0,1]$ with $h=\underset{j=0,1,\dots ,n-1}{\max }(x_{j+1}-x_j)$, and let $V_h=\{u:u\text{ continuous on }[0,1],u{|}_{[x_j,x_{j+1}]}\text{ is linear for each }j,u(0)=u(1)=0\}$. Define the finite element approximation, $u_h\text{ to }u$ based on the approximation space $V_h$. What can be said about $\Vert u-u_h{\Vert }_1$ the error on the Sobolev norm on $H^1(0,1)$?

Solution 6b

Define piecewise linear basis

For our basis of $V_h$, we use the set of hat functions $\{{\varphi }_j{\}}_{j=1}^{n-1}$, i.e., for $j=1,2,\dots ,n-1$

${\varphi }_j(x)=\{\begin{array}{cc}\frac{x-x_{j-1}}{x_j-x_{j-1}}& \text{for }x\in [x_{j-1},x_j]\\ \frac{x_{j+1}-x}{x_{j+1}-x_j}& \text{for }x\in [x_j,x_{j+1}]\\ 0& \text{otherwise}\end{array}$

Define u_h

Since $\{{\varphi }_j{\}}_{j=1}^{n-1}$ is a basis for $V_h$, and $u_h\in V_h$ we have

$u_h=\sum _{j=1}^{n-1}{u}_{hj}{\varphi }_j$.

Discrete Problem

Now, we can write the discrete problem: Find $u_h\in V_h$ such that

$B(u_h,v_h)=F(v_h)$ for all $v_h\in V_h.$

If we consider that $\{{\varphi }_j{\}}_{j=1}^{n-1}$ is a basis of $V_h$ and the linearity of the bilinear form $B$ and the functional $F$, we obtain the equivalent problem:

Find $u_h=(u_{h,1},\dots ,u_{h,n-1})\in {ℝ}^{n-1}$ such that

$\sum _{j=1}^{n-1}{u}_{h,j}\{\underset{0}{\overset{1}{\int }}{{\varphi }_j}^{\prime }{{\varphi }_i}^{\prime }dx+\underset{0}{\overset{1}{\int }}{\varphi }_j{\varphi }_{i}dx\}=\underset{0}{\overset{1}{\int }}f{\varphi }_{i}dx,i=1,\dots ,n-1.$

This last problem can be formulated as a matrix problem as follows:

Find $u_h=(u_{h,1},\dots ,u_{h,n-1})\in {ℝ}^{n-1}$ such that

$A{u}_h=F,$

where $A_{ij}=\underset{0}{\overset{1}{\int }}{{\varphi }_j}^{\prime }{{\varphi }_i}^{\prime }dx+\underset{0}{\overset{1}{\int }}{\varphi }_j{\varphi }_{i}dx$ and $F_j=\underset{0}{\overset{1}{\int }}f{\varphi }_{i}dx$.

Bound error

In general terms, we can use Cea's Lemma to obtain

$\Vert u-u_h{\Vert }_1\le C\underset{v_h\in V_h}{\inf }\Vert u-v_h{\Vert }_1$

In particular, we can consider $v_h$ as the Lagrange interpolant, which we denote by $v_I$. Then,

$\Vert u-u_h{\Vert }_1\le C\Vert u-v_I{\Vert }_1\le Ch\Vert u{\Vert }_{H^2(0,1)}$.

It's easy to prove that the finite element solution is nodally exact. Then it coincides with the Lagrange interpolant, and we have the following punctual estimation:

$\Vert u(x)-u_h(x){\Vert }_{L^{\mathrm{\infty }}([x_{i-1},x_i])}\le \underset{x\in [x_{i-1},x_i]}{\max }\frac{u^{(2)}(x)}{(2)!}{\left(\frac{h}{2}\right)}^2$

Problem 6c

Derive the estimate for $\Vert u-u_h{\Vert }_0$, the error in $L_2(0,1)$. Hint: Let w solve (#) :$\{\begin{array}{c}-w^{″}+w=u-u_h,\\ w(0)=w(1)=0.\end{array}$ We characterize $w$ variationally as $w\in H_0^1,B(w,v)=\int (u-u_h)vdx,\mathrm{\forall }v\in H_0^1$. Let $v=u-u_h$ to get $B(w,u-u_h)=\int (u-u_h{)}^{2}dx=\Vert u-u_h{\Vert }_{L_2}^2.(4)$ Use formula (4) to estimate $\Vert u-u_h{\Vert }_{L_2}$.

Solution 6c

Continuity of Bilinear Form

$\begin{array}{cc}B(u,v)& =\int uv+\int u^{\prime }{v}^{\prime }\\ & \le \Vert u{\Vert }_0\Vert v{\Vert }_0+\Vert u^{\prime }{\Vert }_0\Vert v^{\prime }{\Vert }_0\text{ (Cauchy-Schwarz in L2) }\\ & \le (\Vert u{\Vert }_0^2+\Vert u^{\prime }{\Vert }_0^2{)}^{\frac{1}{2}}(\Vert v{\Vert }_0^2+\Vert v^{\prime }{\Vert }^2{)}^{\frac{1}{2}}\text{ ( Cauchy-Schwarz in R2 ) }\\ & =\Vert u{\Vert }_1\cdot \Vert v{\Vert }_1\end{array}$

Orthogonality of the Error

$B(u-u_h,v_h)=0,\mathrm{\forall }{v}_h\in V_h$.

Bound for L2 norm of w

$\begin{array}{cc}\Vert w{\Vert }_0^2& \le \Vert w{\Vert }_1^2\\ & =B(w,w)\\ & =\int (u-u_h)w\\ & \le \Vert u-u_h{\Vert }_0\Vert w{\Vert }_0\text{ (Cauchy-Schwarz in L2 ) }\end{array}$

Hence,

$(*)\Vert w{\Vert }_0\le \Vert u-u_h{\Vert }_0$

Bound for L2 norm of w

From $(\mathrm{\#})$, we have

$-w^{″}+w=u-u_h$

Then,

$\begin{array}{cc}\Vert w^{″}{\Vert }_0& \le \Vert w{\Vert }_0+\Vert u-u_h{\Vert }_0\text{ (triangle inequality) }\\ & \le \Vert u-u_h{\Vert }_0+\Vert u-u_h{\Vert }_0\text{ (by (*) ) }\\ & =2\Vert u-u_h{\Vert }_0\end{array}$

Bound L2 Error

$\begin{array}{cc}\Vert u-u_h{\Vert }_{L^2}^2& =B(w,u-u_h)\text{ ( from equation (4) ) }\\ & =B(w,u-u_h)-\underset{0}{\underbrace{B(w_h,u-u_h)}}\\ & =B(w-w_h,u-u_h)\\ & \le \Vert w-w_h{\Vert }_{H^1(0,1)}\Vert u-u_h{\Vert }_{H^1(0,1)}\text{ (Continuity of Bilinear Form) }\\ & \le Ch\Vert w{\Vert }_{H^2(0,1)}\Vert u-u_h{\Vert }_{H^1(0,1)}\text{ (Cea Lemma and Polynomial Interpolation Error) }\end{array}$

Finally, from (#), we have that $\Vert w{\Vert }_{H^2(0,1)}\le C\Vert u-u_h{\Vert }_{L_2(0,1)}$. Then,

$\Vert u-u_h{\Vert }_{L^2(0,1)}^2\le Ch\Vert u-u_h{\Vert }_{L_2(0,1)}\Vert u-u_h{\Vert }_{H^1(0,1)}$,

or equivalently,

$\Vert u-u_h{\Vert }_{L^2(0,1)}\le Ch\Vert u-u_h{\Vert }_{H^1(0,1)}\le C{h}^2\Vert u{\Vert }_{H^2(0,1)}$.

Problem 6d

Suppose $\{{\varphi }_1^h,\dots ,{\varphi }_{N_h}^h\}$ is a basis for $V_h\text{ where }{N}_h=\dim V_h\text{, so }{u}_h=\sum _{j=1}^{N_h}{c}_j^h{\varphi }_j^h\text{, for appropriate coefficients }{c}_j^h.$. Show that $\Vert u-u_h{\Vert }_1^2=\Vert u{\Vert }_1^2-C^{T}AC$ where $C=[c_1^h,\dots ,c_{N_h}^h{]}^T\text{ and }A$ is the stiffness matrix.

Solution 6d

We know that

$\begin{array}{cc}\Vert u-u_h{\Vert }_1^2& =⟨u-u_h,u-u_h{⟩}_{H^1(0,1)}\\ & =⟨u-u_h,u{⟩}_{H^1(0,1)}-\underset{0}{\underbrace{⟨u-u_h,u_h{⟩}_{H^1(0,1)}}}\\ & =\Vert u{\Vert }_{H^1(0,1)}^2-⟨u_h,u{⟩}_{H^1(0,1)}\\ & =\Vert u{\Vert }_{H^1(0,1)}^2-\Vert u_h{\Vert }_{H^1(0,1)}^2\end{array}$

where the substitution in the last lines come from the orthogonality of error.

Now, $\Vert u_h{\Vert }_{H^1(0,1)}^2=\sum _{j=1}^{n-1}\sum _{i=1}^{n-1}{u}_{h,j}{u}_{h,i}(\underset{0}{\overset{1}{\int }}{{\varphi }_i}^{\prime }{{\varphi }_j}^{\prime }dx+\underset{0}{\overset{1}{\int }}{\varphi }_i{\varphi }_{j}dx)=\sum _{j=1}^{n-1}\sum _{i=1}^{n-1}{u}_{h,j}{A}_{ij}{u}_{h,i}=C^{t}AC.$

Then, we have obtained

$\Vert u-u_h{\Vert }_1^2=\Vert u{\Vert }_{H^1(0,1)}^2-C^{t}AC.$