爱的深处不正涌流着企图和对方分毫不差的这种不可能实现的热望吗？这种热望不正驱使着人们将不可能由另一极端变成可能，从而引导他们走向那种悲剧性的叛离之路吗？既然相爱的人不可能完全相似，那么不如干脆使他们致力于相互之间毫不相似，使这样的叛离完整地作用于媚态。

# 第一章 相关基础数学知识

Theorem: If the function $f$ is differentiable at $x_0$, then $f$ is continuous at $x_0$.

Notation: $C[a,b]$ stands for the set of continuous function defined in $[a,b]$.

Theorem (Generalized Rolle’s Theorem): Suppose $f\in C[a,b]$ is $n$ times differentiable on $(a,b)$. If $f(x)$ is zero at the $n+1$ distinct numbers $x_0,x_1,...,x_n$ in the $[a,b]$, then a number $c$ in the $(a,b)$ exists with $f^{(n)}(c)=0$

Theorem (Taylor’s Theorem): $f(x)=P_n(x)+R_n(x)$, $P_n(x)=\sum_{k=0}^n\frac{f^{(k)}(x)}{k!}(x-x_0)^{k}$, $R_n(x)=\frac{f^{(n+1)(\xi(x))}}{(n+1)!}(x-x_0)^{n+1}$. $\xi(x)$ is between $x$ and $x_0$.

---

Definition: If $p^*$ is an approximation to $p$, the absolute error is $|p-p^*|$, and the relative error is $\frac{|p-p^*|}{|p|}$, provided that $p\neq 0$.

Definition: The number $p^*$ is said to approximate $p$ to $t$ significant digit if $t$ is the largest nonnegative integer for which $\frac{|p-p^*|}{|p|}<5\times 10^{-t}$.

---

Definition: An algorithm is said to be stable when small changes in the initial data can produce correspondingly small changes in final results. Some algorithm are stable only for certain choices of initial data, these cases are called conditionally stable.

Definition: Suppose that $E_0$ denotes an initial error and $E_n$ represents the magnitude of an error after $n$ subsequent operations.

• If $E_n\approx nCE_0$, where $C$ is a constant independent of $n$, then the growth of error is said to be linear.

• If $E_n\approx C^nE_0$, for some $C>1$, then the growth of error is said to be exponential.

Definition: Suppose $\{\beta_n\}^{\infty}_{n=1}$ is a sequence known to converge to zero, and $\{\alpha_n\}^{\infty}_{n=1}$ converges to a number $\alpha$.

• If a positive constant $K$ exists with $|\alpha_n-\alpha|\leq K|\beta_n|$ for large n, then we say that $\{\alpha\}^\infty_{n=1}$ converges to $\alpha$ with the rate of convergence $O(\beta_n)$, writing $\alpha_n=\alpha+O(\beta_n)$.

Property (Operator errors):

• $\alpha_n+\beta_n=\alpha+\beta+O(\varepsilon_\alpha+\varepsilon_\beta)$.
• $\alpha_n\beta_n=\alpha\beta+\alpha O(\varepsilon_\beta)+\beta O(\varepsilon_\alpha)+O(\varepsilon_\alpha\varepsilon_\beta)$.

Example: $\alpha_n=\frac{n+1}{n^2},\beta_n=\frac{1}{n}$. $|\alpha_n-0|<\frac{n+n}{n^2}=2\times\frac{1}{n}$. So $\alpha_n=0+O(\frac{1}{n})$.

Definition: Suppose

$$\lim_{x\rightarrow 0}G(x)=0,\lim_{x\rightarrow 0}F(x)=L$$

If a positive constant $K$ exists with $|F(x)-L|\leq K|G(x)|$, then we write $F(x)=L+O(G(x))$.

Example: $cos(x)+\frac{1}{2}x^2=1+O(x^4)$.

---

Definition: A mathematical problem is said to be well-posed if a solution:

• exists,
• is unique,
• depends continuously on input data

Otherwise, the problem is ill-posed.

# 第二章 一元方程的近似解

Definition: The system of $m$ nonlinear equations in $n$ unknowns can alternatively be represented by defining a function $f$, mapping $\mathbb{R}^n$ into $\mathbb{R}^m$ by

$$\textbf{f}(\textbf{x})=(f_1(\textbf{x}),...,f_m(\textbf{x}))^T$$

• The function $f_1,f_2,..,f_m$ are the coordinate functions of f.

• The function f is continuous at $x_0\in D$ provided $\lim_{x\rightarrow x_0}\textbf{f}(\textbf{x})=\textbf{f}(\textbf{x}_0)$.

Theorem: Let $f$ be a function from $D\subset \mathbb{R}^n$ into $\mathbb{R}$ and $x_0\in D$. If

$$|\frac{\partial f(\textbf{x})}{\partial x_j}|\leq K,j=1,2,...,n$$

whenever $|x_j-x_{0j}|\leq\delta$, then $f$ is continuous at $\textbf{x}_0$.

# Bisection (Binary search).

Theorem: Suppose that $f\in C[a,b]$ and $f(a)f(b)<0$. The Bisection method generates a sequence $\{p_n\}^\infty_1$ approximating a zero point $p$ of $f$ ($f(p)=0$) with

$$|p_n-p|\leq\frac{b-a}{2^ n}$$

So $p_n=p+O(\frac{1}{2^n})$.

# Fixed Point Method

Definition: Fixed point of given function $g:\mathbb{R}\rightarrow\mathbb{R}$ is value $x^*$ such that $x^*=g(x^*)$

---

For given equation $f(x)=0$, there may be many equivalent fixed point problems $x=g(x)$ with different choices for g.

Example: $f(x)=0\Leftrightarrow g(x)=x$, we can choose $g(x)=x-f(x)$ or $g(x)=x+3f(x)$.

To approximate the fixed point of a function $g(x)$, we choose an initial approximation $p_0$, and generate the sequence by:

$$p_n=g(p_{n-1}),n=1,2,3,...$$

If the sequence converges to $p$ and $g(x)$ is continuous, then we have

$$p=\lim_{n\rightarrow\infty}p_n=\lim_{n\rightarrow\infty}g(p_n)=g(\lim_{n\rightarrow \infty}p_n)=g(p)$$

This technique is called fixed point iteration (or functional iteration).

---

Theorem: If $g(x)\in C[a,b]$ and $g(x)\in[a,b]$ for all $x\in[a,b]$, then $g(x)$ has a fixed point in $[a,b]$. What’s more, if $g'(x)$ exists on $(a,b)$ and a positive constant $k<1$ exists with $|g'(x)|\leq k$ for all $x\in(a,b)$, then the fixed point is unique.

• The proof is easily obtained by proof by contradiction considering the geometric feature.

Corollary: Obviously we have $p_n-p=|g(p_{n-1})-g(p)|=|g'(\xi)||p_{n-1}-p|$ where $\xi\in(a,b)$. So $|p_n-p|\leq k^n|p_0-p|$. If $k<1$, the fixed point is unique and we have $p_n=p+O(k^n)$.

What’s more, $|p_n-p|\leq k^n|p_0-p|\leq k^nmax(p_0-a,b-p_0)$, and

$$|p_n-p_0|\leq |p_n-p_{n-1}|+|p_{n-1}-p_{n-2}|+...+|p_1-p_0|\leq \frac{k^n}{1-k}|p_1-p_0|$$

# Newton’s Method

Suppose that $f\in C^2[a,b]$ and $x^*$ is a solution for $f(x)=0$. Let $\bar{x}\in[a,b]$ e an approximation to $x^*$ such that $f'(\bar{x})\neq 0$ and $|\bar{x}-x^*|$ is very small. Due to the Taylor polynomial:

$$f(x)=f(\bar{x})+(x-\bar{x})f'(\bar{x})+\frac{(x-\bar{x})^2}{2}f''(\xi)$$

where $\xi$ lies between $x$ and $\bar{x}$. So

$$0=f(x^*)\approx f(\bar{x})+(x^*-\bar{x})f'(\bar{x})\\ x^*\approx\bar{x}-\frac{f(\bar{x})}{f'(\bar{x})}$$

Theorem: Let $f\in C^2[a,b]$, if $p\in[a,b]$ is such that $f(p)=0$ and $f'(p)\neq 0$, then there exists a $\delta>0$ such that for any $p_0\in[p-\delta,p+\delta]$, $p_n=p_{n-1}-\frac{f(p_{n-1})}{f'(p_{n-1})}$ converges to $p$.

