# 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}}}{\mbox{,}}\qquad 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)={\begin{bmatrix}1+h{\frac {e^{-x^{2}}}{1+y^{2}}}\\\\.5+h\arctan(x^{2}+y^{2})\end{bmatrix}}\!\,$
The Jacobian of $f(x,y)\!\,$, $J(f)\!\,$, is computed using partial derivatives
$J(f)={\begin{bmatrix}{\frac {h}{1+y^{2}}}e^{-x^{2}}(-2x)&{\frac {-h2ye^{-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{bmatrix}}\!\,$
If $h\!\,$ is sufficiently small and $x,y\!\,$ are restricted to a bounded region $B\!\,$,
$\underbrace {\|J(f)\|} _{C}<1\!\,$
From the mean value theorem, for ${\vec {x_{1}}},{\vec {x_{2}}}\in B\!\,$, there exists ${\vec {\xi }}\!\,$ such that
$f({\vec {x_{1}}})-f({\vec {x_{2}}})=J(f({\vec {\xi }}))({\vec {x_{1}}}-{\vec {x_{2}}})\!\,$
Since $\|J(f)\|\!\,$ is bounded in the region $B\!\,$ give sufficiently small $h\!\,$
$\|f({\vec {x_{1}}})-f({\vec {x_{2}}})\|\leq C\|{\vec {x_{1}}}-{\vec {x_{2}}}\|\!\,$
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)={\begin{bmatrix}x\\\\y\end{bmatrix}}-{\begin{bmatrix}1+h{\frac {e^{-x^{2}}}{1+y^{2}}}\\\\.5+h\arctan(x^{2}+y^{2})\end{bmatrix}}\!\,$
The Jacobian of $g(x,y)\!\,$, $J(g)\!\,$, is computed using partial derivatives
$J(g)={\begin{bmatrix}1-{\frac {h}{1+y^{2}}}e^{-x^{2}}(-2x)&{\frac {-he^{-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{bmatrix}}\!\,$
Then, the Newton method for solving this linear system of equations is given by
${\begin{bmatrix}x_{k+1}\\\\y_{k+1}\end{bmatrix}}={\begin{bmatrix}x_{k}\\\\y_{k}\end{bmatrix}}-\underbrace {{\frac {1}{det(J(g))}}{\begin{bmatrix}1-h\cdot {\frac {2y}{1+(x^{2}+y^{2})^{2}}}&{\frac {he^{-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{bmatrix}}} _{J(g)^{-1}}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{aligned}\|J(g)({\vec {x_{1}}})-J(g)({\vec {x_{2}}})\|&=\|(I-J(f))({\vec {x_{1}}})-(I-J(f))({\vec {x_{2}}})\|\\&=\|{\vec {x_{1}}}-{\vec {x_{2}}}-J(f){\vec {x_{1}}}+J(f){\vec {x_{2}}}\|\\&\leq \|{\vec {x_{1}}}-{\vec {x_{2}}}\|+\|J(f){\vec {x_{2}}}-J(f){\vec {x_{1}}}\|\end{aligned}}\!\,$
and now, using that $J(f)\!\,$ is Lipschitz, we have
$\|J(g)({\vec {x_{1}}})-J(g)({\vec {x_{2}}})\|\leq (1+C)\|({\vec {x_{1}}})-({\vec {x_{2}}})\|\!\,$

#### 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)))\neq 0\!\,$ or equivalently $J(g)\!\,$ is invertible.

#### inverse of Jacobian of g is bounded

$\underbrace {{\frac {1}{det(J(g))}}{\begin{bmatrix}1-h\cdot {\frac {2y}{1+(x^{2}+y^{2})^{2}}}&{\frac {he^{-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{bmatrix}}} _{J(g)^{-1}}\!\,$
Since $J(g)^{-1}\!\,$ exists, ${\mbox{det}}(J(g))\neq 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

$\|J(g)^{-1}(x_{0})*g(x_{0})\|<\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,
$\|(x_{k-1},y_{k+1})-(x^{*},y^{*})\|\leq K\|(x_{k},y_{k})-(x^{*},y^{*})\|^{2}.\!\,$

## Problem 5a

* Outline the derivation of the Adams-Bashforth methods for the numerical solution of the initial value problem $y'=f(x,y){\mbox{, }}y(x_{0})=y_{0}\!\,$.

## Solution 5a

We want to solve the following initial value problem: $y'=f(x,y){\mbox{, }}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})+\int _{t_{n}}^{t_{n+1}}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}+\int _{t_{n}}^{t_{n-1}}\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\left[-{\frac {1}{2}}f(x_{i-1},y_{i-1})+{\frac {3}{2}}f(x_{i},y_{i})\right]\qquad \qquad (1)\!\,$

## Solution 5b

From (a) we have
    $y_{i+1}=y_{i}+\int _{x_{i}}^{x_{i+1}}fdx\!\,$ where $\int _{x_{i}}^{x_{i+1}}fdx\approx \int _{x_{i}}^{x_{i+1}}\left[{\frac {x-x_{i}}{x_{i-1}-x_{i}}}f_{i-1}+{\frac {x-x_{i-1}}{x_{i}-x_{i-1}}}f_{i}\right]dx\!\,$
Then if we let $x=x_{i}+sh\!\,$, where h is the fixed step size we get
${\begin{aligned}dx&=hds\\x-x_{i}&=sh\\x-x_{i-1}&=(1+s)h\\x_{i}-x_{i-1}&=h\end{aligned}}$
    $h\int _{0}^{1}\left[{\frac {sh}{-h}}f_{i-1}+{\frac {(1+s)h}{h}}f_{i}\right]ds=h\left[-{\frac {1}{2}}f_{i-1}+{\frac {3}{2}}f_{i}\right]\!\,$
So we have the Adams-Bashforth method as desired
    $y_{i+1}=y_{i}+h\left[-{\frac {1}{2}}f(x_{i-1},y_{i-1})+{\frac {3}{2}}f(x_{i},y_{i})\right]\!\,$

## 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 $\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\left[-{\frac {1}{2}}y'(t_{i}-h)+{\frac {3}{2}}y'(t_{i})\right]\!\,$

#### Expand Left Hand Side

Expanding about $t_{i}\!\,$, we get
${\begin{array}{|c|c|c|c|}{\mbox{Order}}&y(t_{i}+h)&y(t_{i})&\Sigma \\\hline &&&\\0&y(t_{i})&y(t_{i})&0\\&&&\\1&-y'(t_{i})h&0&-y'(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&{\mathcal {O}}(h^{4})&0&{\mathcal {O}}(h^{4})\\&&&\\\hline \end{array}}$

#### Expand Right Hand side

Also expanding about $t_{i}\!\,$ gives
${\begin{array}{|c|c|c|c|c|}{\mbox{Order}}&{\frac {13}{12}}hx'(t_{i}+2h)&-{\frac {5}{3}}hx'(t_{i}+h)&-{\frac {5}{12}}hx'(t_{i})&\Sigma \\\hline &&&&\\0&0&0&0&0\\&&&&\\1&{\frac {13}{12}}hx'(t_{i})&-{\frac {5}{3}}hx'(t_{i})&-{\frac {5}{12}}hx'(t_{i})&-1\cdot hx'(t_{i})\\&&&&\\2&{\frac {13}{12}}(2h)hx''(t_{i})&-{\frac {5}{3}}hhx''(t_{i})&0&{\frac {1}{2}}h^{2}x''(t_{i})\\&&&&\\3&{\frac {13}{12}}(2h)^{2}{\frac {hx'''(t_{i})}{2}}&-{\frac {5}{3}}h^{2}{\frac {hx'''(t_{i})}{2}}&0&{\frac {4}{3}}h^{3}x'''(t_{i})\\&&&&\\4&{\mathcal {O}}(h^{4})&{\mathcal {O}}(h^{4})&0&{\mathcal {O}}(h^{4})\\\hline \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\leq x\leq 1,\;\;u(0)=u(1)=0.\qquad \qquad (2)$

## Problem 6a

* Give a variational formulation of (2), i.e. express (2) as
* $u\in H,\;\;B(x,y)=F(v),\;\;\forall v\in H\qquad \qquad (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:
$-\int _{0}^{1}u''vdx+\int _{0}^{1}uvdx=\int _{0}^{1}fvdx\!\,$
$\underbrace {\left.-u'v\right|_{0}^{1}} _{0}+\int _{0}^{1}u'v'dx+\int _{0}^{1}uvdx=\int _{0}^{1}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)=\int _{0}^{1}u'v'dx+\int _{0}^{1}uvdx=\int _{0}^{1}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=\max _{j=0,1,\dots ,n-1}(x_{j+1}-x_{j})\!\,$, and let
* $V_{h}=\left\{u:u{\mbox{ continuous on }}[0,1],u|_{[x_{j},x_{j+1}]}{\mbox{ is linear for each }}j,u(0)=u(1)=0\right\}\!\,$.
* Define the finite element approximation, $u_{h}{\mbox{ to }}u\!\,$ based on the approximation space $V_{h}\!\,$. What can be said about $\|u-u_{h}\|_{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 $\{\phi _{j}\}_{j=1}^{n-1}\!\,$, i.e., for $j=1,2,\ldots ,n-1\!\,$
$\phi _{j}(x)={\begin{cases}{\frac {x-x_{j-1}}{x_{j}-x_{j-1}}}&{\mbox{for }}x\in [x_{j-1},x_{j}]\\{\frac {x_{j+1}-x}{x_{j+1}-x_{j}}}&{\mbox{for }}x\in [x_{j},x_{j+1}]\\0&{\mbox{otherwise}}\end{cases}}\!\,$

### Define u_h

Since $\{\phi _{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}\phi _{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 $\{\phi _{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 \mathbb {R} ^{n-1}\!\,$ such that
$\sum _{j=1}^{n-1}u_{h,j}\left\{\int _{0}^{1}\phi _{j}'\phi _{i}'dx+\int _{0}^{1}\phi _{j}\phi _{i}dx\right\}=\int _{0}^{1}f\phi _{i}dx,\quad 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 \mathbb {R} ^{n-1}\!\,$ such that
$Au_{h}=F,\!\,$
where $A_{ij}=\int _{0}^{1}\phi _{j}'\phi _{i}'dx+\int _{0}^{1}\phi _{j}\phi _{i}dx\!\,$ and $F_{j}=\int _{0}^{1}f\phi _{i}dx\!\,$.

### Bound error

In general terms, we can use Cea's Lemma to obtain
$\|u-u_{h}\|_{1}\leq C\inf _{v_{h}\in V_{h}}\|u-v_{h}\|_{1}$
In particular, we can consider $v_{h}\!\,$ as the Lagrange interpolant, which we denote by $v_{I}\!\,$. Then,
$\|u-u_{h}\|_{1}\leq C\|u-v_{I}\|_{1}\leq Ch\|u\|_{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:
$\|u(x)-u_{h}(x)\|_{L^{\infty }([x_{i-1},x_{i}])}\leq \max _{x\in [x_{i-1},x_{i}]}{\frac {u^{(2)}(x)}{(2)!}}\left({\frac {h}{2}}\right)^{2}$

## Problem 6c

* Derive the estimate for $\|u-u_{h}\|_{0}\!\,$, the error in $L_{2}(0,1)\!\,$. Hint: Let w solve
* (#) :${\begin{cases}-w''+w=u-u_{h},\\w(0)=w(1)=0.\end{cases}}$
* We characterize $w\!\,$ variationally as
* $w\in H_{0}^{1},\;\;B(w,v)=\int (u-u_{h})vdx,\forall v\in H_{0}^{1}\!\,$.
* Let $v=u-u_{h}\!\,$ to get
* $B(w,u-u_{h})=\int (u-u_{h})^{2}dx=\|u-u_{h}\|_{L_{2}}^{2}.\qquad \qquad (4)\!\,$
* Use formula (4) to estimate $\|u-u_{h}\|_{L_{2}}\!\,$.

## Solution 6c
### Continuity of Bilinear Form

${\begin{aligned}B(u,v)&=\int uv+\int u'v'\\&\leq \|u\|_{0}\|v\|_{0}+\|u'\|_{0}\|v'\|_{0}{\mbox{ (Cauchy-Schwarz in L2) }}\\&\leq (\|u\|_{0}^{2}+\|u'\|_{0}^{2})^{\frac {1}{2}}(\|v\|_{0}^{2}+\|v'\|^{2})^{\frac {1}{2}}{\mbox{ ( Cauchy-Schwarz in R2 ) }}\\&=\|u\|_{1}\cdot \|v\|_{1}\end{aligned}}\!\,$

### Orthogonality of the Error

$B(u-u_{h},v_{h})=0,\quad \forall v_{h}\in V_{h}\!\,$.

### Bound for L2 norm of w

${\begin{aligned}\|w\|_{0}^{2}&\leq \|w\|_{1}^{2}\\&=B(w,w)\\&=\int (u-u_{h})w\\&\leq \|u-u_{h}\|_{0}\|w\|_{0}{\mbox{ (Cauchy-Schwarz in L2 ) }}\end{aligned}}\!\,$
Hence,
$(*)\qquad \|w\|_{0}\leq \|u-u_{h}\|_{0}\!\,$

### Bound for L2 norm of w

From $(\#)\!\,$, we have
$-w''+w=u-u_{h}\!\,$
Then,
${\begin{aligned}\|w''\|_{0}&\leq \|w\|_{0}+\|u-u_{h}\|_{0}{\mbox{ (triangle inequality) }}\\&\leq \|u-u_{h}\|_{0}+\|u-u_{h}\|_{0}{\mbox{ (by (*) ) }}\\&=2\|u-u_{h}\|_{0}\end{aligned}}\!\,$

### Bound L2 Error

${\begin{aligned}\|u-u_{h}\|_{L^{2}}^{2}&=B(w,u-u_{h}){\mbox{ ( from equation (4) ) }}\\&=B(w,u-u_{h})-\underbrace {B(w_{h},u-u_{h})} _{0}\\&=B(w-w_{h},u-u_{h})\\&\leq \|w-w_{h}\|_{H^{1}(0,1)}\|u-u_{h}\|_{H^{1}(0,1)}{\mbox{ (Continuity of Bilinear Form) }}\\&\leq Ch\|w\|_{H^{2}(0,1)}\|u-u_{h}\|_{H^{1}(0,1)}{\mbox{ (Cea Lemma and Polynomial Interpolation Error) }}\end{aligned}}$
Finally, from (#), we have that $\|w\|_{H^{2}(0,1)}\leq C\|u-u_{h}\|_{L_{2}(0,1)}\!\,$. Then,
$\|u-u_{h}\|_{L^{2}(0,1)}^{2}\leq Ch\|u-u_{h}\|_{L_{2}(0,1)}\|u-u_{h}\|_{H^{1}(0,1)}\!\,$,
or equivalently,
$\|u-u_{h}\|_{L^{2}(0,1)}\leq Ch\|u-u_{h}\|_{H^{1}(0,1)}\leq Ch^{2}\|u\|_{H^{2}(0,1)}\!\,$.

## Problem 6d

* Suppose $\{\phi _{1}^{h},\dots ,\phi _{N_{h}}^{h}\}\!\,$ is a basis for $V_{h}{\mbox{ where }}N_{h}=\dim V_{h}{\mbox{, so }}u_{h}=\sum _{j=1}^{N_{h}}c_{j}^{h}\phi _{j}^{h}{\mbox{, for appropriate coefficients }}c_{j}^{h}.\!\,$. Show that
* $\|u-u_{h}\|_{1}^{2}=\|u\|_{1}^{2}-C^{T}AC\!\,$
* where $C=[c_{1}^{h},\dots ,c_{N_{h}}^{h}]^{T}{\mbox{ and }}A\!\,$ is the stiffness matrix.

## Solution 6d

We know that
${\begin{aligned}\|u-u_{h}\|_{1}^{2}&=\langle u-u_{h},u-u_{h}\rangle _{H^{1}(0,1)}\\&=\langle u-u_{h},u\rangle _{H^{1}(0,1)}-\underbrace {\langle u-u_{h},u_{h}\rangle _{H^{1}(0,1)}} _{0}\\&=\|u\|_{H^{1}(0,1)}^{2}-\langle u_{h},u\rangle _{H^{1}(0,1)}\\&=\|u\|_{H^{1}(0,1)}^{2}-\|u_{h}\|_{H^{1}(0,1)}^{2}\end{aligned}}$
where the substitution in the last lines come from the orthogonality of error.
Now, $\|u_{h}\|_{H^{1}(0,1)}^{2}=\sum _{j=1}^{n-1}\sum _{i=1}^{n-1}u_{h,j}u_{h,i}\left(\int _{0}^{1}\phi _{i}'\phi _{j}'dx+\int _{0}^{1}\phi _{i}\phi _{j}dx\right)=\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
$\|u-u_{h}\|_{1}^{2}=\|u\|_{H^{1}(0,1)}^{2}-C^{t}AC.\!\,$