使用Matlab实现：弦截法、二分法、CG法，求零点、解方程

弦截法

例题

求 $x{e}^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-5
	x = 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^{\ast }|<\frac{1}{2}\times 10^{-5}$ 。

实现

首先通过求端点值，求导，确定左边函数在 $[0,1]$ 有解，且是单根区间。
取 $c=0,d=1$ ，将 $[c,d]$ 对分，设对分点为 $x_0=\frac{1}{2}(c+d)$ ，计算 $f(x_0)$ ，如果 $f(x_0)$ 与 $f(c)$ 同号，……
$d_n-c_n=\frac{1}{2^n}(d-c)$ ，只要n足够大，有根区间长度就会足够小，当长度达到精度要求时，取 $x_n=\frac{1}{2}(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/2
    x = 1/2 * (c + d);
    i ++;
    if f(x) > 0
        c = x;
    elseif f(x) < 0
        d = x;
    else
        break;
    end;
end;

经过18次迭代，$c=x_{17}$ 时，$d-c<\frac{1}{2}\times 10^{-5}$ ，得到原方程的解为：$x^{\ast }=\frac{1}{2}(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-A{x}_0=(0,-1{)}^T，p_0=r_0$

${\alpha }_0=\frac{(r_0,r_0)}{(p_0,A{p}_0)}=\frac{1}{2}$

$x_1=x_0+{\alpha }_0\cdot p_0=(0,-\frac{1}{2}{)}^T$

$r_1=r_0-{\alpha }_0\cdot A\cdot p_0=(\frac{3}{2},0{)}^T$

${\beta }_0=\frac{(r_1,r_1)}{(r_0,r_0)}=\frac{9}{4}$

$p_1=r_1+\beta \cdot p_0=(\frac{3}{2},-\frac{9}{4}{)}^T$

${\alpha }_1=\frac{(r_1,r_1)}{(p_1,A{p}_1)}=\frac{2}{3}$

$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,A{p}_k)=0$ 时，计算终止，有 $x_k=x^{\ast }$ 。

使用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;