• The proof is based on the fixed point iteration $g(x)=x-\frac{f(x)}{f'(x)}$.

# Secant Method

Similar to Newton’s method:

$$p_{n}=p_{n-1}-\frac{f(p_{n-1})}{\frac{f(p_{n-1})-f(p_{n-2})}{p_{n-1}-p_{n-2}}}$$

# Method of False Position

Improved Secant Method. If $f(p_0)f(p_1)<0$, we can chose:

$$p_n=\begin{cases} p_{n-1}-\frac{f(p_{n-1})}{\frac{f(p_{n-1})-f(p_{n-2})}{p_{n-1}-p_{n-2}}}\\ p_{n-1}-\frac{f(p_{n-1})}{\frac{f(p_{n-1})-f(p_{n-3})}{p_{n-1}-p_{n-3}}}\\ \end{cases}$$

dependent on $f(p_{n-1})f(p_{n-2})<0$ or $f(p_{n-1})f(p_{n-3})<0$.

• Line1 across $(p_{n-1},f(p_{n-1}))$, $(p_{n-2},f(p_{n-2}))$.

• Line2 across $(p_{n-1},f(p_{n-1}))$, $(p_{n-3},f(p_{n-3}))$.

We can choose which line to use by where zero point locates: $[p_{n-2},p_{n-1}]$ or $[p_{n-1},p_{n-3}]$.

# Error Analysis for Iteration

Definition: If positive constants $\lambda$ and $\alpha$ exists with

$$\lim_{n\rightarrow\infty}\frac{|p_{n+1}-p|}{|p_n-p|^{\alpha}}=\lambda$$

then $\{p_n\}^\infty_{n=0}$ converges to $p$ of order $\alpha$, with asymptotic error constant $\lambda$.

Higher the order is, more rapidly the sequence converges.

# Fixed point method

Theorem: If $|g'(x)|\leq k<1$ for any $x\in(a,b)$, then the fixed point iteration converges linearly to the unique fixed point.

• Proof:

$$\lim_{n\rightarrow\infty}\frac{|p_{n+1}-p|}{|p_n-p|^1}=\lim_{n\rightarrow\infty}\frac{|g(p_n)-g(p)|}{|p_n-p|}=\lim_{n\rightarrow\infty}|g'(\xi_n)|=|g'(p)|=Constant\\ \xi_n\in(p_n,p)$$

​ So $\alpha=1$。If $g'(p)=0$, then the order might be higher.

What’s more, if $g'(p)=0$ and $|g''(p)|\leq M$, then we have:

$$\lim_{n\rightarrow\infty}\frac{|p_{n+1}-p|}{|p_n-p|^2}=\lim_{n\rightarrow\infty}\frac{|g(p_n)-g(p)|}{|p_n-p|^2}\\ =\lim_{n\rightarrow\infty}\frac{|g(p)+g'(p)(p_{n}-p)+\frac{g''(\xi)}{2}(p_n-p)^2-g(p)|}{|p_n-p|^2}\\ =\lim_{n\rightarrow\infty}\frac{|g''(\xi)|}{2}=\frac{M}{2}$$

So $\alpha=2$ if the $g''(x)$ has a bound.

---

If $|g'(p)|>1$, then the iteration diverges with any starting point other than p.

# From Root to Fix point

For root finding problem $f(x)=0$, we let $g(x)$ be in the form

$$g(x)=x-\phi(x)f(x),\phi(x)\neq 0$$

so $g(x)=x\Rightarrow f(x)=0$.

Since $g'(x)=1-\phi'(x)f(x)-\phi(x)f'(x)$, let $x=p$, $g'(p)=1-\phi(p)f'(p)$, and $g'(p)=0$ if and only if $\phi(x)=\frac{1}{f'(p)}$,

Definition: A solution $p$ of $f(x)=0$ is a zero of multiplicity $m$ of $f(x)$ if for $x\neq p$, we can write

$$f(x)=(x-p)^mq(x)\\ \lim_{x\rightarrow p}q(x)\neq 0$$

when it comes to Newton’s method, $g(x)=x-f(x)/f'(x)$, we have

$$g(x)=x-\frac{(x-p)^mq(x)}{m(x-p)^{m-1}q(x)+(x-p)^mq'(x)}\\ =x-(x-p)\frac{q(x)}{mq(x)+(x-p)q'(x)}\\ \lim_{x\rightarrow p}g'(x)=1-\frac{1}{m}\neq 0$$

注：这里我存疑，因为都是极限意义下，$x=p$ 时$g(x)$ 甚至无定义。

但这是数值计算，实际中$x$ 和$p$ 总有点差值的。

# Aitken’s Δ2 Method

Suppose $\{p_n\}^\infty_n=0$ is a linearly convergent sequence ($\alpha=1$) with limit $p$:

$$\lim_{n\rightarrow\infty}\frac{p_{n+1}-p}{p_n-p}=\lambda\neq 0$$

So when n is large enough:

$$\frac{p_{n+1}-p}{p_n-p}\approx\frac{p_{n+2}-p}{p_{n+1}-p}\\ p\approx p_n-\frac{(p_{n+1}-p_n)^2}{p_{n+2}-2p_{n+1}+p_n}$$

Aitken’s $\Delta^2$ method is to define a new sequence:

$$q_n=p_n-\frac{(p_{n+1}-p_n)^2}{p_{n+2}-2p_{n+1}+p_n}$$

$\{q_n\}^\infty_{n}$ will converge to $p$ more rapidly than $\{p_n\}$。

Definition: The forward difference $\Delta p_n$ of a sequence $\{p_n\}^\infty_n$ is defined by:

$$\Delta p_n=p_{n+1}-p_n$$

So $q_n=p_n-\frac{(\Delta p_n)^2}{\Delta^2p_n}$.

Theorem: $\lim_{n\rightarrow\infty}\frac{q_n-p}{p_n-p}=0$, so $q_n$ converges faster than $p_n$.

# Steffensen’s Method (Aitken’s Method for fixed point iteration)

Repeat:

$$\begin{cases}p_0\\p_1=g(p_0)\\p_2=g(p_1)\\p_3=p_0-\frac{(\Delta p_0)^2}{\Delta^2p_0}\end{cases} \begin{cases}p_3\\p_4=g(p_3)\\p_5=g(p_4)\\p_6=p_3-\frac{(\Delta p_3)^2}{\Delta^2p_3}\end{cases} \begin{cases}p_6\\p_7=g(p_6)\\p_8=g(p_7)\\p_9=p_6-\frac{(\Delta p_6)^2}{\Delta^2p_6}\end{cases}......$$

Theorem: If $g'(p)\neq 1$, and there exists a $\delta>0$ such that $g\in C^3[p-\delta,p+\delta]$, then Steffensen’s method gives a quadratic ($\alpha=2$) convergence for any $p_0\in[p-\delta,p+\delta]$.

# Zeros of Polynomials and Muller’s Method

Find the root of given polynomial:

$$f(x)=x^n+a_1x^{n-1}+...+a_{n-1}x+a_n$$

initially we have three points: $(x_0,f(x_0)),(x_1,f(x_1)),(x_2,f(x_2))$, then we draw a parabola across the three points:

There are two meeting points of x-axis and the parabola, we choose the one which is closer to $x_2$ as $x_3$, and repeat the iteration to draw parabola across $(x_1,f(x_1)),(x_2,f(x_2)),(x_3,f(x_3))$.

When $|x_{n+1}-x_n|<\varepsilon$, the iteration terminates. and we’ve found the root.

# 第三章 差值和多项式拟合

Theorem (Weierstrass Approximation Theorem): Suppose $f$ is defined and continuous on $[a,b]$, for each $\varepsilon>0$, there exists a polynomial $P(x)$ defined on $[a,b]$, with the property that:

$$\forall x\in[a,b],\;|f(x)-P(x)|<\varepsilon$$

# the n-th Lagrange Interpolating Polynomial

$$L_{n,k}(x)=\prod_{i=0,i\neq k}^n\frac{x-x_i}{x_k-x_i}\\ L_{n,k}(x)=\begin{cases}0&x\neq x_k\\1&x=x_k\end{cases} \\ P(x)=\sum_{k=0}^nf(x_k)L_{n,k}(x)$$

Theorem: Suppose $x_0,x_1,..,x_n$ are distinct numbers in the interval $[a,b]$ and $f\in C^{n+1}[a,b]$, then for each $x\in[a,b]$, a number $\xi(x)\in(a,b)$ exists with

$$f(x)=P_n(x)+\frac{f^{(n+1)}(\xi(x))}{(n+1)!}(x-x_0)(x-x_1)...(x-x_n)$$

Where $P_n(x)$ is the Lagrange Interpolating polynomial.

The error is related to $h=|b-a|$. If $|f^{(n+1)}|\leq M_{n+1}$, then:

$$|f(x)-P_1(x)|=\frac{|f''(\xi)|}{2}|(x-x_0)(x-x_1)|\leq\frac{M_2h^2}{8}\\ |f(x)-P_n(x)|\leq \frac{M_{n+1}}{n!}h^n$$

Runge Phenomenon: when $n$ is bigger, the difference between function and polynomial is larger.

# Data Approximation and Neville’s Method

Theorem: Consider $n$ points: $x_1,x_2,...,x_n$, and Lagrange polynomial:

$$P_1(x)=\sum_{k=0,k\neq j}^nL_{n,k}(x)f(x_k)\\ P_2(x)=\sum_{k=0,k\neq i}^nL_{n,k}(x)f(x_k) \\$$

Then $P(x)=\frac{(x-x_j)P_1(x)-(x-x_i)P_2(x)}{x_i-x_j}$.

So the Lagrange Polynomial can be generate recursively.

$$(x_1,x_2,x_3,...,x_{n-1},x_n)\\ \Rightarrow(x_1,x_2,..,x_{n-1})+(x_2,x_3,...,x_n)\\ \Rightarrow (x_1,x_2,...,x_{n-2})+2(x_2,x_3,...,x_{n-1})+(x_3,x_4,...,x_n)\\ \Rightarrow ...$$

# Divided Difference

Definition: The zeroth divided difference of the function $f$ with respect to $x_i$, denoted $f[x_i]$, is simply the value of $f$ at $x_i$:

$$f[x_i]=f_i=f(x_i)$$

while the first divided difference of $f$ is defined as

$$f[x_i,x_{i+1}]=\frac{f[x_{i+1}]-f[x_i]}{x_{i+1}-x_i}$$

the second divided difference is

$$f[x_i,x_{i+1},x_{i+2}]=\frac{f[x_{i+1},x_{i+2}]-f[x_{i},x_{i+1}]}{x_{i+2}-x_i}$$

Theorem: $P_n(x)$ can be rewritten as

$$P_n(x)=f[x_0]+f[x_0,x_1](x-x_0)+f[x_0,x_1,x_2](x-x_0)(x-x_1)(x-x_2)+...$$

which is known as Newton’s interpolatory divided difference formula.

Theorem: a number $\xi$ exists in $(a,b)$ with

$$f[x_0,x_1,...,x_n]=\frac{f^{(n)}(\xi)}{n!}$$

---

Definition: the forward difference is defined by:

$$\Delta f(x_i)=f(x_{i+1})-f(x_i)$$

if $\forall i,x_{i+1}-x_i\equiv h$ and $x=x_0+sh$, then

$$P_n(x)=f[x_0]+f[x_0,x_1]sh+f[x_0,x_1,x_2]s(s-1)h^2+...\\ =\sum_{k=0}^nf[x_0,x_1,...,x_k]C_s^kk!h^k\\ \because f[x_0,x_1,..,x_k]=\frac{1}{k!h^k}\Delta^2f(x_0)\\ \therefore P_n(x)=\sum_{k=0}^nC_s^k\Delta^kf(x_0)$$

The formula is named as Newton Forward Difference Formula.

Definition: the backward difference is defined by:

$$\bigtriangledown f(x_i)=f(x_i)-f(x_{i-1})$$

So we have

$$f[x_n,x_{n-1}]=\frac{1}{h}\bigtriangledown f(x_n)\\ f[x_n,x_{n-1},x_{n-2},...,x_{n-k}]=\frac{1}{k!h^k}\bigtriangledown^kf(x_n)\\ P_n(x)=\sum_{k=0}^n(-1)^lC_{-s}^k\bigtriangledown^k f(x_n)\\ where\;C_{-s}^k=\frac{-s(-s-1)...(-s-k+1)}{k!}$$

The formula is named as Newton Backward Difference Formula.

# Hermite Interpolation

Problem: Given $n+1$ distinct numbers $x_0,x_1,...,x_n$ and non-negative integers $k_0,k_1,...,k_n$, to construct the osculating polynomial approximating a function $f\in C^m[a,b]$ where $m=max(k_0,k_1,...,k_n)$.

The osculating polynomial $H(x)$ requires:

$$f(x_0)=H(x_0),f'(x_0)=H'(x_0),...,f^{(k_0)}(x_0)=H^{(k_0)}(x_0)\\ f(x_1)=H(x_1),f'(x_1)=H'(x_1),...,f^{(k_1)}(x_1)=H^{(k_1)}(x_1)\\ ......\\ f(x_n)=H(x_n),f'(x_n)=H'(x_n),...,f^{(k_n)}(x_n)=H^{(k_n)}(x_n)\\$$

Obviously the degree of this osculating polynomial is at most $K=\sum_{i=1}^nk_i+n$.

since the number of equation to be satisfied is $K+1$.

---

if $k_0=k_1=...=k_n$, then the Hermite polynomial is:

$$H_{2n+1}(x)=\sum_{j=0}^n f(x_j)H_{n,j}(x)+\sum_{j=0}^nf'(x_j)\hat{H}_{n,j}(x)\\ where\;H_{n,j}(x)=[1-2(x-x_j)L'_{n,j}(x_j)]L^2_{n,j}(x)\\ \hat{H}_{n,j}(x)=(x-x_j)L^2_{n,j}(x)$$

the error is

$$f(x)=H_{2n+1}(x)+\frac{\prod_{k=0}^n(x-x_k)^2}{(2n+2)!}f^{(2n+2)}(\xi)$$

# Piecewise Interpolation Linear Polynomial

Instead of constructing a high-degree polynomial to approximate the function, we can construct several segments:

First construct the base function:

$$l_0(x)=\frac{x-x_1}{x_0-x_1},x_0\leq x\leq x_1\\ l_i(x)=\begin{cases}\frac{x-x_{i-1}}{x_i-x_{i-1}}&x_{i-1}\leq x\leq x_i\\\frac{x-x_{i+1}}{x_i-x_{i+1}}&x_{i}\leq x\leq x_{i+1} \end{cases}\\ l_n(x)=\frac{x-x_{n-1}}{x_n-x_{n-1}},x_{n-1}\leq x\leq x_n$$

Then the piecewise interpolation linear polynomial is:

$$P(x)=\sum_{i=0}^nl_i(x)f(x_i)$$

the error bounds is:

$$|f(x)-P(x)|\leq\frac{h^2}{8}M$$

where $h=max(x_{i+1}-x_i)$, $M=max|f''(x)|$.

---

When considering the osculating polynomial and $k_0=k_1=...=k_n=1$, we can construct:

$$P(x)=\sum_{i=0}^n H_i(x)f(x_i)+\hat{H}_i(x)f'(x_i)$$

where

$$H_i(x)=\begin{cases}(1+2\frac{x-x_i}{x_{i-1}-x_i})(\frac{x-x_{i-1}}{x_i-x_{i-1}})^2&x_{i-1}\leq x\leq x_i\\ (1+2\frac{x-x_i}{x_{i+1}-x_i})(\frac{x-x_{i+1}}{x_i-x_{i+1}})^2&x_{i}\leq x\leq x_{i+1} \end{cases}\\ \hat{H}_i(x)=\begin{cases}(x-x_i)(\frac{x-x_{i-1}}{x_i-x_{i-1}})^2&x_{i-1}\leq x\leq x_i\\ (x-x_i)(\frac{x-x_{i+1}}{x_i-x_{i+1}})^2&x_i\leq x\leq x_{i+1} \end{cases}$$

# 第四章 数值微分和数值积分

Formula:

$$f'(x_0)\approx \frac{f(x_0+h)-f(x_0)}{h}\\ f'(x_0)\approx \frac{f(x_0)-f(x_0-h)}{h}\\ f'(x_0)\approx \frac{f(x_0+h)-f(x_0-h)}{2h}$$

# (n+1)-Point Formula

$$\because f(x)=\sum_{k=0}^nf(x_k)L_{n,k}(x)+\frac{(x-x_0)...(x-x_n)}{(n+1)!}f^{(n+1)}(\xi)\\ \therefore f'(x_j)=\sum_{k=0}^nf(x_k)L'_{n,k}(x_j)+\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{k=0,k\neq j}^n(x_j-x_k)$$

note that $f'(x)$ is only meaningful when $x=x_0,x_1,...,x_n$.

when $n=2$ and $x_1=x_0+h,x_2=x_0+2h$:

$$f'(x_0)=\frac{1}{2h}[-3f(x_0)+4f(x_0+h)-f(x_0+2h)]+\frac{h^2}{3}f^{(3)}(\xi)$$

---

By Taylor polynomial, we have

$$f(x_0+h)=f(x_0)+f'(x_0)h+\frac{1}{2}f''(x_0)h^2+\frac{1}{6}f^{(3)}(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_1)h^4\\ f(x_0-h)=f(x_0)-f'(x_0)h+\frac{1}{2}f''(x_0)h^2-\frac{1}{6}f^{(3)}(x_0)h^3+\frac{1}{24}f^{(4)}(\xi_{-1})h^4\\$$

Therefore:

$$f(x_0+h)+f(x_0-h)=2f(x_0)+f''(x_0)h^2+\frac{h^4}{24}[f^{(4)}(\xi_1)+f^{(4)}(\xi_{-1})]$$

Suppose $f^{(4)}$ is continuous, then there exists $\xi\in[x_0-h,x_0+h]$ satisfying

$$f^{(4)}(\xi)=\frac{1}{2}[f^{(4)}(\xi_1)+f^{(4)}(\xi_{-1})]$$

So

$$f''(x_0)=\frac{1}{h^2}[f(x_0-h)-2f(x_0)+f(x_0+h)]-\frac{h^4}{12}f^{(4)}(\xi)$$

# Richardson’s Extrapolation

Suppose that for each number $h\neq 0$ we have a formula $N(h)$ that approximates an unknown value $M$ and that the truncation error involved with the approximation has the form:

$$M-N(h)=K_1h+K_2h^2+K_3h^3+...$$

Then we have

$$M=N(\frac{h}{2})+K_1\frac{h}{2}+K_2\frac{h^2}{4}+...$$

Therefore

$$M=[N(\frac{h}{2})+N(\frac{h}{2})-N(h)]+K_2(\frac{h^2}{2}-h^2)+...\\ =2N(\frac{h}{2})-N(h)+O(h^2)$$

let $N_2(h)=2N(\frac{h}{2})-N(h)$, then

$$\because M=N_2(h)-\frac{K_2}{2}h^2-\frac{3K_3}{4}h^3-...\\ \therefore3M=4N_2(\frac{h}{2})-N_2(h)+O(h^3)$$

---

In general, if $M$ can be written in the form:

$$M=N_1(h)+\sum_{j=1}^{m-1}K_jh^j+O(h^m)$$

then

$$N_j(h)=N_{j-1}(\frac{h}{2})+\frac{N_{j-1}(h/2)-N_{j-1}(h)}{2^{j-1}-1}\\ M=N_j(h)+O(h^j)$$

# Numerical Integration

$$\because f(x)=\sum_{i=0}^nf(x_i)L_{n,i}(x)+\frac{f^{(n+1)}(\xi)}{(n+1)!}\prod_{i=0}^n(x-x_i)\\ \therefore \int_a^bf(x)dx=\sum_{i=0}^na_if(x_i)+\frac{1}{(n+1)!}\int_a^bf^{(n+1)}(\xi)\prod_{i=0}^n(x-x_i)dx\\ where\;a_i=\int_a^bL_{n,i}(x)dx$$

let $I_n(f)=\sum_{i=0}^na_if(x_i)$, $E_n(f)=\frac{1}{(n+1)!}\int_a^bf^{(n+1)}(\xi)\prod_{i=0}^n(x-x_i)dx$.

• when $n=1$, we have the trapezoidal rule:

$$\int_a^bf(x)dx=\frac{b-a}{2}[f(a)+f(b)]-\frac{h^3}{12}f''(\xi)\\ where\;h=b-a$$

• when $n=2$, we have the Simpson’s rule:

$$\int_a^bf(x)dx=\frac{h}{3}[f(a)+4f(\frac{a+b}{2})+f(b)]-\frac{h^5}{90}f^{(4)(\xi)}\\ where\;h=\frac{b-a}{2}$$

# Degree of Accuracy

Definition: The degree of accuracy (precision) of a quadrature formula is the largest positive integer $n$ such that the formula is precise for $P(x)=x^n$.

For example, Simpson’s rule is accurate for $\int x^3dx$ in that $P^{(4)}(x)\equiv 0$, so the error is zero.

# Newton-Cotes formula

Close Newton-Cotes formula:

$$h=\frac{b-a}{n}\\ x_0=a,x_1=a+h,...,x_n=b\\ N(x)=\begin{cases} \sum_{i=0}^na_if(x_i)+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_{-1}^{n+1}t^2(t-1)...(t-n)dt&n\;is\;even\\ \sum_{i=0}^na_if(x_i)+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_{-1}^{n+1}t(t-1)...(t-n)dt&n\;is\;odd \end{cases}\\ where\;a_i=\int_a^bL_{n,i}(x)dx$$

• when $n=0$ we have the trapezoidal rule.

• when $n=1$ we have the Simpson’s rule.

---

Open Newton-Cotes formula:

$$h=\frac{b-a}{n+2}\\ x_0=a+h,x_1=a+2h,...,x_n=b-h\\ N(x)=\begin{cases} \sum_{i=0}^na_if(x_i)+\frac{h^{n+3}f^{(n+2)}(\xi)}{(n+2)!}\int_{-1}^{n+1}t^2(t-1)...(t-n)dt&n\;is\;even\\ \sum_{i=0}^na_if(x_i)+\frac{h^{n+2}f^{(n+1)}(\xi)}{(n+1)!}\int_{-1}^{n+1}t(t-1)...(t-n)dt&n\;is\;odd \end{cases}\\ where\;a_i=\int_a^bL_{n,i}(x)dx$$

• when $n=0$, we have the Midpoint rule:

$$\int_a^bf(x)dx=2hf(\frac{a+b}{2})+\frac{h^3}{3}f''(\xi)\\ where\;h = (b-a)/2$$

• when $n=1$:

$$\int_a^bf(x)dx=\frac{3h}{2}[f(a+h)+f(a+2h)]+\frac{3h^3}{4}f''(\xi)\\ where\;h=(b-a)/3$$

The difference between open Newton-Cotes formula and close Newton-Cotes formula is the choice of $x_0,...,x_n$. For close case, $x_0=a$ and $x_n=b$, $h=\frac{b-a}{n}$ and $x_i=x_0+ih$. For open case, $x_0=a+h,x_n=b-h,h=\frac{b-a}{n+2}$ and $x_i=x_0+ih$. So the difference is whether $a,b$ are involved.

# Composite Rule

The idea is to divide the interval $[a,b]$ into $n$ subintervals:

$$\int_a^bf(x)dx=\sum_{j=0}^{n-1}\int_{x_j}^{x_{j+1}}f(x)dx\\ \approx\sum_{j=0}^{n-1}\frac{x_{j+1}-x_j}{6}[f(x_j)+4f(\frac{x_j+x_{j+1}}{2})+f(x_{j+1})]$$

The formula above is the composite Simpson’s rule. The composite trapezoidal rule is just similar to Simpson’s rule. And we have composite mid point rule and so on.

# Adaptive Quadrature Method

The basic idea of adaptive quadrature method is to estimate the error of Simpson’s rule, if the error is beyond tolerance, then we divide the interval into two separate subintervals and apply the Simpson’s rule each.

If we view the interval as a whole:

$$\int_a^bf(x)dx=\frac{b-a}{3}[f(a)+4f(\frac{a+b}{2})+f(b))]-\frac{(b-a)^5}{90}f^{(4)}(\xi)\\ \xi\in(a,b)$$

If we divide the interval evenly:

$$\int_a^bf(x)dx=\int_a^{\frac{a+b}{2}}f(x)dx+\int_{\frac{a+b}{2}}^bf(x)dx\\ =\frac{b-a}{6}[f(a)+4f(\frac{3a+b}{4})+f(\frac{a+b}{2})]\\ +\frac{b-a}{6}[f(\frac{a+b}{2})+4f(\frac{a+3b}{4})+f(b)]\\ -\frac{(b-a)^5}{90*2^5}[f^{(4)}(\xi_1)+f^{(4)}(\xi_2)]\\ \xi_1\in(a,\frac{a+b}{2}),\xi_2\in(\frac{a+b}{2},b)$$

---

$$\because \frac{f^{(4)}(\xi_1)+f^{(4)}(\xi_2)}{2}=f^{(4)}(\hat{\xi})\\ \hat\xi\in(\xi_1,\xi_2)\in(a,b)\\ \therefore \int_a^bf(x)dx=\frac{b-a}{6}[f(a)+4f(\frac{3a+b}{4})+2f(\frac{a+b}{2})+4f(\frac{a+3b}{4})+f(b)]\\ -\frac{(b-a)^5}{90*2^4}f^{(4)}(\hat{\xi})$$

The error estimation is derived by assuming that $\xi\approx\hat{\xi}$, or $f^{(4)}(\xi)\approx f^{(4)}(\hat{\xi})$.

The success of the technique depends on the accuracy of this assumption! (If $b-a$ is small enough, then $\xi\approx\hat{\xi}$)

So we have the estimation of error:

$$\frac{b-a}{6}[f(a)+4f(\frac{3a+b}{4})+2f(\frac{a+b}{2})+4f(\frac{a+3b}{4})+f(b)]\\ -\frac{(b-a)^5}{90*2^4}f^{(4)}(\hat{\xi})\approx \frac{b-a}{3}[f(a)+4f(\frac{a+b}{2})+f(b))]-\frac{(b-a)^5}{90}f^{(4)}(\xi)\\ \Rightarrow \frac{(b-a)^5}{90}f^{(4)}(\xi)\approx \frac{16}{15}[...]$$

