第2章 方程求根实验

2.1 实验目的

1、理解迭代法原理，并利用迭代法求非线性方程的根
2、熟悉MATLAB软件种非线性方程（组）的求解命令及用法
3、使用MATLAB解决一些非线性方程问题

2.2 迭代方法

一、理论知识

1、迭代格式
2、零点定理
3、压缩映射原理
4、迭代收敛判别

---

---

二、迭代算法设计

1、二分法

二分法原理：将区间二分，根据零点定理，判断根在某个分段内，再进行二分，依次类推，重复进行，直到满足精度为之。

具体算法不赘述

例2-1 使用二分法求方程$x^3+x-1=0在[0,1]$内的近似根（误差$<10^{-5}$）。

function y = exam2_1(m,n,er)
%二分迭代程序[m，n]表示迭代初始区域，er表示误差
%y为向量，其分量表示在迭代过程中的各个区间中点
syms x xk
a = m;b = n;k = 0;
ff = x^3 + x - 1;
while b - a > er
       xk = (a + b)/2;
       fx = subs(ff,x,xk);
       fa = subs(ff,x,a);
       k = k +1;
       if fx == 0
           y(k) = xk;
           break;
       elseif fa*fx<0
           b = xk;
       else a = xk;
       end
       y(k) = xk;
end
plot(y,'.-');
grid on

在命令行窗口下执行：

>> ab = exam2_1(0,1,1e-5);
>> vpa(ab,8)

可以在迭代区间中点数列ab

ans =
[ 0.5, 0.75, 0.625, 0.6875, 0.65625, 0.671875, 0.6796875, 0.68359375, 0.68164062, 0.68261719, 0.68212891, 0.68237305, 0.68225098, 0.68231201, 0.68234253, 0.68232727, 0.6823349]

图2-1 二分迭代法（区间中点数列）

---

2、简单迭代法

将方程$f(x)=0$，经过简单变形处理，按照式子$x_{n+1}=\phi (x_n)$设置迭代格式，若迭代收敛，即收敛到不动点$x^{\ast }$。

例2-2 设计不同迭代算法求方程$x^2=2$的正根

$方法1：x^2=2<=>x=\frac{2}{x}$

$方法2：x^2=2<=>x=\frac{2}{x}<=>x+x=x+\frac{2}{x}<=>\frac{1}{2}(x+\frac{2}{x})$

$方法3：x^2=2<=>x^3=2x=>x^3=2x-x^2+2$

$方法4：x^2=2<=>x^2+x=2+x<=>x=\frac{2+x}{1+x}<=>x=1+\frac{1}{1+x}$

选用初值x=1，同时对以上三种格式进行迭代计算

function [y1,y2,y3,y4] = exam2_2(n)
%四种迭代方式求x^2 - 2 = 0 的根；n表示迭代次数
y1 = ones(1,n);y2 = ones(1,n);y3 = ones(1,n);y4 = ones(1,n);
for i = 1:n
    y1(i + 1) = 2/y1(i);
    y2(i + 1) = 0.5*(y2(i) + 2/y2(i));
    y3(i + 1) = (2 * y3(i) - y3(i)^2 + 2)^(1/3);
    y4(i + 1) = 1 + 1/(1 + y4(i));
end

序号	1	2	3	4	5	6	7	8	9	10	11
方法1	1	2	1	2	1	2	1	2	1	2	1
方法2	1	1.5	1.416666667	1.414215686	1.414213562	1.414213562	1.414213562	1.414213562	1.414213562	1.414213562	1.414213562
方法3	1	1.44224957	1.410200216	1.414764791	1.414137399	1.414224077	1.414212111	1.414213763	1.414213535	1.414213566	1.414213562
方法4	1	1.5	1.4	1.416666667	1.413793103	1.414285714	1.414201183	1.414215686	1.414213198	1.414213625	1.414213552

例2-3 设置两种不同迭代格式求解以下方程：

代码无自由发挥

$$(1)x^3-x^2-x-1=0;(2)5x-e^x=0。$$

---

3、牛顿法

详细定义略，可以百度自行了解，或观看例题解法

例2-4 使用牛顿迭代求方程$x^3+x-1=0$再$[0,1]$内的近似根。

function x = exam2_4(xx,n)
%牛顿迭代法，xx为迭代初值，n为迭代步数
x = zeros(1,n+1);
x(1) = xx;
for i = 1:n
    x(i+1) = x(i) - (x(i)^3 + x(i) - 1)/(3 * x(i)^2 + 1);
end

命令窗口下执行：

>> format long
>> x=exam2_4(1,10)

