Bisection/Newton/Fixed Point

Bisection 二分法求零点

我们使用二分法求零点，根据函数的连续性，假如在一个函数在interval[a,b]上连续，而且在这个interval的两端取得不同的符号，那么显然我们在这个区间上至少有一个零点，如果函数是单调的，那么在这个interval上的零点唯一。

我们可以用二分法进行求解，需要给定interval的两端[a,b]的值，以及tolerance的值，通过区间的长度和tolerance来计算需要迭代的次数n，然后进行n次迭代后输出。

Newton 牛顿法和modified牛顿法

二分法相当于只使用了函数interval两端的信息，但是我们可以有更快的选择。牛顿法是一种收敛速度比二分法更快的方法，可以描述为在x^(k)上做一条与原函数相切的直线，将这条直线与x轴的焦点做为下次迭代的x^(k+1)。

和二分法不同，尽管牛顿法更快，但是牛顿法并不是全局收敛的，仅仅对于零点a附近某个区间内的点是收敛的，并且我们不能判断给定一个点x₀以及对应的函数值导数值等，去判断能否从x₀收敛到a。

我们既然提到了牛顿法是收敛更快的算法，那么要快多少呢？二分法的收敛是线性的，因为我们可以视error为e^(k)每次迭代都变为e^(k+1)都是原来的一半。而牛顿法，如果零点a是simple的（即f’(a)不等于0），那么牛顿法是quadratically收敛的。

如果零点a不是simple的，且m是a的multiplicity，如果想让牛顿法收敛，我们需要确保x⁽⁰⁾的选择是合理的，且在interval上没有其他x使得f’(x)=0，此时牛顿法的收敛为线性的，如果我们使用modified牛顿法（m对应multiplicity的值），此时modified牛顿法的收敛依然是quadratically的。

如果f在interval[a,b]上是单调的，且f(a)和f(b)异号，那么对任何[a,b]上的x₀，假如x₀满足f(x₀)*f’’(x₀)>0，那么牛顿法在x₀上收敛。

Fixed Point定点法

定点法：x^(k+1)=phi(x^(k))

根据上图不难得出结论，假如a是phi的一个fixed point(即有a=phi(a))，lim e^(k+1)/e^(k) = phi’(a),我们希望error在这个迭代的过程中是不断递减的，那么必须有phi’(a)<1。那么我们根据a的关于phi的导数，有如下三种情况

1. |phi’(a)|<1 :存在一个r，对于任意的(a-0,a+0)中的x0定点法收敛（注意，此处的定点法收敛指的是对于合适的x₀定点法收敛）
2. |phi’(a)|>1 :存在至少1个x₀不收敛到a（即定点法不收敛）
3. |phi’(a)|=1 :可能收敛也可能不收敛

对于phi，我们可以根据phi在fixed point a上的导数来判断收敛的速度，如果phi⁽ⁱ⁾=0, i=1,2…p-1，且phi^(p)!=0，那么定点法x^(k+1)=phi(x^(k))的order为p，即lim e^(k+1)/(e^(k))^p=constant。

Convergence Order

判断convergence order，本质上就是判断constant=e^(k+1)/(e^(k))^p中p的值，如果p为1，则为线性的，如果p为2，则为quadratically的，以此类推。通过求log可以很快的得出结论并画出图像。

LU分解

Ax=b, A是一个n*n矩阵，b是一个n维向量，求解x。

Backward and forward substitution

A=LU，L是一个下三角矩阵，而U是一个上三角矩阵，将这个式子和上面的式子连立起来，能得到Ly=b,Ux=y。

这样做的好处是在正常情况下，为了求解x的值，需要令b右乘A的逆矩阵，即x=A^-1b，但是求逆矩阵是非常复杂的，需要很高的computing power，为了避免求逆矩阵，我们就采用LU分解。

L(U)是下(上)三角矩阵，对于L使用forward substitution，因为显然从方程的角度去看这个问题，我们有L_1,1 y₁=b₁，从而可以一一解出y的值。 对于U使用backward substitution，解答的方法同上，但是由于我们是从U最下面的元素开始解出x的值的U_n,n x_n=y_n，所以将这个方法称为backward substitution。

LU分解存在的条件

不是每一个A都可以分解成LU的形式，那么什么样的矩阵可以分解为LU的形式呢？显然对于一个不可逆矩阵，是无法分解成LU的形式的，但是对于一个可逆矩阵，也可能出现无法分解成LU的情况。

1. 如果A是一个row/column对角线dominant的矩阵（即a对角线上的元素大于同一row/column中所有其他元素的绝对值之和），矩阵A能够LU分解。
2. A是一个symmetric positve definite矩阵，则A能够进行LU分解。