The method is to divide the interval if the estimation of error is beyond tolerance.

# Gaussian Quadrature

Definition: Legendre polynomials are a collection \

• for each $n$, $P_n(x)$ is a polynomial of degree $n$.
• $\int_{-1}^1P(x)P_n(x)dx=0$ if $P(x)$ is a polynomial of degree $< n$.

The first few Legendre polynomials are

$$P_0=1\\ P_1(x)=x\\ P_2(x)=x^2-\frac{1}{3}\\ P_3(x)=x^3-\frac{3}{5}x\\ ...$$

How to derive Legendre polynomial is quite simple:

$$P_n(x)=x^n+a_1x^{n-1}+...+a_{n-1}x+a_n\\ \begin{cases}\int_{-1}^1P_n(x)=0\\ \int_{-1}^1xP_n(x)=0\\ \int_{-1}^1x^2P_n(x)=0\\ ...\\ \int_{-1}^1x^{n-1}P_n(x)=0\\ \end{cases}$$

fine, that seems not that simple. We need to solve the equation system.

Theorem: Suppose that $x_1,x_2,..,x_n$ are the roots of the $nth$ Legendre polynomial $P_n(x)$, and that for each $i=1,2,...,n$, the numbers $c_i$ are defined by

$$c_i=\int_{-1}^1\prod_{j=1,j\neq i}^n\frac{x-x_j}{x_i-x_j}dx$$