x =

  列 1 至 5

   1.000000000000000   0.750000000000000   0.686046511627907   0.682339582597314   0.682327803946513

  列 6 至 10

   0.682327803828019   0.682327803828019   0.682327803828019   0.682327803828019   0.682327803828019

  列 11

   0.682327803828019

可以看出当前牛顿迭代收敛速度快于二分法。

缺点：

1、导数计算比较困难

2、在根的领域内导数值非常小

3、方程有多重根，牛顿迭代法可能失效

---

4、弦截法

详细定义略，可以百度自行了解，或观看例题解法

例2-5 使用两点弦截法求方程$x^3-3x-1=0在x_0=2$内的近似根（误差$<10^{-5},取x_0=2,x_1=1.9$）。

function [y,k] = exam2_5(er,n,xa,xb)
%两点弦截法，er表示误差，n表示步数，xa、xb表示迭代初值
%y为向量，其分量表示在迭代过程中的各个点，k表示实际迭代步数
syms x xk
x0 = xa;x1 = xb;
ff = x^3 - 3*x -1;
y(1) = xa;y(2) = xb;
k = 2;
while abs(x1 - x0) > er & k<n
    fx1 = subs(ff,x,x1);
    fx0 = subs(ff,x,x0);
    x2 = x1 - fx1 * (x1 -x0)/(fx1 - fx0);
    k = k+1;
    y(k) = x2;
    x0 = x1;
    x1 = x2;
end

命令窗口执行：

>> format long
>> [y,k] = exam2_5(1e-5,100,2,1.9)
y =
  列 1 至 5

   2.000000000000000   1.900000000000000   1.881093935790725   1.879411060169918   1.879385274283925

  列 6

   1.879385241572444
   
k =
     6

上述迭代过程标明：使用两点弦截法，迭代k-2=4步，即可达到误差限内近似解x=1.879385241572444

---

---

三、迭代算法加速

1、松弛法

直接上例题

例2-6 已知：$(1)x^3-x^2-x-1=0$$(2)5x-e^x=0$

function [y1,y2,z1,z2] = exam2_6(n)
%松弛加速迭代方式；n表示迭代次数；
%y1，y2表示有关方程x^3-x^2-x-1=0的迭代；
%z1，z2表示有关方程5x-exp(x)=0的迭代
y1 = ones(1,n+1);y2 = ones(1,n+1);z1 = ones(1,n+1);z2 = ones(1,n+1);
y1(1) = 3;y2(1) = 3;
for i = 1 : n
    y1(i+1) = y1(i)^3 - y1(i)^2 -1;
    y2(i+1) = (-2*y2(i)^3 + y2(i)^2 -1)/(-3*y2(i)^2 + 2*y2(i) + 1);
    z1(i+1) = exp(z1(i))/5;
    z2(i+1) = (1 - z2(i))*exp(z2(i))/(5 - exp(z2(i)));
end

命令窗口执行：

>> [y1,y2,z1,z2] = exam2_6(7)

序号	1	2	3	4	5	6	7	8	…
y1	3	17	4623	98781980237	9.64E+32	8.96E+98	7.18E+296	NaN	…
y2	3	2.3	1.951703992	1.848473773	1.839355695	1.839286759	1.839286755	1.839286755	…
z1	1	0.543656366	0.344458539	0.282245118	0.265220746	0.260743747	0.259579007	0.259276841	…
z2	1	0	0.25	0.259156526	0.259171102	0.259171102	0.259171102	0.259171102	…

---

2、埃特金$(Aitken)$法

直接上例题

例2-7 已知$x=2$为方程$x^3-5x+2=0$的一个根，设计一个埃特金加速迭代格式，并进行数值验证。

function [y1,y2] = exam2_7(xx,n)
%Aitken加速迭代方式；xx表示迭代初值；n表示迭代次数；
%y1 = (x^3 + 2)/5；y2相对应的Aitken加速迭代
y1 = ones(1,n+1);
y2 = ones(1,n+1);
z1 = ones(1,n+1);
z2 = ones(1,n+1);
y1(1) = xx;
y2(1) = xx;
syms x
ff = (x^3 + 2)/5;
for i = 1:n
    y1(i+1) = subs(ff,x,y1(i));
end
for i = 1:n
    a = subs(ff,x,y2(i));
    b = subs(ff,x,a);
    y2(i+1) = b - (b-a)^2/(b-2*a + y2(i));
end

命令窗口执行：

>> [y1,y2] = exam2_7(2.5,5)

序号	1	2	3	4	5	6
y1	2.5	3.525	9.160065625	154.118363	7.32E+05	7.85E+16
y2	2.5	2.272101943	2.102171217	2.018093791	2.000651733	2.000000873