第七讲 MATLAB中求方程的近似根(解)

1、第七讲MATLAB中求方程的近似根(解)教学目的：学习matlab中求根命令，了解代数方程求根求解的四种方法，即图解法、准解析法、数值方法以及迭代方法，掌握对分法、迭代法、牛顿切法线求方程近似根的基本过程；掌握求代数方程（组）的解的求解命令教学重点：求方程近似解的几种迭代方法，代数方程（组）的解的求解命令的使用方法.利用所学的编程知识，结合具体的实例，编制程序进行近似求根掌握相关的代数方程（组）的求解命令及使用技巧教学难点：方程的近似求解和非线性方程（组）的求解一、问题背景和实验目的求代数方程的根是最常见的数学问题之一（这里称为代数方程，主要是想和后面的微分方程区别开为简明起见，在本实验的以下

2、叙述中，把代数方程简称为方程），当是一次多项式时，称为线性方程，否则称之为非线性方程当是非线性方程时，由于的多样性，尚无一般的解析解法可使用，但如果对任意的精度要求，能求出方程的近似根，则可以认为求根的计算问题已经解决，至少能满足实际要求同时对于多未知量非线性方程(组)而言，简单的迭代法也是可以做出来的，但在这里我们介绍相关的命令来求解，不用迭代方法求解通过本实验，达到下面目的：1. 了解对分法、迭代法、牛顿切线法求方程近似根的基本过程；2. 求代数方程（组）的解首先，我们先介绍几种近似求根有关的方法：1 对分法对分法思想：将区域不断对分，判断根在某个分段内，再对该段对分，依此类推，直到满足精

3、度为止对分法适用于求有根区间内的单实根或奇重实根设在上连续，即 ，或，则根据连续函数的介值定理，在内至少存在一点 ，使下面的方法可以求出该根：(1) 令，计算；(2) 若，则是的根，停止计算，输出结果若 ，则令，若，则令，；，有、以及相应的(3) 若 (为预先给定的精度要求)，退出计算，输出结果；反之，返回(1)，重复(1)，(2)，(3)以上方法可得到每次缩小一半的区间序列，在中含有方程的根当区间长很小时，取其中点为根的近似值，显然有以上公式可用于估计对分次数分析以上过程不难知道，对分法的收敛速度与公比为的等比级数相同由于，可知大约对分10次，近似根的精度可提高三位小数对分法的收敛速度较慢，

4、它常用来试探实根的分布区间，或求根的近似值2. 迭代法a) 松弛法：由方程构造一个等价方程则迭代方程是：，其中松弛法的加速效果是明显的 (见附录4)，甚至不收敛的迭代函数经加速后也能获得收敛b) Altken方法：松弛法要先计算，在使用中有时不方便，为此发展出以下的 Altken 公式： ；，这就是Altken 公式，它的加速效果也是十分明显的，它同样可使不收敛的迭代格式获得收敛(见附录5)3. 牛顿(Newton)法（牛顿切线法）是非线性方程其迭代公式为： 即为牛顿法公式.牛顿法的缺点是：(1)对重根收敛很慢；(2)对初值要求较严，要求相当接近真值因此，常用其他方法确定初值，再用牛顿法提高精