If $P(x)$ is any polynomial of degree $<2n$, then

$$\int_{-1}^1P(x)dx=\sum_{i=1}^nc_iP(x_i)$$

# 第五章 线性常系数微分方程

Problem：$\frac{dy}{dt}=f(t,y),a\leq t\leq b$. subject to $y(a)=\alpha$.

Definition: A function $f(t,y)$ is said to satisfy a Lipschitz Condition in the variable $y$ on a set $D\subset\mathbb{R}^2$, if a constant $L>0$ exists with the property that

$$|f(t,y_1)-f(t,y_2)|\leq L|y_1-y_2|$$

whenever $(t,y_1),(t,y_2)\in D$. The constant $L$ is called a Lipschitz Constant for $f$.

Definition: A set $D\subset\mathbb{R}^2$ is said to be convex if whenever

$$(t_1,y_1),(t_2,y_2)\in D$$

then for any $\lambda\in[0,1]$, $((1-\lambda)t_1+\lambda t_2,(1-\lambda)y_1,\lambda y_2)\in D$.

Theorem: Suppose $f(t,y)$ is defined on a convex set $D\subset\mathbb{R}^2$, If a constant $L$ exists with

$$|\frac{\partial f}{\partial y}(t,y)|\leq L$$

for all $(t,y)\in D$, Then $f$ satisfies a Lipschitz condition on $D$ in the variable $y$ with Lipschitz constant $L$.

