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

Solution 4a
Methods {ii} and {iii} are appropriate fixed point iterations for $$ x + \ln(x) = 0 \!\,$$

{i} is not a fixed point iteration
To be a suitable fixed point iteration, the range of each iteration must be within the domain of the next iteration. In the case of {i}, $$ x_{k+1}=-\ln x_k \!\,$$ can return values in $$(-\infty,\infty) \!\,$$, but can only take x-values in $$ (0,\infty) \!\,$$.

{ii} and {iii} are fixed point iterations
Let $$ f(x):= e^{-x} \!\,$$ and $$ g(x):=\frac{x+e^{-x}}{2} \!\,$$

It is clear that

$$ f:(-\infty,\infty) \to (0,\infty) \!\,$$

and

$$ g:(-\infty,\infty) \to (0,\infty) \!\,$$

Also notice that given domain $$ D = (0,1) $$

$$ |f(x)-f(y)|\leq |f'(\eta)| |x-y| \!\,$$ for some $$ \eta \in D  \!\,$$

where $$ \max_{\eta \in D} |f'(\eta)| = \max_{\eta \in D} |e^{-x}| \leq 1 \!\,$$

and

$$ |g(x)-g(y)|\leq |g'(\eta)| |x-y| \!\,$$ for some $$ \eta \in D  \!\,$$

where $$ \max_{\eta \in D} |g'(\eta)| = \max_{\eta \in D} |\frac12 (1-e^{-x})| \approx 0.316060279 \!\,$$

That is, f,g are both contractions.

Solution 4b
For the interval (0,1), {iii} is a better fixed point iteration than {ii} since its Lipschitz constant is smaller.

Solution 4c
Newton's method gives an iterative formula that has quadratic convergence, compared to {iii}'s linear convergence.

Finding the weak form
Let $$V = H_0^1(0,1):=\{ v \in H^1 : v(0) = v(1) = 0 \} \!\,$$.

Multiplying by a test function $$ v \in V \!\,$$ and integrating by parts we obtain


 * $$ \int_0^1 (-u''v) dx = \int_0^1 u'v' dx - \left. u'v \right |_0^1 = \int_0^1 u'v' dx \!\,$$

Then the Variational formulation is: Find $$ u \in V \!\,$$ such that


 * $$ a(u,v) = \int_0^1 u'v' dx + \int_0^1 e^uv dx = l(v) = 0 \!\,$$ for all $$v \in V \!\,$$

Define piecewise linear basis functions
Let's consider a partition of the interval $$[0,1] \!\,$$,


 * $$\mathcal{\tau}_h: 0=x_0< x_1 < \dots < x_{m+1}=1 \!\,$$.

As our finite element space, we take


 * $$V_h = \{ v \in V : \left. v \right |_{[x_{i-1},x_{i}]} \in \mathcal{P}_1([x_{i-1},x_i]), i = 1,\dots, m+1, v(0)=v(1) = 0 \} \!\,$$.

For our basis of $$V_h \!\,$$, we use the hat functions $$\{ \phi_j \}_{j=1}^{m} \!\,$$, i.e., for $$j=1,2,\ldots,m \!\,$$

$$\phi_j(x)= \begin{cases} \frac{x-x_{j-1}}{h} & \mbox{for } x \in [x_{j-1},x_j] \\ \frac{x_{j+1}-x}{h} & \mbox{for } x \in [x_j, x_{j+1}] \\ 0           & \mbox{otherwise} \end{cases} \!\,$$

Therefore,

$$ \phi_j'(x)=\begin{cases} \frac1h &\mbox{for } x \in [x_{j-1},x_j] \\ -\frac1h &\mbox{for } x \in [x_j,x_{j+1}] \\ 0   &\mbox{otherwise} \end{cases} \!\,$$

Thus the discrete formulation reads: Find $$u_h \in V_h \!\,$$ such that


 * $$ a(u_h,v_h) = \int_0^1 u_h'v_h' dx + \int_0^1 e^{u_h}v_h dx = 0, \forall v_h \in V_h \!\,$$.

Since $$\{ \phi_j \}_{j=1}^{m} \!\,$$ is a basis for $$V_h \!\,$$, we have


 * $$ u_h=\sum_{j=1}^n u_{hj}\phi_j \!\,$$.

Then, we obtain the following equivalent discrete problem: Find $$ u_h=(u_{h1},\ldots, u_{hm} )^t \!\,$$ such that


 * $$ a(u_h,\phi_i) = \sum_{j=1}^m u_{hj} \int_{x_{i-1}}^{x_{i+1}} \phi_j'\phi_i' dx + \int_{x_{i-1}}^{x_{i+1}} e^{\sum_{j=1}^m u_{hj}\phi_j}\phi_i dx = 0 \mbox{ for } i=1,\dots,m\!\,$$.

Rewrite problem in matrix form
The equivalent discrete problem for $$i=1,2,\ldots,m \!\,$$

$$ \sum_{j=1}^m u_{hj}\int_{x_{i-1}}^{x_{i+1}} \phi_j'\phi_i' dx + \int_{x_{i-1}}^{x_{i+1}} e^{\sum_{j=1}^m u_{hj}\phi_j}\phi_i dx =0 \!\,$$

can be rewritten in matrix form as follows:

$$ \underbrace{ \begin{bmatrix}

\frac2h & -\frac1h   &         & \ldots \\ \ddots & \ddots     & \ddots   &   \\ &  -\frac1h & \frac2h  & -\frac1h  \\ &           &  -\frac1h& \frac 2h \end{bmatrix}}_A \underbrace{ \begin{bmatrix} u_{h1} \\ \vdots \\ \vdots \\ u_{hm} \end{bmatrix}}_{U_h}+\underbrace{h \begin{bmatrix} e^{u_{h1}} \\ e^{u_{h2}} \\ \vdots \\ e^{u_{hm}} \end{bmatrix}}_{F_h(U_h)}

\!\,$$

The first integral can be computed as follows:

$$ \begin{align} \int_0^1 \phi_j'\phi_i' &=   \begin{cases} \frac2h & \mbox{for } i=j \\ -\frac1h & \mbox{for } |i-j| =1 \\ 0    & \mbox{otherwise} \end{cases} \end{align} \!\,$$

The second integral can be approximated using the trapezoidal rule i.e.

$$ \begin{align} \int_{x_{i-1}}^{x_{i+1}} e^{u_h} \phi_i(x) &= \int_{x_{i-1}}^{x_i}e^{u_h}\phi_i(x)+\int_{x_i}^{x_{i+1}}e^{u_h}\phi_i(x)\\ &\approx \left(\frac{e^{u_h(x_{i-1})}\overbrace{\phi_i(x_{i-1})}^0+e^{u_h(x_i)}\phi_i(x_i)}{2}\right)h + \left(\frac{e^{u_h(x_{i})}\phi_i(x_{i})+e^{u_h(x_{i+1})}\overbrace{\phi_i(x_{i+1})}^0}{2}\right)h \\ &= h e^{u_h(x_i)}\overbrace{\phi_i(x_i)}^1 \\ &= h e^{u_h(x_i)}\\ &= h e^{\sum_{j=1}^m u_{hj} \phi_j(x_i)}\\ &= h e^{u_{hi}} \end{align} \!\,$$

Notice that the boundary conditions impose that $$ u_{h0}=u_{h(m+1)}=0 \!\,$$.

Local Order of Accuracy
Note that $$x_i \approx x(t_i) \!\,$$. Also, denote the uniform step size $$\Delta t \!\,$$ as h. Hence,

$$t_{i+1}=t_i + h \!\,$$

$$t_{i+2}=t_i+2h \!\,$$

Therefore, the given equation may be written as

$$x(t_i+2h)-3x(t_i+h)+2x(t_i)\approx h\left[\frac{13}{12}x'(t_i+2h)-\frac53x'(t_i+h)-\frac{5}{12}x'(t_i)\right] \!\,$$

Expand Left Hand Side Using Taylor Expansion
Expanding about $$ t_i \!\, $$, we get

$$

\begin{array}{|c|c|c|c|c|} \mbox{Order} & x(t_i+2h) & -3x(t_i+h) & 2x(t_i) & \Sigma \\ \hline & & & & \\ 0& x(t_i) & -3x(t_i) & 2x(t_i) & 0 \\ & & & & \\ 1& 2x'(t_i)h & -3x'(t_i)h & 0 & -x'(t_i)h \\ & & & & \\ 2& \frac{4}{2!}x(t_i)h^2 & -\frac{3}{2!}x(t_i)h^2 & 0 & \frac{1}{2}x''(t_i)h^2 \\ & & & & \\ 3& \frac{8}{3!}x(t_i)h^3 & -\frac{3}{3!}x(t_i)h^3 & 0 & \frac{5}{6}x'''(t_i)h^3 \\ & & & & \\ 4& \mathcal{O}(h^4) & \mathcal{O}(h^4) & 0 & \mathcal{O}(h^4) \\ & & & & \\ \hline \end{array}

$$

Expand Right Hand Side Using Taylor Expansion
Also expanding about $$ t_i \!\,$$ gives

$$

\begin{array}{|c|c|c|c|c|} \mbox{Order} & \frac{13}{12}hx'(t_i+2h) & -\frac53hx'(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} h x'(t_i) & -\frac{5}{12}hx'(t_i) & -1\cdot hx'(t_i) \\ &&&&\\ 2& \frac{13}{12}(2h) hx(t_i) & -\frac{5}{3} h h x(t_i) & 0  & \frac12 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 &  \frac43 h^3 x'''(t_i)\\ &&&&\\ 4& \mathcal{O}(h^4)& \mathcal{O}(h^4) &0 & \mathcal{O}(h^4)\\ \hline \end{array}

$$

Compare Terms to Determine Order
Taking the difference of the left and right hand side Taylor expansions show that the two step equation is of order two since the terms match up to order 2.

$$x_{i+2}-3x_{i+1}+2x_i - \Delta t \cdot \left [ \frac{13}{12} f(t_{i+2},x_{i+2})-\frac53 f(t_{i+1},x_{i+1})-\frac{5}{12}f(t_i,x_i) \right ] = \mathcal{O}(\Delta t^3) \!\,$$

Notice that the order 3 term on the left hand side differs from the order 3 term on the right hand side. ($$ \frac56 \neq \frac43 \!\,$$)

Stability Condition
Our two-step equation is stable if the roots of the equation

$$\lambda^2 - 3\lambda + 2 = 0 \!\,$$

Satisfy $$|\lambda_j| \leq 1 \!\,$$, and if $$|\lambda_j| = 1 \!\,$$, then it must be a simple root.

Since $$\lambda^2 - 3\lambda + 2 = (\lambda - 1)(\lambda - 2) = 0 \!\,$$ yields the roots 1 and 2, the two-step equation is not stable.

Convergence
A multi-step method is convergent if and only if it is stable and consistent. Our two-step equation is not stable, hence not guaranteed to converge.