弦截法

例题

求 $xe^x - 1 = 0$ 的根，取初值$x_0 = 0.5, x_1 = 0.6。$

实现

其数学迭代公式为：
$x_{k + 1} = x_k - \frac{f(x_k)} {f(x_k) - f(x_{k - 1})} \cdot (x_k - x_{k - 1})$

使用Matlab实现：

f = inline('x * e^x - 1');
x = 0;
xx = 0.5;
xxx = 0.6;
i=0;
while abs(xx-xxx) >= 1e-5x = xx;xx = xxx;xxx = xx - f(xx) / (f(xx) - f(x)) * (xx - x);i = i + 1;
end

经过迭代4次，得到$x_5 = 0.56714$。
如果需要看到更多小数位数，可以使用：

format long g

二分法

例题

用二分法求方程 $2e^{-x} - sinx = 0$ 在区间 $[0, 1]$ 内的根，要求 $| x_k - x^* | < \frac12 \times 10^{-5}$ 。

实现

首先通过求端点值，求导，确定左边函数在 $[0, 1]$ 有解，且是单根区间。
取 $c = 0, d = 1$ ，将 $[c, d]$ 对分，设对分点为 $x_0 = \frac12 (c + d)$ ，计算 $f(x_0)$ ，如果 $f(x_0)$ 与 $f(c)$ 同号，……
$d_n - c_n = \frac1{2^n} (d - c)$ ，只要n足够大，有根区间长度就会足够小，当长度达到精度要求时，取 $x_n = \frac12 (c_n + d_n)$ 作为方程的根的近似值。

使用Matlab实现：

f = inline('2 * e^(-x) - sin(x)');
c = 0;
d = 1;
x = 1/2 * (c + d);
i = 0;
while d - c >= 1e-5/2x = 1/2 * (c + d);i ++;if f(x) > 0c = x;elseif f(x) < 0d = x;elsebreak;end;
end;

经过18次迭代，$c = x_{17}$ 时，$d - c < \frac12 \times 10^{-5}$ ，得到原方程的解为：$x^* = \frac12 (c_{18} + d_{18}) = 0.92103$ 。

CG法

例题

取 $x_0 = (0, 0)^T，用CG法求解：$
$\left[\begin{array}{cc} 6 & 3\\ 3 & 2\\ \end{array}\right] \cdot \left[\begin{array}{c} x_1\\ x_2\\ \end{array}\right] = \left[\begin{array}{c} 0\\ -1\\ \end{array}\right]$

实现

方程组化为：$A \cdot x = b$ 。

取 $x_0 = (0, 0)^T，r_0 = b - Ax_0 = (0, -1)^T，p_0 = r_0$

$\alpha_0 = \frac{(r_0, r_0)}{(p_0, Ap_0)} = \frac12$

$x_1 = x_0 + \alpha_0 \cdot p_0 = (0, -\frac12)^T$

$r_1 = r_0 - \alpha_0 \cdot A \cdot p_0 = (\frac32, 0)^T$

$\beta_0 = \frac{(r_1, r_1)}{(r_0, r_0)} = \frac94$

$p_1 = r_1 + \beta \cdot p_0 = (\frac32, -\frac94)^T$

$\alpha_1 = \frac{(r_1, r_1)}{(p_1, Ap_1)} = \frac23$

$x_2 = x_1 + \alpha_1 \cdot p_1 = (1, -2)^T$

$r_2 = r_1 - \alpha_1 \cdot A \cdot p_1 = (0, 0)^T$

在$r_k = 0 或 (p_k, Ap_k) = 0$ 时，计算终止，有 $x_k = x^*$ 。

使用Matlab实现：

x = [0;0];
A = [6,3;3,2];
b = [0;-1];
r = b - A * x;
p = r;
i = 0;
while !all(r == [0;0]) && !all(dot(p, A * p) == [0;0])alpha = dot(r, r) / dot(p, A * p);xx = x + alpha * p;rr = r - alpha * A * p;beta = dot(rr, rr) / dot(r, r);pp = rr + beta * p;i ++;r = rr;x = xx;p = pp;
end;