Definition: The initial problem

$$\begin{cases}y'(t)=f(t,y)\\y(a)=\alpha\end{cases}$$

is a well-posed problem if

• a unique solution $y(t)$ exists,

• for any $\varepsilon>0$, there exists a positive constant $k(\varepsilon)$ with the property that, whenever $|\varepsilon_0|<\varepsilon$ and $\delta(t)$ is continuous with $|\delta(t)|<\varepsilon$ on $[a,b]$, a unique solution $z(t)$ to the problem

  $$\begin{cases}z'(t)=f(t,z)+\delta(t)\\z(a)=\alpha+\varepsilon_0\end{cases}$$

  exists with $|z(t)-y(t)|<k(\varepsilon)\varepsilon$.

Theorem: Suppose that $D=\{(t,y)|a\leq t\leq b,-\infty\leq y\leq\infty\}$ and $f(t,y)$ is continuous on $D$. If $f$ satisfies a Lipschitz condition on $D$ in the variable $y$, then problem:

$$\begin{cases}y'(t)=f(t,y)\\y(a)=\alpha\end{cases}$$

has a unique solution $y(t)$ and is well-posed.

# Euler’s Method

To solve a well-posed initial-value problem:

$$\begin{cases}y'(t)=f(t,y)\\y(a)=\alpha\end{cases}$$

First, we construct the mesh points in the interval $[a,b]$ with equal space:

$$a=t_0<t_1<...<t_N=b$$

with equal stepi size $h=(b-a)/N$.

Due to Taylor formula, we have:

$$y(t_{j+1})=y(t_j)+hy'(t_j)+\frac{h^2}{2}y''(\xi_j)\\ =y(t_j)+hf(t_j,y(t_j))+\frac{h^2}{2}y''(\xi_j)$$

So we can generate $y(t_1),y(t_2),...,y(t_N)$ based on $y(t_0)=\alpha$.

Then we can construct Lagrange polynomial with $(t_0,y(t_0)),...,(t_N,y(t_N))$.

Lemma: If $s$ and $t$ are positive real number, $\{a_j\}_{j=0}^k$ is a sequence satisfying $a_0\geq-\frac{t}{s}$ and

$$a_{j+1}\leq (1+s)a)j+t,j=0,1,...,k$$

then

$$a_{j+1}\leq e^{(j+1)s}(a_0+\frac{t}{s})-\frac{t}{s}$$

Theorem: Suppose $f$ is continuous and satisfies a Lipschitz Condition on $D$, and that $|y''(t)|\leq M$. $y^*(t)$ is the solution for the problem, and $y(t_j)$ is generated by Euler’s Method, then:

$$|y^*(t_j)-y(t_j)|\leq\frac{hM}{2L}[e^{L(t_j-a)}-1]$$

---

# Improved Euler Method (Predict-Correct)

we can generate $y(t_0),...,y(t_n)$ using Euler’s Method, and correct them using trapezoidal rule:

$$\frac{y(t_{i+1})-y(t_i)}{h}\approx\frac{y'(t_{i+1})+y'(t_i)}{2}\\ \Rightarrow \hat{y}(t_{i+1})=\frac{h}{2}[f(t_{i+1},y(t_{i+1}))+f(t_i,\hat{y}(t_i)]+\hat{y}(t_i)$$

where $y(t_{i+1})$ is generated by Euler’s Method, and $\hat{y}(t_{i+1})$ is the corrected value.

# Higher-Order Taylor Methods

Definition: The local truncation error in Euler’s Method is defined by

$$\tau_{i+1}(h)=\frac{y(t_{i+1})-(y(t_i)+hf(t_i,y_i))}{h}\\ =\frac{h}{2}y''(\xi_j)$$

Taylor Method of order $n$:

$$y(t_0)=\alpha\\ y(t_{i+1})=y(t_{i})+hT^{(n)}(t_i,y(t_i))\\ where\;T^{(n)}(t_i,y(t_i))=f(t_i,y(t_i))+\frac{h}{2}f'(t_i,y(t_i))+...+\frac{h^{(n-1)}}{n!}f^{(n-1)}(t_i,y(t_i))$$

Euler’s Method is the Taylor method of order 1.

# Runge-Kutta Method

$$\begin{cases}y'(t)=f(t,y)\\y(a)=\alpha\end{cases}\\ a=t_0<t_1<...<t_N=b$$

The Runge-Kutta Method is by computing:

$$y(t_{n+1})=y(t_n)+\sum_{i=1}^Nc_iK_i$$

where

$$K_1=hf(t_n,y_n)\\ K_i=hf(t_n+\alpha_ih,y_n+\sum_{j=1}^{i-1}b_{ij}K_j)$$

The coefficients are generated by Taylor’s Theorem and error analysis.

when $N=4$:

$$y_0=y(a)=\alpha\\ K_1=hf(t_i,y_i)\\ K_2=hf(t_i+\frac{h}{2},y_i+\frac{1}{2}K_1)\\ K_3=hf(t_i+\frac{h}{2},y_i+\frac{1}{2}K_2)\\ K_4=hf(t_{i+1},y_i+K_3)\\ y(t_{i+1})=y(t_i)+\frac{1}{6}(K_1+2K_2+2K_3+K_4)$$

# Adams Fourth-Order Predictor-Corrector

First we need to compute $y(t_0),y(t_1),y(t_2),y(t_3)$ by Runge-Kutta method.

Then we can combine the implicit and the explicit method:

• Use Four-step Adams-Bashforth method as predictor:

  $$\hat{y}(t_{n+1})=y(t_n)+\frac{h}{24}[55f(t_n,y(t_{n}))\\ -59f(t_{n-1},y(t_{n-1}))\\ +37f(t_{n-2},y(t_{n-2}))\\ -9f(t_{n-3},y(t_{n-3}))]$$

  this is predictor for $y(t_{n+1})$ using explicit method.

• Then use Three-step Adams-Moulton Method as a corrector:

  $$y(t_{n+1})=y(t_n)+\frac{h}{24}[ 9f(t_{n+1},\hat{y}(t_{n+1}))\\ +19f(t_n,y(t_n))\\ -5f(t_{n-1},y(t_{n-1}))\\ +f(t_{n-2},y(t_{n-2})) ]$$

  this is corrector for $y(t_{n+1})$ using implicit method while the $y(t_{n+1})$ on the right side of the equation is replaced by the predictor $\hat{y}(t_{n+1})$.

# 第六章 线性方程组

The Linear System of Equations （LSEs）

# Gaussian Eliminiation

• Maximal Row Pivoting Technique:

  choose $max|a_{jk}|,k\leq j\leq n$ as the pivot element.

  choose the maximum integer on the k-th column.

• Partial Pivoting Technique:

  choose $max|a_{ij}|,k\leq i,j\leq n$ as the pivot element. which means that the pivot element is not necessarily on the k-th column.

• Scaled Partial Pivoting Technique:

  compute the max integer for each row:

  $$s_i=max_{1\leq j\leq n}|a_{ij}|$$

  Then choose the number p satisifying:

  $$\frac{a_{p1}}{s_p}=max_{1\leq i\leq n}\frac{a_{i1}}{s_i}$$

  then $a_{p1}$ as the pivot element.

# Matrix Factorization

$$\textbf{Ax=b}\\ A=\textbf{LU}$$

where

$$\textbf{L}=\begin{pmatrix}1&0&0&...&0\\ l_{21}&1&0&...&0\\ l_{31}&l_{32}&1&...&0\\ ...&...&...&...&...\\ l_{n1}&l_{n2}&l_{n3}&...&1 \end{pmatrix}\\ \textbf{U}=\begin{pmatrix} u_{11}&u_{12}&u_{13}&...&u_{1n}\\ 0&u_{22}&u_{23}&...&u_{2n}\\ 0&0&u_{33}&...&u_{3n}\\ ...&...&...&...&...\\ 0&0&0&...&u_{nn} \end{pmatrix}\\$$

Then the LSEs can be written as $\textbf{LUx=b}$. and we can solve:

$$\begin{cases} \textbf{Ly=b}\\ \textbf{Ux=y} \end{cases}$$

The process of factorization is essentially the same as Gaussian Elimination.

Consider Each step of Gaussian Elimination:

• swap row i and row j
• multiple row i by x then add to row j.

The two operations above can be described as matrix. e.g:

$$M=\begin{pmatrix} 0&1&0&0\\ 1&0&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{pmatrix}$$

then $B=MA$ is the result of swapping the first two rows of $A$.

$$M=\begin{pmatrix} 1&0&0&0\\ 2&1&0&0\\ 0&0&1&0\\ 0&0&0&1 \end{pmatrix}$$

then $B=MA$ is the result of adding 2*row1 to row2.

So we have:

$$M_{m}M_{m-1}...M_1A=U$$

which is actually Gaussian Elimination.

So $A=LU=M_1^{-1}...M_m^{-1}U$.

# Doolittle Factorization

The Method is to apply the LU factorization directly.

First select

$$\begin{cases} u_{1j}=a_{1j}&1\leq j\leq n\\ l_{i1}=a_{i1}/u_{11}&1\leq i\leq n \end{cases}$$

If $a_{11}=0$, then the factorization is impossible. (Not unique solution)

Then compute

$$\begin{cases} u_{ij}=a_{ij}-\sum_{k=1}^{i-1}l_{ik}u_{kj}&j\geq i\\ l_{ij}=(a_{ij}-\sum_{k=1}^{j-1}l_{ik}u_{kj})/u_{jj}&i\geq j+1 \end{cases}$$

Doolittle factorization works only when we don’t need any row swaps in Gaussian Elimination.

# Special types of matrices

Definition: The $n\times n$ matrix $A$ is said to be strictly diagonally dominant when

$$\forall i\in[1,n].\quad|a_{ii}|>\sum_{j=1,j\neq i}^n|a_{ij}|$$

Theorem: A strictly diagonally dominant matrix $A$ is nonsingular.

In this case, Gaussian elimination can be performed on any linear system of the form $\textbf{Ax}=\textbf{b}$ to obtain its unique solution without row or column interchanges.

Definition: A matrix $A$ is positive definite if it is symmetric and for every column vector $x$ holds $x^TAx>0$.

Theorem: If $A$ is an $n\times n$ positive definite matrix, then

• $A$ is nonsingular
• $a_{ii}>0$ for all $1\leq i\leq n$.
• a_{ij}^2

Theorem: A symmetric matrix $A$ is positive definite if and only if each of its leading principal submatrices has a positive determinant.

Corollary: The matrix $A$ is positive definite if and only if $A$ can be factored in the form $LDL^T$, where $L$ is lower triangular with 1s on its diagonal and $D$ is a diagonal matrix with positive diagonal entries.

Corollary: The matrix $A$ is positive definite if and only if $A$ can be factored in the form $LL^T$, where $L$ is lower triangular with nonzero diagonal entries.

# Choleski’s Method

The method is to factor a positive definite matrix $A=LL^T$.

First set $l_{11}=\sqrt{a_{11}}$.

For $ j =2,…,n$, set $l_{j1}=a_{1j}/l_{11}$.

Then for $i=2,...,n-1$, $j\geq i+1$:

$$l_{ii}=\sqrt{a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2}\\ l_{ji}=(a_{ij}-\sum_{k=1}^{i-1}l_{ik}l_{jk})/l_{ii}$$

Finally set $l_{nn}=\sqrt{a_{nn}-\sum_{k=1}^{n-1}l_{nk}^2}$.

---

If we need to factor $A=LDL^T$：

First set $d_{11}=a_{11}$.

For $j=2,...,n$, set $l_{j1}=a_{1j}/d_{11}$.

For $i=2,...,n$, $j\geq i+1$

$$d_{ii}=a_{ii}-\sum_{j=1}^{i-1}l_{ij}^2d_{jj}\\ l_{ji}=(a_{ij}-\sum_{k=1}^{i-1}l_{ik}l_{jk}d_{kk})/d_{ii}$$

# 第七章 线性方程组的迭代

# Norms of Vectors and Matrices

Definition: A vector norm on $\mathbb{R}^n$ is a function: $||\cdot||$, from $\mathbb{R}^n$ into $\mathbb{R}$ with the following properties:

• $||\textbf{x}||\geq 0$ for all $\textbf{x}\in\mathbb{R}^n$.
• $||\textbf{x}||=0$ if and only if $\textbf{x}=(0,0,...,0)^T$.
• $||\alpha\textbf{x}||=|\alpha|\cdot||\textbf{x}||$ for all $\alpha\in\mathbb{R}$ and $\textbf{x}\in\mathbb{R}^n$.
• $||\textbf{x+y}||\leq ||\textbf{x}||+||\textbf{y}||$ for all $\textbf{x,y}\in\mathbb{R}^n$.

Definition: The $l_2$ and $l_{\infty}$ norm for vectors are defined by

$$||\textbf{x}||_2=\sqrt{\sum_{i=1}^n|x_i|^2}\\ ||\textbf{x}||_{\infty}=\max_{1\leq i\leq n}|x_i|$$

$l_2$ norm is called Euclidean Norm.

Definition: the $l_2$ and $l_\infty$ distances between $\textbf{x}$ and $\textbf{y}$ are defined by

$$||\textbf{x}-\textbf{y}||_2=\sqrt{\sum_{i=1}^n|x_i-y_i|^2}\\ ||\textbf{x}-\textbf{y}||_\infty=\max_{1\leq i\leq n}|x_i-y_i|$$

Definition: A sequence $\{\textbf{x}\}_0^\infty$ of vectors is said to converge to $\textbf{x}$ with respect to the norm $||\cdot||$ if :

$$\lim_{n\rightarrow \infty}||\textbf{x}_n-\textbf{x}||=0$$

Theorem: $||\textbf{x}||\leq ||\textbf{x}||_2\leq\sqrt{n}||x||_\infty$

---

Definition: A matrix norm on the set of all $n\times n$ matrices is a real valued function satisfying:

• $||\textbf{A}||\geq 0$
• $||\textbf{A}||=0$ if and only if $\textbf{A}$ is all zero
• $||\alpha\textbf{A}||=\alpha||\textbf{A}||$
• $||\textbf{A}+\textbf{B}||\leq||\textbf{A}||+||\textbf{B}||$
• $||\textbf{AB}||\leq||\textbf{A}||\cdot||\textbf{B}||$

An example of the matrix norm is:

$$||\textbf{A}||=\max_{||\textbf{x}||=1}||\textbf{Ax}||$$

This is called the natural matrix norm.

Corollary: $||\textbf{Ax}||\leq||\textbf{A}||\cdot||\textbf{x}||$. when the norm is natural matrix norm.

Definition: The $l_2$ and $l_{\infty}$ norm for matrices are defined by

$$||\textbf{A}||_\infty=\max_{||\textbf{x}||_\infty=1}||\textbf{Ax}||_\infty\\ ||\textbf{A}||_2=\max_{||\textbf{x}||_2=1}||\textbf{Ax}||_2$$

Theorem: $||\textbf{A}||_\infty=\max_{1\leq i\leq n}\sum_{j=1}^n|a_{ij}|$.

# Eigenvalues and Eigenvectors

Definition: If $\textbf{A}$ is a square matrix, the characteristic polynomial of $\textbf{A}$ is defined by:

$$p(\textbf{A})=det(\textbf{A}-\lambda\textbf{I})$$

Definition: The spectral radius $\rho(\textbf{A})$ of a matrix $\textbf{A}$ is defined by

$$\rho(\textbf{A})=max|\lambda|$$

where $\lambda$ is an eigenvalue of $\textbf{A}$.

Theorem: $||\textbf{A}||_2=\sqrt{\rho(A^TA)}$ and $\rho(\textbf{A})\leq||\textbf{A}||$ for any natural norm.

Definition: we say a matrix $A$ is convergent if

$$\lim_{n\rightarrow\infty}A^n=0$$

Theorem: The following statements are equivalent

• $\textbf{A}$ is a convergent matrix
• $\lim_{n\rightarrow\infty}||\textbf{A}^n||=0$ for some natural norm.
• $\lim_{n\rightarrow\infty}||\textbf{A}^n||=0$ for all natural norms
• $\rho(\textbf{A})<1$
• $\lim_{n\rightarrow\infty}\textbf{A}^n\textbf{x}=0$ for every $\textbf{x}$.

# Iterative techniques for solving linear system

First we can split matrix $A$ into three parts:

• diagonal matrix \textbf
• strictly lower-triangular part \textbf
• strictly upper-triangular part $\textbf{U}$.

So

$$\textbf{Ax}=\textbf{b}\Leftrightarrow \textbf{Dx}=(\textbf{L}+\textbf{U})\textbf{x}+\textbf{b}\\ \therefore \textbf{x}_{n+1}=\textbf{D}^{-1}(\textbf{L}+\textbf{U})\textbf{x}_n+\textbf{D}^{-1}\textbf{b}$$

And this is the matrix form of the Jacobi iterative technique. We can repeat this iteration.

---

Gauss-Seidel method:

$$\textbf{x}_{n+1}=(\textbf{D}-\textbf{L})^{-1}\textbf{U}\textbf{x}_n+(\textbf{D}-\textbf{L})^{-1}\textbf{b}$$

---

Lemma: If the spectral radius $\rho(\textbf{T})$ satisfies $\rho(\textbf{T})<1$, then

$$(\textbf{I}-\textbf{T})^{-1}=\textbf{I}+\textbf{T}+\textbf{T}^2+...$$

Theorem: Iteration $\textbf{x}_{n+1}=\textbf{Tx}_n+\textbf{c}$ converges to the unique solution of $\textbf{x}=\textbf{Tx}+\textbf{c}$ if and only if $\rho(\textbf{T})<1$.

Corollary: If $||\textbf{T}||<1$ for any one natural matrix norm and $\textbf{c}$ is a given vector, then the iteration $\textbf{x}_{n+1}=\textbf{Tx}_n+\textbf{c}$ converges, and the following error bound hold:

• $||\textbf{x}_n-\textbf{x}||\leq ||\textbf{T}||^n||\textbf{x}_0-\textbf{x}||$

# SOR Technique

$$\omega\textbf{Ax}=\omega\textbf{b}\\ \Rightarrow(\textbf{D}-\omega\textbf{L})\textbf{x}_{n+1}=[(1-\omega)\textbf{D}+\omega\textbf{U}]\textbf{x}_n+\omega\textbf{b}\\ \textbf{T}_\omega=(\textbf{D}-\omega\textbf{L})^{-1}[(1-\omega)\textbf{D}+\omega\textbf{U}]\\ \textbf{c}_\omega=\omega(\textbf{D}-\omega\textbf{L})^{-1}\textbf{b}\\ \therefore \textbf{x}_{n+1}=\textbf{T}_\omega\textbf{x}_n+\textbf{c}_\omega$$

Theorem(Kahan): $\rho(\textbf{T}_\omega)\geq|\omega-1|$.

So the SOR method can converge if and only if $0<\omega<2$.

the optimal choice of $\omega$ for the SOR method is

$$\omega=\frac{2}{1+\sqrt{1-[\rho(\textbf{D}^{-1}(\textbf{L}+\textbf{U}))]^2}}$$

and we have $\rho(\textbf{T}_\omega)=\omega-1$.

# 第八章 近似理论

# Discrete Least Squares Approximation

• Minimax Rule:

  $$E_1=\min_{a_0,a_0}\max_{1\leq i\leq 10}|y_i-(a_1x_i+x_0)|$$

• Absolute Deviation Rule:

  $$E_2=\min_{a_0,a_1}\sum_{i=1}^{10}|y_i-(a_1x_i+a_0)|$$

• Least Square Rule:

  $$E = \min_{a_0,a_1}\sum_{i=1}^{10}[y_i-(a_1x_i+a_0)]^2$$

  we have

  $$\begin{cases}\frac{\partial}{\partial a_0}E=\sum_{i=1}^m2(y_i-a_1x_i-a_0)(-1)=0\\ \frac{\partial}{\partial a_1}=\sum_{i=1}^m2(y_i-a_1x_i-a_0)(-x_i)=0\end{cases}$$

---

​ The genera form of discrete least square rule:

$$\min_{a_0,a_1,...,a_n}E=\sum_{i=1}^m(y_i-\sum_{k=0}^na_kx_i^k)^2\\ \frac{\partial E}{\partial a_j}=0\Rightarrow\\ \sum_{k=0}^na_k\sum_{i=1}^m x_i^{j+k}=\sum_{i=1}^my_ix_i^j,j=0,1,...,n$$

​ Let

$$\textbf{R}=\begin{pmatrix}x_1^0&x_1^1&...&x_1^{n-1}&x_2^n\\ x_2^0&x_2^{1}&...&x_2^{n-1}&x_2^n\\ ...&...&...&...&...\\ x_m^0&x_m^{1}&...&x_m^{n-1}&x_m^n \end{pmatrix}, \textbf{a}=\begin{pmatrix}a_0\\a_1\\...\\a_n\end{pmatrix}, \textbf{y}=\begin{pmatrix}y_0\\y_1\\...\\y_m\end{pmatrix}$$

​ Then the equations can be written as

$$\textbf{R}^T\textbf{Ra}=\textbf{R}^T\textbf{y}$$

---

If we can compute the function, then the continuous least square rule is:

$$\min_{a_0,...,a_n}E=\int_a^b[f(x)-P_n(x)]^2dx=\int_{a}^b[f(x)-\sum_{k=0}^na_kx^k]^2dx$$

Definition: The set of functions $\{\phi_0,\phi_1,...,\phi_n\}$ is said to be linearly independent on $[a,b]$ if for all $x\in [a,b]$:

$$c_0\phi_0(x)+c_1\phi_1(x)+...+c_n\phi_n(x)=0\\ \Rightarrow c_0=c_1=...=c_n=0$$

Theorem: If $\phi_j(x)$ is a polynomial of degree j, then $\{\phi_0,...,\phi_n\}$ is linearly independent on any interval.

# Weight function and Orthogonal Set of functions

In addition to the discrete least square rule, the error can be

$$E=\int_a^bw(x)[f(x)-\sum_{k=0}^na_k\phi_k(x)]^2dx$$

where $w(x)$ is the weight function. The normal function can be written as

$$0=\frac{\partial E}{\partial a_j}\Rightarrow\\ \sum_{k=0}^na_k\int_a^bw(x)\phi_k(x)\phi_j(x)dx=\int_a^bw(x)f(x)\phi_j(x)dx,j=0,1,...,n$$

we can choose the function set $\{\phi_0,...,\phi_n\}$ satisfying that:

$$\int_a^bw(x)\phi_k(x)\phi_j(x)dx=\begin{cases}0&j\neq k\\ \alpha_j>0&j=k\end{cases}$$

so the coefficients can be easily solved:

$$a_j=\frac{1}{\alpha_j}\int_a^bw(x)f(x)\phi_j(x)dx$$

Definition: $\{\phi_0,..,\phi_n\}$ is said to be an orthogonal set of functions with respect to the weight function $w(x)$ if

$$\int_a^bw(x)\phi_j(x)\phi_k(x)dx=\begin{cases}0&j\neq k \\ \alpha_k>0&j=k\end{cases}$$

Theorem(Gram-Schmidt orthogonalize):

$$B_k=\frac{\int_a^bxw(x)[\phi_{k-1}(x)]^2dx}{\int_a^bw(x)[\phi_{k-1}(x)]^2dx},C_k=\frac{\int_a^bxw(x)\phi_{k-1}(x)\phi_{k-2}(x)dx}{\int_a^bw(x)[\phi_{k-2}(x)]^2dx}\\ \begin{cases}\phi_0(x)=1\\ \phi_1(x)=x-B_1 \\ ...\\ \phi_k(x)=(x-B_k)\phi_{k-1}(x)-C_k\phi_{k-2}(x) \end{cases}$$

Then we generate a set of polynomial functions.

when $w(x)=1$, then we can generate Legendre Polynomial.

# Chebyshev Polynomials

$$T_n(x)=cos[n\;arccos (x)],x\in[-1,1]\\ T_0=1\\ T_n(\theta(x))=cos(n\theta),\theta=arccos(x)\\ T_{n+1}(\theta)=cos[(n+1)\theta]=2cos(n\theta)cos\theta-T_{n-1}(\theta)\\ T_0=1,T_1=x,\\ T_{n+1}(x)=2xT_n(x)-T_{n-1}(x)$$

Theorem: The Chebyshev polynomials is orthogonal with respect to the weight function

$$w(x)=\frac{1}{\sqrt{1-x^2}}$$

Theorem: Zero points:

$$\bar{x}_k=cos(\frac{2k-1}{2n}\pi),k=1,2,...,n\\ T_n(\bar{x}_k)=0$$

extrema points:

$$\bar{x}_k'=cos(\frac{k\pi}{n}),k=0,1,...,n\\ T_n(\bar{x}_k')=(-1)^k,T_n'(\bar{x}_k')=0$$