5、度以下是本实验中的几个具体的实验具体实验1：对分法先作图观察方程：的实根的分布区间，再利用对分法在这些区间上分别求出根的近似值程序如下：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('x3-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)

6、 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; endy=x(i);p=a;x;b'function u=f(x)u=x3-3*x+1;endend图形如下:结果为：1.5321具体实验2：普通迭代法采用迭代过程：求方程在 0.5 附近的根，精确到第 4 位小数构造等价方程：用迭代公式： ，具体实验3：迭代法的加速1松弛迭代法，迭代公式为

7、 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;endx另外有程序可以参考，详见参见附录4具体实验4：迭代法的加速2Altken迭代法迭代公式为：，%（符号计算）syms x fx gx;gx=(x3+1)/3;fx=x3-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, '

8、;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);enddisp(num2str(k), ' ', num2str(x), ' ', num2str(u), ' ', num2str(v)%（数值计算）function y,p=a

9、lthken() % 求方根的迭代程序 clc,format long e, x(1)=6; i=1;p=;ezplot('x3-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,'*') 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, form

10、atfunction u=phi(x) u=(x3+1)/3; end function u=f(x) u=x3+1-3*x; end end具体实验5：牛顿法用牛顿法计算方程在-2到2之间的三个根提示：，迭代公式：function y,p=newton() % 求方根的迭代程序clc,format long e, x(1)=6; i=1; p=; ezplot('x3-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,'*'

11、), 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*x2-3; end function u=f(x) u=x3+1-3*x; endend结果：结果为： 1.5321 进一步思考：用迭代法求3的平方根. 迭代公式为. 编写M函数文件My_sqrt.m, 求3正的平方根. 要求误差小于仅要求写出源程序试使用以上介绍的迭代法来相互比较参考程序：function y=my_sqrt(a) % 求方根的迭代程序 if

12、 nargin=1|isa(a,'double') , error('输入数字为一个正数!'),endif a<0, error('输入数字为正数!'), endif a>0 format long e, x(1)=0; x(2)=1; i=1; while abs(x(i+1)-x(i)>=10(-5)i=i+1;x(i+1)=1/2*(x(i)+a/(x(i)+eps);end y=x(i+1);i,format end现在我们简单介绍图解法如何来求解一元方程和二元方程的根：例：exp(-3*t)*sin(4*t+2)+4*

13、exp(-0.5*t)*cos(2*t)=0.5>>ezplot('exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5',0 5)>>hold on, line(0,5,0,0)验证：t=3.5203>>syms x; t=3.5203; vpa(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans =例：x2*exp(-x*y2/2)+exp(-x/2)*sin(x*y)=0y2 *cos(y+x2) +x2*exp(x+y)=0>>

14、ezplot('x2*exp(-x*y2/2)+exp(-x/2)*sin(x*y)')>> hold on ezplot('y2 *cos(y+x2) +x2*exp(x+y)')具体的结果请大家自己下来运行二、 关于直接利用函数（命令）求解方程及简介 (1) solve('f(x)')，f(x)为一个具体的表达式(2) roots(A)，A为某个多项式按x降幂排列的系数矩阵(3) fzero('f(x)', x0)，f(x)为一个具体的表达式，x0为一个具体的数值(4) linsolve(A,b)，A为一方程组的系

15、数矩阵，b为方程组右端的常数矩阵1单变量的多项式方程求根：命令格式：roots(A)例：x3-6*(x2)-72*x-27=0;>>p=1 -6 -72 -27>>r=roots(p)r=12.1229-5.7345-0.38842. 多项式型方程的准解析解法命令格式： x,=solve(eqn1,eqn2, )例：x2+y2-1=0 0.75*x3-y+0.9=0>>syms x y;>> x,y=solve('x2+y2-1=0', '75*x3/100-y+9/10=0')检验： >>eval(&

16、#39;x.2+y.2-1'), eval('75*x.3/100-y+9/10')具体结果就请大家下来自己运行3. 线性方程组的求解例：求线性方程组的解，已知m=1 2 3 4 5;2 3 4 5 6;3 4 5 6 7 8;4 5 6 7 8 ;5 6 7 8 0，b=1;2;3;4;5 for i=1:5for j=1:5m(i, j)=i+j-1;endendm(5, 5)=0;b=1:5' linsolve(m, b)4. 非线性方程数值求解(1)单变量非线性方程求解在MATLAB中提供了一个fzero函数，可以用来求单变量非线性方程的根该函数的调用格

17、式为： z=fzero('fname',x0,tol,trace)其中fname是待求根的函数文件名，x0为搜索的起点一个函数可能有多个根，但fzero函数只给出离x0最近的那个根tol控制结果的相对精度，缺省时取tol=eps，trace指定迭代信息是否在运算中显示，为1时显示，为0时不显示，缺省时取trace=0例: 求f(x)=x-10x+2=0在x0=0.5附近的根步骤如下：(a) 建立函数文件funx.m function fx=funx(x)fx=x-10.x+2;(b)调用fzero函数求根 z=fzero('funx',0.5)z = 0.375

18、8(2)非线性方程组的求解对于非线性方程组F(X)=0，用fsolve函数求其数值解fsolve函数的调用格式为： X=fsolve('fun',X0,option)其中X为返回的解，fun是用于定义需求解的非线性方程组的函数文件名，X0是求根过程的初值，option为最优化工具箱的选项设定最优化工具箱提供了20多个选项，用户可以使用optimset命令将它们显示出来如果想改变其中某个选项，则可以调用optimset()函数来完成例如，Display选项决定函数调用时中间结果的显示方式，其中off为不显示，iter表示每步都显示，final只显示最终结果optimset(Dis

19、play,off)将设定Display选项为off例: 求下列非线性方程组在(0.5,0.5) 附近的数值解(a) 建立函数文件myfun.mfunction q=myfun(p)x=p(1);y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y);(b) 在给定的初值x0=0.5,y0=0.5下，调用fsolve函数求方程的根x=fsolve('myfun',0.5,0.5',optimset('Display','off')x = 0.6354 0.3734将求

20、得的解代回原方程，可以检验结果是否正确，命令如下：q=myfun(x)q = 1.0e-009 * 0.2375 0.2957可见得到了较高精度的结果精品案例：螺旋线与平面的交点问题： 螺旋线与平面相交的情况多种多样, 根据螺旋线与平面方程的不同可以相交, 也可以不相交. 在相交的情况下, 可以交于一点, 也可以交于好多点. 对于各种相交的情况, 要求其交点的坐标并不是一件容易的事. 本次实验就以此为背景讨论下面的具体问题：已知螺旋线的参数方程为 . 平面的方程为：. 求该螺旋线与平面的交点. 要求：1）求出所有交点的坐标；2）在同一图形窗口画出螺旋线、平面和交点. 实验过程： 1.1 问题分

21、析 可以采用多种方法求螺旋线与平面的交点坐标, 包括fsolve等. 先对方程化简,减少变量个数,使用图解方法求方程的根.再分别画出螺旋线,平面,及其交点.1.2 算法描述与分析先对方程化简,减少变量个数,再利用fsolve, 选择适当的初值, 求其数值解;再分别会出图形；最后对图形作出必要的修饰. 1.3 源程序及注释将螺旋线的参数方程代入平面方程后可得： 等价变形得 ：建立下面M文件intersect_point.m%使用图解法求交点，并且三维图%画图确定解的个数和大概位置theta=0:0.01:8*pi;y1=4*(cos(theta)+sin(theta);y2=2-0.5*thet

22、a;plot(theta,y1,theta,y2) %画出两个函数的图形%画螺旋线%theta=0:pi/100:8*pi;x=4*cos(theta);y=4*sin(theta);z=theta;figure %新建图形窗口plot3(x,y,z) %画含有参数的空间曲线hold on%透明的画平面%x1=-5:0.1:5; %取值和螺旋线的范围-4,4有关.y1=x1;X1 Y1=meshgrid(x1,y1);%网格化，画曲面Z1=4-2*X1-2*Y1;surf(X1,Y1,Z1) %或者使用mesh(X1,Y1,Z1)shading flatalpha(0.5) %设置透明度alp

23、ha('z') %设置透明度方向%求交点坐标,为避免变量混淆和覆盖，这里用t代替theta%i=1for n=2,5,9,11 %根据画图确定解的大概位置作为初值 t(i)=fsolve(inline('4*cos(t)+4*sin(t)+0.5 *t-2'),n)%选择不同初值求交点 x0(i)=4*cos(t(i); y0(i)=4*sin(t(i); z0(i)=t(i); i=i+1;endplot3(x0,y0,z0,'ro')1.4 测试结果（写清输入输出情况）从图形可见在 内与三角曲线有4个交点. 交点坐标为： theta的数值解为

24、：t=2.1961 5.3759 9.1078 11.1023四个交点的近似坐标为：x0 =-2.3413 2.4635 -3.8007 0.4261 y0 =3.2432 -3.1514 1.2468 -3.9772 z0 =2.1961 5.3759 9.1078 11.10231.5 调试和运行程序过程中产生的问题及采取的措施求交点的时候会出现重根和漏根的情形,通过选择适当的初值避免了上述情况.1.6 对算法和程序的讨论、分析, 改进设想及其它经验教训solve函数只能求解一个数值解,不能全部求出;用fsolve函数好; 为了满足更好的视觉效果,可以对图形进行进一步的修饰.习题1.已知多

25、项式2.解方程组:(1) (2)3.求解方程: 4.求解多项式方程 5.求下列代数方程（组）的解：(1) (2) 6.选择适当的迭代过程，分别使用：（1）普通迭代法；（2）与之相应的松弛迭代法和 Altken 迭代法求解方程 在 1.4 附近的根，精确到4位小数，请注意迭代次数的变化7.分别用对分法、普通迭代法、松弛迭代法、Altken 迭代法、牛顿切法线等5种方法，求方程 的正的近似根，(建议取 时间许可的话，可进一步考虑 的情况)五、附录为供近似求根的算法附录1：对分法程序（fulu1.m）syms x fx; a=0;b=1;fx=x3-3*x+1;x=(a+b)/2;k=0;ffx=s

26、ubs(fx, 'x', x);if ffx=0; disp('the root is:', num2str(x)else disp('k ak bk f(xk)')while abs(ffx)>0.0001 & a<b; disp(num2str(k), ' ', num2str(a), ' ', num2str(b), ' ', num2str(ffx) fa=subs(fx, 'x', a);ffx=subs(fx, 'x', x);if f

27、a*ffx<0 b=x; else a=x; end k=k+1;x=(a+b)/2;end disp(num2str(k), ' ', num2str(a), ' ', num2str(b), ' ', num2str(ffx)end注：实验时，可将第 2 行的 a、b 改为其它区间端点进行其它实验附录2：普通迭代法（fulu2.m）syms x fx gx; gx=(x3+1)/3;fx=x3-3*x+1; disp('k x f(x)')x=0.5;k=0; ffx=subs(fx, 'x', x);w

28、hile abs(ffx)>0.0001; disp(num2str(k), ' ', num2str(x), ' ', num2str(ffx); x=subs(gx, 'x', x);ffx=subs(fx, 'x', x);k=k+1;enddisp(num2str(k), ' ', num2str(x), ' ', num2str(ffx)附录3：收敛/发散判断（fulu3.m）syms x g1 g2 g3 dg1 dg2 dg3;x1=0.347;x2=1.53;x3=-1.88;

29、 g1=(x3+1)/3;dg1=diff(g1, 'x');g2=1/(3-x2);dg2=diff(g2, 'x');g3=(3*x-1)(1/3);dg3=diff(g3, 'x');disp('1 ', num2str(abs(subs(dg1, 'x', x1), ' ', . num2str(abs(subs(dg1, 'x', x2), ' ', num2str(abs(subs(dg1, 'x', x3)disp('2 '

30、;, num2str(abs(subs(dg2, 'x', x1), ' ', . num2str(abs(subs(dg2, 'x', x2), ' ', num2str(abs(subs(dg2, 'x', x3)disp('3 ', num2str(abs(subs(dg3, 'x', x1), ' ', . num2str(abs(subs(dg3, 'x', x2), ' ', num2str(abs(subs(dg3, &#

31、39;x', x3)附录4：松弛迭代法（fulu4.m）syms fx gx x dgx;gx=(x3+1)/3;fx=x3-3*x+1;dgx=diff(gx, 'x');x=0.5;k=0;ggx=subs(gx, 'x', x);ffx=subs(fx, 'x', x);dgxx=subs(dgx, 'x', x);disp('k x w')while abs(ffx)>0.0001; w=1/(1-dgxx); disp(num2str(k), ' ', num2str(x),

32、 ' ', num2str(w) x=(1-w)*x+w*ggx;k=k+1;ggx=subs(gx, 'x', x);ffx=subs(fx, 'x', x);dgxx=subs(dgx, 'x', x);end disp(num2str(k), ' ', num2str(x), ' ', num2str(w)附录5： Altken 迭代法（fulu5.m） syms x fx gx; gx=(x3+1)/3;fx=x3-3*x+1;disp('k x x1 x2')x=0.5;k=0;ffx=subs(fx, 'x', x);while abs(ffx)>0.0001;