1 第七讲 MATLAB中求方程的近似根(解)
教学目的：学习matlab中求根命令，了解代数方程求根求解的四种方法，即图解法、准解析法、数值方法以及迭代方法，掌握对分法、迭代法、牛顿切法线求方程近似根的基本过程；掌握求代数方程（组）的解的求解命令．
教学重点：求方程近似解的几种迭代方法，代数方程（组）的解的求解命令的使用方法.利用所学的编程知识，结合具体的实例，编制程序进行近似求根．掌握相关的代数方程（组）的求解命令及使用技巧．
教学难点：方程的近似求解和非线性方程（组）的求解．
一、问题背景和实验目的
求代数方程0)(xf的根是最常见的数学问题之一（这里称为代数方程，主要是想和后面的微分方程区别开．为简明起见，在本实验的以下叙述中，把代数方程简称为方程），当)(xf是一次多项式时，称0)(xf为线性方程，否则称之为非线性方程．
当0)(xf是非线性方程时，由于)(xf的多样性，尚无一般的解析解法可使用，但如果对任意的精度要求，能求出方程的近似根，则可以认为求根的计算问题已经解决，至少能满足实际要求．同时对于多未知量非线性方程(组)而言，简单的迭代法也是可以做出来的，但在这里我们介绍相关的命令来求解，不用迭代方法求解．
通过本实验，达到下面目的：
1. 了解对分法、迭代法、牛顿切线法求方程近似根的基本过程；
2. 求代数方程（组）的解．
首先，我们先介绍几种近似求根有关的方法：
1． 对分法
对分法思想：将区域不断对分，判断根在某个分段内，再对该段对分，依此类推，直到满足精度为止．对分法适用于求有根区间内的单实根或奇重实根．
设)(xf在],[ba上连续，0)()(bfaf，即 ()0fa，()0fb或()0fa，()0fb．则根据连续函数的介值定理，在),(ba内至少存在一点 ，使()0f．
下面的方法可以求出该根： 2 (1) 令0()/2xab，计算0()fx；
(2) 若0()0fx，则0x是()0fx的根，停止计算，输出结果0xx．
若 0()()0fafx，则令1aa，10bx，若0()()0fafx，则令10ax，1bb；111()/2xab．……，有ka、kb以及相应的()/2kkkxab．
(3) 若()kfx (为预先给定的精度要求)，退出计算，输出结果()/2kkkxab；
反之，返回(1)，重复(1)，(2)，(3)．
以上方法可得到每次缩小一半的区间序列{[,]}kkab，在(,)kkab中含有方程的根．
当区间长kkba很小时，取其中点()/2kkkxab为根的近似值，显然有
2111()/2()/(2)()/2kkkkkkxbababa
以上公式可用于估计对分次数k．
分析以上过程不难知道，对分法的收敛速度与公比为12的等比级数相同．由于1021024，可知大约对分10次，近似根的精度可提高三位小数．对分法的收敛速度较慢，它常用来试探实根的分布区间，或求根的近似值．
2. 迭代法
a) 松弛法：由方程()0fx构造一个等价方程
()xx．
则迭代方程是：
1(1)()kkkkkxxx，1/(1'())kkx，其中'()1x．
松弛法的加速效果是明显的 (见附录4)，甚至不收敛的迭代函数经加速后也能获得收敛．
b) Altken方法：
松弛法要先计算'()kx，在使用中有时不方便，为此发展出以下的 Altken 公式：
(1)()kkxx ；
(2)(1)()kkxx；
(2)(2)(1)2(2)(1)1()/(2)kkkkkkkxxxxxxx，,2,1,0k
这就是Altken 公式，它的加速效果也是十分明显的，它同样可使不收敛的迭代格式获得收敛(见附录5)．
3. 牛顿(Newton)法（牛顿切线法）
()0fx是非线性方程
其迭代公式为： 3 1(()/'())kkkkxxfxfx ,2,1,0k
即为牛顿法公式.
牛顿法的缺点是：(1)对重根收敛很慢；(2)对初值0x要求较严，要求0x相当接近真值*x．
因此，常用其他方法确定初值0x，再用牛顿法提高精度．
以下是本实验中的几个具体的实验
具体实验1：对分法
先作图观察方程：3310xx的实根的分布区间，再利用对分法在这些区间上分别求出根的近似值．程序如下：
function [y,p]=erfen()
clc, x=[];a=[];b=[]; a(1)=1;b(1)=2; i=1;x(i)=(a(i)+b(i))/2; e=abs(f(x(i)));
ezplot('x^3-3*x+1',[a(1),b(1)]);hold on, plot([a(i),b(i)],[0,0])
while e>10^(-5)
plot([a(i),a(i)],[0,100],[x(i) x(i)],[0 100],[b(i) b(i)],[0 100]),pause(0.5)
if f(a(i))*f(x(i))<0
a(i+1)=a(i);b(i+1)=x(i);x(i+1)=(a(i+1)+b(i+1))/2;
else
a(i+1)=x(i);b(i+1)=b(i);x(i+1)=(a(i+1)+b(i+1))/2;
end
e=abs(f(x(i)));i=i+1;
end
y=x(i);p=[a;x;b]'
function u=f(x)
u=x^3-3*x+1;
end
end
图形如下:
结果为：1.5321 4 具体实验2：普通迭代法
采用迭代过程：1()kkxx求方程3310xx在 0.5 附近的根，精确到第 4 位小数．
构造等价方程：3(1)/3xx
用迭代公式： 31(1)/3kkxx，,2,1,0k
具体实验3：迭代法的加速1——松弛迭代法
3()(1)/3xx，2()'xx，21/(1)kkx
迭代公式为
31(1)(1)/3kkkkkxxx
clc;x=[];w=[]; x(1)=1;w(1)=1/(1-x(1));
for i=1:10
w(i)=1/(1- x(i)); x(i+1)=(1-w(i))*x(i)+ w(i)*(x(i)^3+1)/3;
end
x
另外有程序可以参考，详见参见附录4．
具体实验4：迭代法的加速2——Altken迭代法
迭代公式为：
(1)3(1)/3kkxx，(2)(1)3(1)/3kkxx
(2)(2)(1)2(2)(1)1()/(2)kkkkkkkxxxxxxx，,2,1,0k
%（符号计算）
syms x fx gx;
gx=(x^3+1)/3;fx=x^3-3*x+1; disp('k x x1 x2')
x=0.5;k=0; ffx=subs(fx, 'x', x);
while abs(ffx)>0.0001;
u=subs(gx, 'x', x);v=subs(gx, 'x', u);
disp([num2str(k), ' ', num2str(x), ' ', num2str(u), ' ', num2str(v)])
x=v-(v-u)^2/(v-2*u+x);k=k+1;ffx=subs(fx, 'x', x);
end
disp([num2str(k), ' ', num2str(x), ' ', num2str(u), ' ', num2str(v)])
%（数值计算）
function [y,p]=althken()
% 求方根的迭代程序
clc,format long e, x(1)=6; i=1;p=[];ezplot('x^3-3*x+1',[x(1)-9,x(1)+1]);hold on
plot([x(1)-20,x(1)+2],[0,0])
while abs(f(x(i)))>=10^(-5)
plot(x(i),0,'*') 5 t1=phi(x(i));t2=phi(t1); x(i+1)=t2-(t2-t1)^2/(t2-2*t1+x(i)+eps);
p=[p;[i, x(i),t1,t2]]; i=i+1; pause(0.1)
end
p,y=x(i), i, format
function u=phi(x)
u=(x^3+1)/3;
end
function u=f(x)
u=x^3+1-3*x;
end
end
具体实验5：牛顿法
用牛顿法计算方程3310xx在-2到2之间的三个根．
提示：3()31fxxx，2'()33fxx迭代公式：
2321(31)/(33)kkkkkxxxxx
function [y,p]=newton()
% 求方根的迭代程序
clc,format long e, x(1)=6; i=1; p=[]; ezplot('x^3-3*x+1',[x(1)-9,x(1)+1]);hold on
plot([x(1)-20,x(1)+2],[0,0])
while abs(f(x(i)))>=10^(-5)
plot(x(i),0,'*'), x(i+1)=x(i)-f(x(i))/(df(x(i))+eps);
p=[p;[i, x(i)]]; i=i+1; pause(0.1)
end
format short, p,y=x(i), i,
function u=df(x)
u=3*x^2-3;
end
function u=f(x)
u=x^3+1-3*x;
end
end
结果：