---

What’s more, we can normalize Chebyshev Polynomial too ones with leading coefficient 1:

$$\tilde{T}_n(x)=\frac{1}{2^{n-1}}T_n(x)$$

let $\tilde{\prod}_n$ denotes all degree=n polynomials with leading coefficient 1, then:

Theorem:

$$\forall P_n(x)\in\tilde{\prod}_n,\frac{1}{2^{n-1}}=\max_{x\in[-1,1]}|\tilde{T}_n(x)|\leq\max_{x\in[-1,1]}|P_n(x)|$$

equality can occur only if $P_n\equiv\tilde{T}_n$.

Corollary: If $P(x)$ is the interpolating polynomial of degree at most $n$ with nodes at the roots of $T_{n+1}(x)$, then

$$\max_{x\in[-1,1]}|f(x)-P(x)|\leq\frac{1}{2^n(n+1)!}\max_{x\in[-1,1]}|f^{(n+1)}(x)|$$

# 第九章 特征值与特征向量

Definition: left eigenvector is defined by

$$\textbf{y}^T\textbf{A}=\lambda\textbf{y}^T$$

Definition: Spectrum of $\textbf{A}$ is the set of eigenvalues of $\textbf{A}$, denoted by $\lambda(\textbf{A})$.

Definition: Spectrum radius of $\textbf{A}$ is:

$$\rho(\textbf{A})=max\{|\lambda|,\lambda\in\lambda(\textbf{A})\}$$

Theorem: If $\lambda_1,...,\lambda_k$ are distinct eigenvalues with associated eigenvectors

$$\textbf{x}^{(1)},...,\textbf{x}^{(k)}$$

then the set of eigenvectors is linearly independent.

Definition: orthogonal requires $\textbf{v}_i\cdot\textbf{v}_j=c\delta_{ij}$, when $c=1$, it’s orthonormal.

---

Theorem (Gerschgorin Circle Theorem): $\textbf{A}$ is an $n\times n$ matrix, and $\mathbb{R}_i$ denotes the circle in the complex plane with center $a_{ii}$:

$$\mathbb{R}_i=\{z\in\mathcal{C}||z-a_{ii}|\leq\sum_{j=1,j\neq i}^n|a_{ij}|\}$$

Then the eigenvalues of $\textbf{A}$ are contained with $\mathbb{R}=\cup_{i=1}^n\mathbb{R}_i$.