判断A是否是正定矩阵，我们可以通过eig(A)>0来判断，eig(A)表示A的特征值，而正定矩阵的特征值是全正的。

以上的两种情形是A能够进行LU分解的特殊情况，一个更general的条件是，给定矩阵A，其LU分解存在且唯一的充要条件是A对角线上的所有方阵A_i都是非奇异矩阵(行列式不等于0)

General LU Decomposition

LU分解是将一个矩阵分解为L*U，其中L的对角线上的元素为1。具体的操作过程可以视为将其进行不断的transform，不断的进行每个row之间的数乘和加减。

对于特殊的A进行特殊的LU分解-Cholesky分解

如果A是一个symmetric positive definite的矩阵，那么A可以分解为H*H^T，根据其对称的形式，可以很简单的用矩阵表示出H中的每个元素。

对于特殊的A进行特殊的LU分解-Thomas Algorithm

如果A是一个对称矩阵，并且只有三列对角线上的元素，那么我们可以用类似上面的general LU decomposition的方法进行分解，但是每一个value的表示形式会变得更加简便。

建立一个对角线矩阵

a = [1:10]’;

c = [10:19]’;

e = [102:111]’;

n = length(a);

A = spdiags([e a c], [−1 0 1], n, n)

LU with pivoting

Pivoting是一种改变改变矩阵的行的顺序的方法，使得一个不能进行LU分解的矩阵可以进行LU分解。表示为Ax=b,PA=LU。

infinity norm of the error: norm(x-x1,inf)

infinity norm of the normalized residual error: norm(b-A*x1, inf)/norm(b, inf)

Condition Number

P-norm: ||A||_p=(SIGMA|x_i|^p)^1/p

Condition number: K_{p(A)=||A||p*||A^-1||p}

condition number: cond(A,n)

K_{p(A)是衡量sensitivity的一个指标，即对于b的微小的扰动，x的扰动程度。如果Kp(A)越小，表示A越well-conditioned。}

Iterative Methods of Linear system

Ax=b, A是一个n*n矩阵，b是一个n维向量，求解x。

Stationary methods

x^(k+1) = B*x^(k)+g

s.t. x^* = Bx^\+g

我们设置stopping criterion为

1. ||x^(k+1)-x^(k)||<=tolerance
2. ||b-Ax^(k+1)||<=tolerance

这种方法收敛的充要条件是ρ(B)<1，ρ(B)=max(|eig(B)|)。

同时，为了估计当前迭代状态下的误差error ||x^k-x^*||，我们可以通过||x^(k+1)-x^k||来约束这个误差。

目标：构建一个A=P-(P-A)。

Jacobi method & Gauss Seidel method

Richardson method

不难看出，其实Jacobi method和Gauss Seidel method都是Richardson method的特殊化的形式，就像是牛顿法本质上是定点法的特殊演绎一样。

其中α的最佳选择α_opt=2/(λ_min+λ_min)，λ_min/λ_min表示P^-1A的最小/大的特征值。

Gradient method* 此节不在考试范围内

Interpolation

interpolation的直译结果为【插值】，但是其实在这门课所涉及到的地方，将其理解为【拟合】这个概念会更好理解一些，从原函数f(x)上取得一系列的点，然后根据这些点来设计一个新的函数f*(x)来近似这个原函数f(x)。

拉格朗日多项式插值

对于任意n个nodes{x_{i,yi}，存在唯一的一个多项式f(x)，其degree小于等于n，且对于任何xi，都有f(xi)=yi。}

根据上面的公式，为了获取f(x)的表达形式，我们设计φ_k(x_j)=1 if j=k, 0 else。然后求和φ_k(x_j) f(x_k)，得到f*(x)。

分段线性插值

这种插值的方式就很好理解了，我们将函数分为很多段，然后将每一段用一条直线连接起来，相比于上面，虽然这样得到的f*(x)不是光滑的，但是也是原函数的一种近似的方式。

Spline插值

与分段式线性插值有一些类似，我们将函数分为很多段，但是我们还希望让我们的新的f*(x)是光滑的，在上面的近似方式中，不光滑的点仅仅存在在两个interval的连接处，那么我们在spline的方法中，对每个interval [x_i,x_j]都使用一个三次函数（注意！不再是上面的线性/一次函数了）进行近似，我们可以调整这个函数，使其不仅经过x_i和x_j对应的点，同时在这两个点上的导数与前一个以及后一个interval是相同的。

Integral 积分

midpoint rule和trapezoidal rule都是order 1的，这意味着它们可以近似线性函数而没有任何误差，而Simpson rule的degree of exactness等于3，这意味着order低于3的function的积分值是准确的。