Problem 4 Edit
Consider the system
x
=
1
+
h
e
−
x
2
1
+
y
2
,
y
=
.5
+
h
arctan
(
x
2
+
y
2
)
{\displaystyle x=1+h{\frac {e^{-x^{2}}}{1+y^{2}}}{\mbox{,}}\qquad y=.5+h\arctan(x^{2}+y^{2})\!\,}
.
Problem 4a Edit
Show that if the parameter
h
>
0
{\displaystyle h>0\!\,}
is chosen sufficiently small, then this system has a unique solution
(
x
∗
,
y
∗
)
{\displaystyle (x^{*},y^{*})\!\,}
within some rectangular region.
Solution 4a Edit
The system of equations may be expressed in matrix notation.
f
(
x
,
y
)
=
[
1
+
h
e
−
x
2
1
+
y
2
.5
+
h
arctan
(
x
2
+
y
2
)
]
{\displaystyle 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
)
{\displaystyle f(x,y)\!\,}
,
J
(
f
)
{\displaystyle J(f)\!\,}
, is computed using partial deriatives
J
(
f
)
=
[
h
1
+
y
2
e
−
x
2
(
−
2
x
)
−
h
2
y
e
−
x
2
(
1
+
y
2
)
2
h
⋅
2
x
1
+
(
x
2
+
y
2
)
2
h
⋅
2
y
1
+
(
x
2
+
y
2
)
2
]
{\displaystyle 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
{\displaystyle h\!\,}
is sufficiently small and
x
,
y
{\displaystyle x,y\!\,}
are restricted to a bounded region
B
{\displaystyle B\!\,}
,
‖
J
(
f
)
‖
⏟
C
<
1
{\displaystyle \underbrace {\|J(f)\|} _{C}<1\!\,}
From the mean value theorem, for
x
1
→
,
x
2
→
∈
B
{\displaystyle {\vec {x_{1}}},{\vec {x_{2}}}\in B\!\,}
, there exists
ξ
→
{\displaystyle {\vec {\xi }}\!\,}
such that
f
(
x
1
→
)
−
f
(
x
2
→
)
=
J
(
f
(
ξ
→
)
)
(
x
1
→
−
x
2
→
)
{\displaystyle f({\vec {x_{1}}})-f({\vec {x_{2}}})=J(f({\vec {\xi }}))({\vec {x_{1}}}-{\vec {x_{2}}})\!\,}
Since
‖
J
(
f
)
‖
{\displaystyle \|J(f)\|\!\,}
is bounded in the region
B
{\displaystyle B\!\,}
give sufficiently small
h
{\displaystyle h\!\,}
‖
f
(
x
1
→
)
−
f
(
x
2
→
)
‖
≤
C
‖
x
1
→
−
x
2
→
‖
{\displaystyle \|f({\vec {x_{1}}})-f({\vec {x_{2}}})\|\leq C\|{\vec {x_{1}}}-{\vec {x_{2}}}\|\!\,}
Therefore,
f
{\displaystyle f\!\,}
is a contraction and from the contraction mapping theorem there exists an unique fixed point (solution) in a rectangular region.
Problem 4b Edit
Derive a fixed point iteration scheme for solving the system and show that it converges.
{\displaystyle \!\,}
Solution 4b Edit
Use Newton Method Edit
To solve this problem, we can use the Newton Method. In fact, we want to find the zeros of
g
(
x
,
y
)
=
[
x
y
]
−
[
1
+
h
e
−
x
2
1
+
y
2
.5
+
h
arctan
(
x
2
+
y
2
)
]
{\displaystyle 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
)
{\displaystyle g(x,y)\!\,}
,
J
(
g
)
{\displaystyle J(g)\!\,}
, is computed using partial deriatives
J
(
g
)
=
[
1
−
h
1
+
y
2
e
−
x
2
(
−
2
x
)
−
h
e
−
x
2
1
+
y
2
h
⋅
2
x
1
+
(
x
2
+
y
2
)
2
1
−
h
⋅
2
y
1
+
(
x
2
+
y
2
)
2
]
{\displaystyle 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
[
x
k
+
1
y
k
+
1
]
=
[
x
k
y
k
]
−
1
d
e
t
(
J
(
g
)
)
[
1
−
h
⋅
2
y
1
+
(
x
2
+
y
2
)
2
h
e
−
x
2
1
+
y
2
−
h
⋅
2
x
1
+
(
x
2
+
y
2
)
2
1
−
h
1
+
y
2
e
−
x
2
(
−
2
x
)
]
⏟
J
(
g
)
−
1
g
(
x
k
,
y
k
)
{\displaystyle {\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 Edit
Jacobian of g is Lipschitz Edit
We want to show that
J
(
g
)
{\displaystyle J(g)\!\,}
is a Lipschitz function. In fact,
‖
J
(
g
)
(
x
1
→
)
−
J
(
g
)
(
x
2
→
)
‖
=
‖
(
I
−
J
(
f
)
)
(
x
1
→
)
−
(
I
−
J
(
f
)
)
(
x
2
→
)
‖
=
‖
x
1
→
−
x
2
→
−
J
(
f
)
x
1
→
+
J
(
f
)
x
2
→
‖
≤
‖
x
1
→
−
x
2
→
‖
+
‖
J
(
f
)
x
2
→
−
J
(
f
)
x
1
→
‖
{\displaystyle {\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
)
{\displaystyle J(f)\!\,}
is Lipschitz, we have
‖
J
(
g
)
(
x
1
→
)
−
J
(
g
)
(
x
2
→
)
‖
≤
(
1
+
C
)
‖
(
x
1
→
)
−
(
x
2
→
)
‖
{\displaystyle \|J(g)({\vec {x_{1}}})-J(g)({\vec {x_{2}}})\|\leq (1+C)\|({\vec {x_{1}}})-({\vec {x_{2}}})\|\!\,}
Jacobian of g is invertible Edit
Since
J
(
f
)
{\displaystyle J(f)\!\,}
is a contraction, the spectral radius of the Jacobian of
f
{\displaystyle f\!\,}
is less than 1 i.e.
ρ
(
J
(
f
)
)
<
1
{\displaystyle \rho (J(f))<1\!\,}
.
On the other hand, we know that the eigenvalues of
J
(
g
)
{\displaystyle J(g)\!\,}
are
λ
(
J
(
g
)
)
=
1
−
λ
(
J
(
f
)
)
{\displaystyle \lambda (J(g))=1-\lambda (J(f))\!\,}
.
Then, it follows that
d
e
t
(
I
−
(
J
(
g
)
)
)
≠
0
{\displaystyle det(I-(J(g)))\neq 0\!\,}
or equivalently
J
(
g
)
{\displaystyle J(g)\!\,}
is invertible.
inverse of Jacobian of g is bounded Edit
1
d
e
t
(
J
(
g
)
)
[
1
−
h
⋅
2
y
1
+
(
x
2
+
y
2
)
2
h
e
−
x
2
1
+
y
2
−
h
⋅
2
x
1
+
(
x
2
+
y
2
)
2
1
−
h
1
+
y
2
e
−
x
2
(
−
2
x
)
]
⏟
J
(
g
)
−
1
{\displaystyle \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
{\displaystyle J(g)^{-1}\!\,}
exists,
det
(
J
(
g
)
)
≠
0
{\displaystyle {\mbox{det}}(J(g))\neq 0\!\,}
. Given a bounded region (bounded
x
,
y
{\displaystyle 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 Edit
‖
J
(
g
)
−
1
(
x
0
)
∗
g
(
x
0
)
‖
<
∞
{\displaystyle \|J(g)^{-1}(x_{0})*g(x_{0})\|<\infty \!\,}
since
J
(
g
)
−
1
{\displaystyle J(g)^{-1}\!\,}
and
g
(
x
)
{\displaystyle g(x)\!\,}
is bounded.
Then, using a good enough approximation
(
x
0
,
y
0
)
{\displaystyle (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
∗
)
‖
≤
K
‖
(
x
k
,
y
k
)
−
(
x
∗
,
y
∗
)
‖
2
.
{\displaystyle \|(x_{k-1},y_{k+1})-(x^{*},y^{*})\|\leq K\|(x_{k},y_{k})-(x^{*},y^{*})\|^{2}.\!\,}
Problem 5a Edit
Outline the derivation of the Adams-Bashforth methods for the numerical solution of the initial value problem
y
′
=
f
(
x
,
y
)
,
y
(
x
0
)
=
y
0
{\displaystyle y'=f(x,y){\mbox{, }}y(x_{0})=y_{0}\!\,}
.
Solution 5a Edit
We want to solve the following initial value problem:
y
′
=
f
(
x
,
y
)
,
y
(
x
0
)
=
y
0
{\displaystyle y'=f(x,y){\mbox{, }}y(x_{0})=y_{0}\!\,}
.
First, we integrate this expression over
[
t
n
,
t
n
+
1
]
{\displaystyle [t_{n},t_{n+1}]\!\,}
, to obtain
y
(
t
n
+
1
)
=
y
(
t
n
)
+
∫
t
n
t
n
+
1
f
(
t
,
y
(
t
)
)
d
t
{\displaystyle 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
)
)
{\displaystyle f(t,y(t))\!\,}
using an appropriate interpolation polynomial of degree
p
{\displaystyle p\!\,}
at
t
n
,
t
n
−
1
,
.
.
.
,
t
n
−
p
{\displaystyle t_{n},t_{n-1},...,t_{n-p}\!\,}
.
This idea generates the Adams-Bashforth methods.
y
n
+
1
=
y
n
+
∫
t
n
t
n
−
1
∑
k
=
0
p
g
(
t
n
−
k
)
l
k
(
t
)
d
t
{\displaystyle 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
{\displaystyle y_{n}\!\,}
denotes the approximated solution,
g
(
t
)
≡
f
(
t
,
y
)
{\displaystyle g(t)\equiv f(t,y)\!\,}
and
l
k
(
t
)
{\displaystyle l_{k}(t)\!\,}
denotes the associated Lagrange polynomial.
Problem 5b Edit
Derive the Adams-Bashforth formula
y
i
+
1
=
y
i
+
h
[
−
1
2
f
(
x
i
−
1
,
y
i
−
1
)
+
3
2
f
(
x
i
,
y
i
)
]
(
1
)
{\displaystyle 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 Edit
From (a) we have
y
i
+
1
=
y
i
+
∫
x
i
x
i
+
1
f
d
x
{\displaystyle y_{i+1}=y_{i}+\int _{x_{i}}^{x_{i+1}}fdx\!\,}
where
∫
x
i
x
i
+
1
f
d
x
≈
∫
x
i
x
i
+
1
[
x
−
x
i
x
i
−
1
−
x
i
f
i
−
1
+
x
−
x
i
−
1
x
i
−
x
i
−
1
f
i
]
d
x
{\displaystyle \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
+
s
h
{\displaystyle x=x_{i}+sh\!\,}
, where h is the fixed step size we get
d
x
=
h
d
s
x
−
x
i
=
s
h
x
−
x
i
−
1
=
(
1
+
s
)
h
x
i
−
x
i
−
1
=
h
{\displaystyle {\begin{aligned}dx&=hds\\x-x_{i}&=sh\\x-x_{i-1}&=(1+s)h\\x_{i}-x_{i-1}&=h\end{aligned}}}
h
∫
0
1
[
s
h
−
h
f
i
−
1
+
(
1
+
s
)
h
h
f
i
]
d
s
=
h
[
−
1
2
f
i
−
1
+
3
2
f
i
]
{\displaystyle 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
[
−
1
2
f
(
x
i
−
1
,
y
i
−
1
)
+
3
2
f
(
x
i
,
y
i
)
]
{\displaystyle 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 Edit
Analyze the method (1). To be specific, find the local truncation error, prove convergence and find the order of convergence.
Solution 5c Edit
Find Local Truncation Error Using Taylor Expansion Edit
Note that
y
i
≈
y
(
t
i
)
{\displaystyle y_{i}\approx y(t_{i})\!\,}
. Also, denote the uniform step size
Δ
t
{\displaystyle \Delta t\!\,}
as h. Hence,
t
i
−
1
=
t
i
−
h
{\displaystyle t_{i-1}=t_{i}-h\!\,}
t
i
+
1
=
t
i
+
h
{\displaystyle t_{i+1}=t_{i}+h\!\,}
Therefore, the given equation may be written as
y
(
t
i
+
h
)
−
y
(
t
i
)
≈
h
[
−
1
2
y
′
(
t
i
−
h
)
+
3
2
y
′
(
t
i
)
]
{\displaystyle 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 Edit
Expanding about
t
i
{\displaystyle t_{i}\!\,}
, we get
Order
y
(
t
i
+
h
)
y
(
t
i
)
Σ
0
y
(
t
i
)
y
(
t
i
)
0
1
−
y
′
(
t
i
)
h
0
−
y
′
(
t
i
)
h
2
1
2
!
y
″
(
t
i
)
h
2
0
1
2
y
″
(
t
i
)
h
2
3
−
1
3
!
y
‴
(
t
i
)
h
3
0
−
1
6
y
‴
(
t
i
)
h
3
4
O
(
h
4
)
0
O
(
h
4
)
{\displaystyle {\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 Edit
Also expanding about
t
i
{\displaystyle t_{i}\!\,}
gives
Order
13
12
h
x
′
(
t
i
+
2
h
)
−
5
3
h
x
′
(
t
i
+
h
)
−
5
12
h
x
′
(
t
i
)
Σ
0
0
0
0
0
1
13
12
h
x
′
(
t
i
)
−
5
3
h
x
′
(
t
i
)
−
5
12
h
x
′
(
t
i
)
−
1
⋅
h
x
′
(
t
i
)
2
13
12
(
2
h
)
h
x
″
(
t
i
)
−
5
3
h
h
x
″
(
t
i
)
0
1
2
h
2
x
″
(
t
i
)
3
13
12
(
2
h
)
2
h
x
‴
(
t
i
)
2
−
5
3
h
2
h
x
‴
(
t
i
)
2
0
4
3
h
3
x
‴
(
t
i
)
4
O
(
h
4
)
O
(
h
4
)
0
O
(
h
4
)
{\displaystyle {\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 Edit
A method is convergent if and only if it is both stable and consistent.
Stability Edit
It is easy to show that the method is zero stable as it satisfies the root condition. So stable.
Consistency Edit
Truncation error is of order 2. Truncation error tends to zero as h tends to zeros. So the method is consistent.
Order of Convergence Edit
Dahlquist principle: consistency + stability = convergent. Order of convergence will be 2.
Problem 6 Edit
Consider the problem
−
u
″
+
u
=
f
(
x
)
,
0
≤
x
≤
1
,
u
(
0
)
=
u
(
1
)
=
0.
(
2
)
{\displaystyle -u''+u=f(x),\;\;0\leq x\leq 1,\;\;u(0)=u(1)=0.\qquad \qquad (2)}
Problem 6a Edit
Give a variational formulation of (2), i.e. express (2) as
u
∈
H
,
B
(
x
,
y
)
=
F
(
v
)
,
∀
v
∈
H
(
3
)
{\displaystyle 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 Edit
Multiplying (2) by a test function and using integration by parts we obtain:
−
∫
0
1
u
″
v
d
x
+
∫
0
1
u
v
d
x
=
∫
0
1
f
v
d
x
{\displaystyle -\int _{0}^{1}u''vdx+\int _{0}^{1}uvdx=\int _{0}^{1}fvdx\!\,}
−
u
′
v
|
0
1
⏟
0
+
∫
0
1
u
′
v
′
d
x
+
∫
0
1
u
v
d
x
=
∫
0
1
f
v
d
x
{\displaystyle \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
∈
H
{\displaystyle u\in H\!\,}
such that
B
(
u
,
v
)
=
∫
0
1
u
′
v
′
d
x
+
∫
0
1
u
v
d
x
=
∫
0
1
f
v
d
x
=
F
(
v
)
{\displaystyle B(u,v)=\int _{0}^{1}u'v'dx+\int _{0}^{1}uvdx=\int _{0}^{1}fvdx=F(v)\!\,}
for all
v
∈
H
,
{\displaystyle v\in H,\!\,}
where
H
:=
H
0
1
(
0
,
1
)
=
{
v
:
v
∈
H
1
,
v
(
0
)
=
v
(
1
)
=
0
}
{\displaystyle H:=H_{0}^{1}(0,1)=\{v:v\in H^{1},v(0)=v(1)=0\}\!\,}
.
Problem 6b Edit
Let
0
=
x
0
<
x
1
<
⋯
<
x
n
=
1
{\displaystyle 0=x_{0}<x_{1}<\cdots <x_{n}=1\!\,}
be a mesh on
[
0
,
1
]
{\displaystyle [0,1]\!\,}
with
h
=
max
j
=
0
,
1
,
…
,
n
−
1
(
x
j
+
1
−
x
j
)
{\displaystyle h=\max _{j=0,1,\dots ,n-1}(x_{j+1}-x_{j})\!\,}
, and let
V
h
=
{
u
:
u
continuous on
[
0
,
1
]
,
u
|
[
x
j
,
x
j
+
1
]
is linear for each
j
,
u
(
0
)
=
u
(
1
)
=
0
}
{\displaystyle 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
to
u
{\displaystyle u_{h}{\mbox{ to }}u\!\,}
based on the approximation space
V
h
{\displaystyle V_{h}\!\,}
. What can be said about
‖
u
−
u
h
‖
1
{\displaystyle \|u-u_{h}\|_{1}\!\,}
the error on the Sobolev norm on
H
1
(
0
,
1
)
{\displaystyle H^{1}(0,1)\!\,}
?
Solution 6b Edit
Define piecewise linear basis Edit
For our basis of
V
h
{\displaystyle V_{h}\!\,}
, we use the set of hat functions
{
ϕ
j
}
j
=
1
n
−
1
{\displaystyle \{\phi _{j}\}_{j=1}^{n-1}\!\,}
, i.e., for
j
=
1
,
2
,
…
,
n
−
1
{\displaystyle j=1,2,\ldots ,n-1\!\,}
ϕ
j
(
x
)
=
{
x
−
x
j
−
1
x
j
−
x
j
−
1
for
x
∈
[
x
j
−
1
,
x
j
]
x
j
+
1
−
x
x
j
+
1
−
x
j
for
x
∈
[
x
j
,
x
j
+
1
]
0
otherwise
{\displaystyle \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 Edit
Since
{
ϕ
j
}
j
=
1
n
−
1
{\displaystyle \{\phi _{j}\}_{j=1}^{n-1}\!\,}
is a basis for
V
h
{\displaystyle V_{h}\!\,}
, and
u
h
∈
V
h
{\displaystyle u_{h}\in V_{h}\!\,}
we have
u
h
=
∑
j
=
1
n
−
1
u
h
j
ϕ
j
{\displaystyle u_{h}=\sum _{j=1}^{n-1}u_{hj}\phi _{j}\!\,}
.Discrete Problem Edit
Now, we can write the discrete problem: Find
u
h
∈
V
h
{\displaystyle u_{h}\in V_{h}\!\,}
such that
B
(
u
h
,
v
h
)
=
F
(
v
h
)
{\displaystyle B(u_{h},v_{h})=F(v_{h})\!\,}
for all
v
h
∈
V
h
.
{\displaystyle v_{h}\in V_{h}.\!\,}
If we consider that
{
ϕ
j
}
j
=
1
n
−
1
{\displaystyle \{\phi _{j}\}_{j=1}^{n-1}\!\,}
is a basis of
V
h
{\displaystyle V_{h}\!\,}
and the linearity of the bilinear form
B
{\displaystyle B\!\,}
and the functional
F
{\displaystyle F\!\,}
, we obtain the equivalent problem:
Find
u
h
=
(
u
h
,
1
,
…
,
u
h
,
n
−
1
)
∈
R
n
−
1
{\displaystyle u_{h}=(u_{h,1},\dots ,u_{h,{n-1}})\in \mathbb {R} ^{n-1}\!\,}
such that
∑
j
=
1
n
−
1
u
h
,
j
{
∫
0
1
ϕ
j
′
ϕ
i
′
d
x
+
∫
0
1
ϕ
j
ϕ
i
d
x
}
=
∫
0
1
f
ϕ
i
d
x
,
i
=
1
,
…
,
n
−
1.
{\displaystyle \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
,
…
,
u
h
,
n
−
1
)
∈
R
n
−
1
{\displaystyle u_{h}=(u_{h,1},\dots ,u_{h,{n-1}})\in \mathbb {R} ^{n-1}\!\,}
such that
A
u
h
=
F
,
{\displaystyle Au_{h}=F,\!\,}
where
A
i
j
=
∫
0
1
ϕ
j
′
ϕ
i
′
d
x
+
∫
0
1
ϕ
j
ϕ
i
d
x
{\displaystyle A_{ij}=\int _{0}^{1}\phi _{j}'\phi _{i}'dx+\int _{0}^{1}\phi _{j}\phi _{i}dx\!\,}
and
F
j
=
∫
0
1
f
ϕ
i
d
x
{\displaystyle F_{j}=\int _{0}^{1}f\phi _{i}dx\!\,}
.
Bound error Edit
In general terms, we can use Cea's Lemma to obtain
‖
u
−
u
h
‖
1
≤
C
inf
v
h
∈
V
h
‖
u
−
v
h
‖
1
{\displaystyle \|u-u_{h}\|_{1}\leq C\inf _{v_{h}\in V_{h}}\|u-v_{h}\|_{1}}
In particular, we can consider
v
h
{\displaystyle v_{h}\!\,}
as the Lagrange interpolant, which we denote by
v
I
{\displaystyle v_{I}\!\,}
. Then,
‖
u
−
u
h
‖
1
≤
C
‖
u
−
v
I
‖
1
≤
C
h
‖
u
‖
H
2
(
0
,
1
)
{\displaystyle \|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
∞
(
[
x
i
−
1
,
x
i
]
)
≤
max
x
∈
[
x
i
−
1
,
x
i
]
u
(
2
)
(
x
)
(
2
)
!
(
h
2
)
2
{\displaystyle \|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 Edit
Derive the estimate for
‖
u
−
u
h
‖
0
{\displaystyle \|u-u_{h}\|_{0}\!\,}
, the error in
L
2
(
0
,
1
)
{\displaystyle L_{2}(0,1)\!\,}
. Hint: Let w solve
(#) :
{
−
w
″
+
w
=
u
−
u
h
,
w
(
0
)
=
w
(
1
)
=
0.
{\displaystyle {\begin{cases}-w''+w=u-u_{h},\\w(0)=w(1)=0.\end{cases}}}
We characterize
w
{\displaystyle w\!\,}
variationally as
w
∈
H
0
1
,
B
(
w
,
v
)
=
∫
(
u
−
u
h
)
v
d
x
,
∀
v
∈
H
0
1
{\displaystyle w\in H_{0}^{1},\;\;B(w,v)=\int (u-u_{h})vdx,\forall v\in H_{0}^{1}\!\,}
.Let
v
=
u
−
u
h
{\displaystyle v=u-u_{h}\!\,}
to get
B
(
w
,
u
−
u
h
)
=
∫
(
u
−
u
h
)
2
d
x
=
‖
u
−
u
h
‖
L
2
2
.
(
4
)
{\displaystyle 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
{\displaystyle \|u-u_{h}\|_{L_{2}}\!\,}
.
Solution 6c Edit
Continuity of Bilinear Form Edit
B
(
u
,
v
)
=
∫
u
v
+
∫
u
′
v
′
≤
‖
u
‖
0
‖
v
‖
0
+
‖
u
′
‖
0
‖
v
′
‖
0
(Cauchy-Schwartz in L2)
≤
(
‖
u
‖
0
2
+
‖
u
′
‖
0
2
)
1
2
(
‖
v
‖
0
2
+
‖
v
′
‖
2
)
1
2
( Cauchy-Schwartz in R2 )
=
‖
u
‖
1
⋅
‖
v
‖
1
{\displaystyle {\begin{aligned}B(u,v)&=\int uv+\int u'v'\\&\leq \|u\|_{0}\|v\|_{0}+\|u'\|_{0}\|v'\|_{0}{\mbox{ (Cauchy-Schwartz in L2) }}\\&\leq (\|u\|_{0}^{2}+\|u'\|_{0}^{2})^{\frac {1}{2}}(\|v\|_{0}^{2}+\|v'\|^{2})^{\frac {1}{2}}{\mbox{ ( Cauchy-Schwartz in R2 ) }}\\&=\|u\|_{1}\cdot \|v\|_{1}\end{aligned}}\!\,}
Orthogonality of the Error Edit
B
(
u
−
u
h
,
v
h
)
=
0
,
∀
v
h
∈
V
h
{\displaystyle B(u-u_{h},v_{h})=0,\quad \forall v_{h}\in V_{h}\!\,}
.
Bound for L2 norm of w Edit
‖
w
‖
0
2
≤
‖
w
‖
1
2
=
B
(
w
,
w
)
=
∫
(
u
−
u
h
)
w
≤
‖
u
−
u
h
‖
0
‖
w
‖
0
(Cauchy-Schwartz in L2 )
{\displaystyle {\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-Schwartz in L2 ) }}\end{aligned}}\!\,}
Hence,
(
∗
)
‖
w
‖
0
≤
‖
u
−
u
h
‖
0
{\displaystyle (*)\qquad \|w\|_{0}\leq \|u-u_{h}\|_{0}\!\,}
Bound for L2 norm of w Edit
From
(
#
)
{\displaystyle (\#)\!\,}
, we have
−
w
″
+
w
=
u
−
u
h
{\displaystyle -w''+w=u-u_{h}\!\,}
Then,
‖
w
″
‖
0
≤
‖
w
‖
0
+
‖
u
−
u
h
‖
0
(triangle inequality)
≤
‖
u
−
u
h
‖
0
+
‖
u
−
u
h
‖
0
(by (*) )
=
2
‖
u
−
u
h
‖
0
{\displaystyle {\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 Edit
‖
u
−
u
h
‖
L
2
2
=
B
(
w
,
u
−
u
h
)
( from equation (4) )
=
B
(
w
,
u
−
u
h
)
−
B
(
w
h
,
u
−
u
h
)
⏟
0
=
B
(
w
−
w
h
,
u
−
u
h
)
≤
‖
w
−
w
h
‖
H
1
(
0
,
1
)
‖
u
−
u
h
‖
H
1
(
0
,
1
)
(Continuity of Bilinear Form)
≤
C
h
‖
w
‖
H
2
(
0
,
1
)
‖
u
−
u
h
‖
H
1
(
0
,
1
)
(Cea Lemma and Polynomial Interpolation Error)
{\displaystyle {\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
)
≤
C
‖
u
−
u
h
‖
L
2
(
0
,
1
)
{\displaystyle \|w\|_{H^{2}(0,1)}\leq C\|u-u_{h}\|_{L_{2}(0,1)}\!\,}
. Then,
‖
u
−
u
h
‖
L
2
(
0
,
1
)
2
≤
C
h
‖
u
−
u
h
‖
L
2
(
0
,
1
)
‖
u
−
u
h
‖
H
1
(
0
,
1
)
{\displaystyle \|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
)
≤
C
h
‖
u
−
u
h
‖
H
1
(
0
,
1
)
≤
C
h
2
‖
u
‖
H
2
(
0
,
1
)
{\displaystyle \|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 Edit
Suppose
{
ϕ
1
h
,
…
,
ϕ
N
h
h
}
{\displaystyle \{\phi _{1}^{h},\dots ,\phi _{N_{h}}^{h}\}\!\,}
is a basis for
V
h
where
N
h
=
dim
V
h
, so
u
h
=
∑
j
=
1
N
h
c
j
h
ϕ
j
h
, for appropriate coefficients
c
j
h
.
{\displaystyle 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
A
C
{\displaystyle \|u-u_{h}\|_{1}^{2}=\|u\|_{1}^{2}-C^{T}AC\!\,}
where
C
=
[
c
1
h
,
…
,
c
N
h
h
]
T
and
A
{\displaystyle C=[c_{1}^{h},\dots ,c_{N_{h}}^{h}]^{T}{\mbox{ and }}A\!\,}
is the stiffness matrix.
Solution 6d Edit
We know that
‖
u
−
u
h
‖
1
2
=
⟨
u
−
u
h
,
u
−
u
h
⟩
H
1
(
0
,
1
)
=
⟨
u
−
u
h
,
u
⟩
H
1
(
0
,
1
)
−
⟨
u
−
u
h
,
u
h
⟩
H
1
(
0
,
1
)
⏟
0
=
‖
u
‖
H
1
(
0
,
1
)
2
−
⟨
u
h
,
u
⟩
H
1
(
0
,
1
)
=
‖
u
‖
H
1
(
0
,
1
)
2
−
‖
u
h
‖
H
1
(
0
,
1
)
2
{\displaystyle {\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
=
∑
j
=
1
n
−
1
∑
i
=
1
n
−
1
u
h
,
j
u
h
,
i
(
∫
0
1
ϕ
i
′
ϕ
j
′
d
x
+
∫
0
1
ϕ
i
ϕ
j
d
x
)
=
∑
j
=
1
n
−
1
∑
i
=
1
n
−
1
u
h
,
j
A
i
j
u
h
,
i
=
C
t
A
C
.
{\displaystyle \|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
A
C
.
{\displaystyle \|u-u_{h}\|_{1}^{2}=\|u\|_{H^{1}(0,1)}^{2}-C^{t}AC.\!\,}