Moreover, the union of any $k$ of these circles that do not intersect the remaining $(n-k)$ circles contains precisely $k$ (counting multiplicities) of the eigenvalues.

# Computing Eigenvalues and Eigenvectors

Iterative Power Method

Assume that the matrix $\textbf{A}$ has n eigenvalues

$$|\lambda_1|\geq|\lambda_2|\geq|\lambda_3|\geq...\geq|\lambda_n|$$

with an associated collection of linearly independent eigenvectors

$$\{\textbf{v}^{(1)},...,\textbf{v}^{(n)}\}$$

If $\textbf{x}$ is any vector in $\mathbb{R}^n$, then constants $\beta_1,...,\beta_n$ exist with

$$\textbf{x}=\sum_{j=1}^n\beta_j\textbf{v}^{(j)}$$

So we have

$$\textbf{A}^k\textbf{x}=\lambda_1^k\sum_{j=1}^n\beta_j(\frac{\lambda_j}{\lambda_1})^k\textbf{v}^{(j)}\\ =\lambda_1^k(\beta_1\textbf{v}^{(1)}+\sum_{j=2}^n\beta_j(\frac{\lambda_j}{\lambda_1})^k\textbf{v}^{(j)})$$

$$\because \lim_{k\rightarrow\infty}(\frac{\lambda_j}{\lambda_1})^k=0\\ \therefore\lim_{k\rightarrow\infty}\textbf{A}^k\textbf{x}=\lim_{k\rightarrow\infty}\lambda_1^k\beta_1\textbf{v}^{(1)}$$

However, the limit may converge to zero or diverge.

# The power method

• Choose an arbitrary unit vector $\textbf{x}^{(0)}$ relative to $||\cdot||_{\infty}$.

• Suppose a component $x_{p0}^{(0)}$ of $\textbf{x}^{(0)}$ with

  $$x_{p0}^{(0)}=1=||\textbf{x}^{(0)}||_\infty$$

• let $\textbf{y}^{(1)}=\textbf{Ax}^{(0)}$, and define

  $$\mu^{(1)}=y_{p0}^{(1)}=\frac{y_{p0}^{(1)}}{x_{p0}^{(0)}} =\frac{\beta_1\lambda_1v_{p0}^{(1)}+\sum_{j=2}^n\beta_j\lambda_jv_{p0}^{(j)}}{\beta_1v_{p0}^{(1)}+\sum_{j=2}^n\beta_jv_{p0}^{(j)}}\\ =\lambda_1[\frac{\beta_1v_{p0}^{(1)}+\sum_{j=2}^n\beta_j(\lambda_j/\lambda_1)v_{p0}^{(j)}}{\beta_1v_{p0}^{(1)}+\sum_{j=2}^n\beta_jv_{p0}^{(j)}}]$$

  where $\textbf{v}^{(i)}$ is the eigenvectors and $\beta_j$ is coefficient satisfying:

  $$\textbf{x}^{(0)}=\sum_{j=1}^n\beta_j\textbf{v}^{(j)}$$

• let $p1$ be the least integer such that $y_{p1}^{(1)}=||\textbf{y}^{(1)}||_\infty$.

  Define $\textbf{x}^{(1)}$ by

  $$\textbf{x}^{(1)}=\frac{\textbf{Ax}^{(0)}}{y_{p1}^{(1)}}$$

• Similarly, define

  $$\textbf{y}^{(m)}=\textbf{Ax}^{(m-1)}\\ \mu^{(m)}=y_{pm}^{(m)}=\lambda_1[\frac{\beta_1v_{p_{m-1}}^{(1)}+\sum_{j=2}^n\beta_j(\lambda_j/\lambda_1)^mv_{p_{m-1}}^{(j)}}{\beta_1v_{p_{m-1}}^{(1)}+\sum_{j=2}^n\beta_j(\lambda_j/\lambda_1)^{m-1}v_{p_{m-1}}^{(j)}}]$$

  let $p_m$ be the smallest integer with $|y_{pm}^{(m)}|=||\textbf{y}^{(m)}||_\infty$,

  $$\textbf{x}^{(m)}=\frac{\textbf{y}^{(m)}}{y_{pm}^{(m)}}=\frac{\textbf{A}^m\textbf{x}^{(0)}}{\prod_{k=1}^my_{pk}^{(k)}}$$

• $$\lim_{m\rightarrow\infty}\mu^{(m)}=\lambda_1$$

So the process of iteration is actually computing

$$\textbf{y}^{(m)}\rightarrow p_m\rightarrow \mu^{(m)}\rightarrow \textbf{y}^{(m+1)}\rightarrow...$$

if $\mu^{(m)}=0$, then output $\textbf{x}^{(m-1)}$ and 0 immediately.

# Symmetric Power Method

If $\textbf{A}$ is symmetric:

• Choose an arbitrary unit vector $||\textbf{x}||_2^{(0)}$.
• Let $\textbf{y}^{(m)}=\textbf{Ax}^{(m-1)}$.
• Set $\mu^{(m)}=(\textbf{x}^{(m-1)})^T\textbf{y}^{(m)}$.
• let $\textbf{x}^{(m)}=\textbf{y}/||\textbf{y}||_2$.
• Repeat computing $\mu$.

# Inverse Power method

Noticing that the two methods above can only compute the biggest eigenvalue.

And we can compute the eigenvalue which is most close to $q$ by applying the method above to $(\textbf{A}-q\textbf{I})^{-1}$ since its eigenvalues can be sorted as:

$$\frac{1}{\lambda_1-q},\frac{1}{\lambda_2-q},...,\frac{1}{\lambda_n-q}$$

# Rayleigh Quotient Iteration

$$\because \textbf{x}^T\textbf{Ax}=\lambda\textbf{x}^T\textbf{x}\\ \therefore \lambda=\frac{\textbf{x}^T\textbf{Ax}}{||\textbf{x}||_2^2}$$

It’s useful to accelerate the convergence of a method.

By computing:

$$q^{(m)}=\frac{(\textbf{x}^{(m-1)})^T\textbf{Ax}^{(m-1)}}{||\textbf{x}^{(m-1)}||_2^2}$$

and solve the linear system $(\textbf{A}-q^{(m)}\textbf{I})\textbf{y}^{(m)}=\textbf{x}^{(m-1)}$ to get $\textbf{y}^{(m)}$.

and we can converge to $\frac{1}{\lambda-q}$ more quickly.

# QR Factorization

# Householder transformation

Definition: Householder transformation has form

$$\textbf{H}=\textbf{I}-2\frac{\textbf{vv}^T}{\textbf{v}^T\textbf{v}}$$

for nonzero vector $\textbf{v}$.

Property: $\textbf{H}=\textbf{H}^T=\textbf{H}^{-1}$.

Given a vector $\textbf{a}$, we want to find a vector $\textbf{v}$ so that

$$\textbf{Ha}=\begin{pmatrix}\alpha\\0\\0\\...\\0\end{pmatrix}=\alpha\textbf{e}_1$$

the vector is

$$\textbf{v}=(\textbf{a}-\alpha\textbf{e}_1)\frac{\textbf{v}^T\textbf{v}}{2\textbf{v}^T\textbf{a}}$$

and $\alpha=-sign(a_1)||\textbf{a}||_2$, where $a_1$ is the first element of $\textbf{a}$.

---

The process of QR factorization:

Consider a matrix \textbf{A}=\begin{pmatrix}2&2\\1&4\\2&6\end

• The first column of it is $\textbf{u}_1$, then choose \textbf{v}=\textbf{u}_1-(-||\textbf{v}_1||_2\textbf{e}_1)=\begin{pmatrix}5\\1\\2\end

• Then $\textbf{H}=\textbf{I}-2\frac{\textbf{vv}^T}{\textbf{v}^T\textbf{v}}$

  $$\textbf{HA}=\begin{pmatrix}-3&-6.66\\0&2.26\\0&2.53\end{pmatrix}$$

• choose the second column of $\textbf{HA}$ and repeat the process.

# Givens Rotations

Consider a matrix $\textbf{A}$, if we want it to be up triangle matrix.

so we need to change all $a_{ij}(i>j)$ to zero.

consider the rotate operator:

$$\textbf{R}=\begin{pmatrix}1&0&0&0&0\\ 0&c&0&s&0\\ 0&0&1&0&0\\ 0&-s&0&c&0\\ 0&0&0&0&1 \end{pmatrix}$$

It will only work on the second and the fourth row of the matrix.

$$\textbf{RA}=\begin{pmatrix}a_{11}&a_{12}&a_{13}&a_{14}&a_{15}\\ ca_{21}+sa_{41}&ca_{22}+sa_{42}&ca_{23}+sa_{43}&ca_{24}+sa_{44}&ca_{25}+sa_{45}\\ a_{31}&a_{32}&a_{33}&a_{34}&a_{35}\\ ca_{41}-sa_{21}&ca_{42}-sa_{22}&ca_{43}-sa_{23}&ca_{44}-sa_{24}&ca_{45}-sa_{25}\\ a_{51}&a_{52}&a_{53}&a_{54}&a_{55}\\ \end{pmatrix}$$

So the basic idea is that, by applying rotate operators, we can change $\textbf{A}$'s first column into $\alpha\textbf{e}_1$, then change the second row,…

since the first column is all zero when $i>1$, so changing the second row won’t influence the first column.