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

1第七讲MATLAB学习matlab准解析法、数值方法以及迭代方法，把握对分法、迭代法、牛顿切法线求方程近似根的根本过程；把握求代数方程〔组〕的解的求解命令．教学重点：求方程近似解的几种迭代方法，代数方程〔组〕的解的求解命令的使用程〔组〕的求解命令及使用技巧．教学难点：方程的近似求解和非线性方程〔组〕的求解．一、问题背景和试验目的f(x)0的根是最常见的数学问题之一〔这里称为代数方程，主要是想和后面的微分方程区分开．为简明起见，在本试验的以下表达中，把代数方程简称为方程f(x)f(x)0为线性方程，否则称之为非线性方程．f(x)0f(x)的多样性，尚无一般的解析解法可使用，但假设能满足实际要求．同时对于多未知量非线性方程(组)而言，简洁的迭代法也是可以做出来的，但在这里我们介绍相关的命令来求解，不用迭代方法求解．通过本试验，到达下面目的：了解对分法、迭代法、牛顿切线法求方程近似根的根本过程；求代数方程〔组〕的解．首先，我们先介绍几种近似求根有关的方法：对分法到满足精度为止．对分法适用于求有根区间内的单实根或奇重实根．f(x)在[a,b上连续，f(af(b)0，即f(a)0，f(b)0f(a)0，f(b)0依据连续函数的介值定理，在(ab内至少存在一点f()0．下面的方法可以求出该根：x0

(ab/2f(x；0f(x0

)0x0

f(x)0xx．0假设f(af(x0

)0a1

a，b1

xf(af(x0

)0a1

x，b0

b；x(a1

b/2．……，有a1

、bxk

(ak

b)/2．kf(xk

) 为预先给定的精度要求)，退出计算，输出结果xk

(ak

b)/2；k反之，返回(1)，重复(1)，(2)，(3)．以上方法可得到每次缩小一半的区间序列{[abk k

]}，在(abk k

中含有方程的根．当区间长bk

axk

(ak

b/2为根的近似值，明显有kx (bk k

a)/2(bk

a

)/(22)

(ba)/2k1以上公式可用于估量对分次数k．分析以上过程不难知道，对分法的收敛速度与公比为

1的等比级数一样．由于22101024，可知大约对分10它常用来摸索实根的分布区间，或求根的近似值．迭代法f(x)0构造一个等价方程x(x)．则迭代方程是：x

(1)x(x)，1/(1”(x

))，其中”(x)1．

k k k k k k敛．Altken松弛法要先计算”(xk

，在使用中有时不便利，为此进展出以下的Altken公式：x(1)k

(x) ；kx(2)k

(x(1))；kx x(2)(x(2)x(1))2/(x(2)

2x(1)x

) k1 k k

k k kAltken5)．牛顿(Newton)法〔牛顿切线法〕f(x)0是非线性方程其迭代公式为：

x x(f(x)/f”(x)) k0,1,2,k k k牛顿法的缺点是：(1)对重根收敛很慢；(2)对初值x要求较严，要求x相当接近真值0 0x*．因此，常用其他方法确定初值x，再用牛顿法提高精度．0以下是本试验中的几个具体的试验1：对分法x33x10的实根的分布区间，再利用对分法在这些区间上分别求出根的近似值．程序如下：function[y,p]=erfenclc,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)]);holdon, whilee>10^(-5)plot([a(i),a(i)],[0,100],[x(i)x(i)],[0100],[b(i)b(i)],[0100]),pause(0.5)iff(a(i))*f(x(i))<0a(i+1)=a(i);b(i+1)=x(i);x(i+1)=(a(i+1)+b(i+1))/2;elseend

a(i+1)=x(i);b(i+1)=b(i);x(i+1)=(a(i+1)+b(i+1))/2;e=abs(f(x(i)));i=i+1;endy=x(i);p=[a;x;b]”functionu=f(x)u=x^3-3*x+1;endend图形如下:1.532132：一般迭代法承受迭代过程：x (x)求方程x33x10在0.5四周的根，准确到第4位小k数．x(x31)/3用迭代公式：x (x31)/3，k0,1,2,k具体试验3：迭代法的加速1——松弛迭代法迭代公式为

(x)(x31)/3，”(x)x2，k

1/(1x2)kxk1

(1)xk k

(xk

31)/3clc;x=[];w=[];x(1)=1;w(1)=1/(1-x(1));fori=1:10w(i)=1/(1-x(i)); x(i+1)=(1-w(i))*x(i)+w(i)*(x(i)^3+1)/3;endx另外有程序可以参考，详见参见附录4．具体试验4：迭代法的加速2——Altken迭代法迭代公式为：

x(1)

(x31)/3，x(2)

(x(1)31)/3k k k kx x(2)(x(2)x(1))2/(x(2)

2x(1)x

〕symsxfxgx;

k1 k k k

k k kgx=(x^3+1)/3;fx=x^3-3*x+1;disp(”k x x1 x=0.5;k=0;ffx=subs(fx,”x”,x);whileabs(ffx)>0.0001;u=subs(gx,”x”,x);v=subs(gx,”x”,u);disp([num2str(k),” ”,num2str(x),” ”,num2str(u),” 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]=althken%求方根的迭代程序x(1)=6;i=1;p=[];ezplot(”x^3-3*x+1”,[x(1)-9,x(1)+1]);holdonplot([x(1)-20,x(1)+2],[0,0])whileabs(f(x(i)))>=10^(-5)plot(x(i),0,”*”)4end

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)p,y=x(i),i,formatfunctionu=phi(x)u=(x^3+1)/3;endfunctionu=f(x)u=x^3+1-3*x;endend5：牛顿法x33x10在-22f(x)x33x1f”(x)3x23迭代公式：function[y,p]=newton%求方根的迭代程序

x

x2(xk

33xk

1)/(3x23)ki=1;p=[];ezplot(”x^3-3*x+1”,[x(1)-9,x(1)+1]);holdonplot([x(1)-20,x(1)+2],[0,0])whileabs(f(x(i)))>=10^(-5)plot(x(i),0,”*”), p=[p;[i,x(i)]]; i=i+1; pause(0.1)endformatshort,p,y=x(i),i,functionu=df(x)u=3*x^2-3;endfunctionu=f(x)u=x^3+1-3*x;endend结果：5结果为：1.5321※进一步思考：用迭代法求3的平方根.迭代公式为x (x3/x)/2.编写M函数n1 n nMy_sqrt.m,3x.要求误差小于105．仅要求写出源程序．试使用以上介绍的迭代法来相互比较参考程序：functiony=my_sqrt(a)%求方根的迭代程序ifnargin~=1|~isa(a,”double”), 输入数字为一个正数!”),endifa<0, error(”输入数字为正数!”),endifa>0formatlonge, x(1)=0; x(2)=1;i=1;whileabs(x(i+1)-x(i))>=10^(-5)i=i+1;x(i+1)=1/2*(x(i)+a/(x(i)+eps));endy=x(i+1);i,formatend

现在我们简洁介绍图解法如何来求解一元方程和二元方程的根：例：exp(-3*t)*sin(4*t+2)+4*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”,[05])>>holdon,line([0,5],[0,0])t=3.5203>>symsx;t=3.5203;vpa(exp(-3*t)*sin(4*t+2)+4*exp(-0.5*t)*cos(2*t)-0.5)ans=610-.43167073997540938989914138801396e-4x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)=0y^2*cos(y+x^2)+x^2*exp(x+y)=0>>ezplot(”x^2*exp(-x*y^2/2)+exp(-x/2)*sin(x*y)”)>>holdonezplot(”y^2*cos(y+x^2)+x^2*exp(x+y)”)具体的结果请大家自己下来运行〔命令〕求解方程及简介solve(”f(x)”)，f(x)为一个具体的表达式．roots(A)，Axfzero(”f(x)”,x0)，f(x)为一个具体的表达式，x0linsolve(A,b)，A，b单变量的多项式方程求根：命令格式：roots(A)例：x^3-6*(x^2)-72*x-27=0;>>p=[1-6-72-27]>>r=roots(p)r=12.1229-5.7345-0.3884多项式型方程的准解析解法命令格式：[x,…]=solve(eqn1,eqn2,…)例：x^2+y^2-1=00.75*x^3-y+0.9=0>>symsxy;>>[x,y]=solve(”x^2+y^2-1=0”,”75*x^3/100-y+9/10=0”)检验：>>[eval(”x.^2+y.^2-1”),eval(”75*x.^3/100-y+9/10”)]具体结果就请大家下来自己运行线性方程组的求解例：求线性方程组mxb的解，m=[12345;23456;345678;45678;56780]，b=[1;2;3;4;5]fori=1:5forj=1:5m(i,j)=i+j-1;endendm(5,5)=0;b=[1:5]”;linsolve(m,b)非线性方程数值求解用格式为：z=fzero(”fname”,x0,tol,trace)其中fnamex0toltol=eps，trace指10trace=0．例:f(x)=x-10x+2=0x0=0.5步骤如下：functionfx=funx(x)fx=x-10.^x+2;z=fzero(”funx”,0.5)z= 0.3758非线性方程组的求解．fsolveX=fsolve(”fun”,X0,option)X020户可以使用optimset命令将它们显示出来．假设想转变其中某个选项，则可以调用optimset函数来完成．例如，Display‘off’‘iter’‘final’optimset(‘Display’,‘off’)Display‘off’．例:求以下非线性方程组在(0.5,0.5)functionq=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);fsolvex=fsolve(”myfun”,[0.5,0.5]”,optimset(”Display”,”off”))x= 0.63540.3734q=myfun(x)q= 1.0e-009*0.2375 0.2957可见得到了较高精度的结果．精品案例：螺旋线与平面的交点问题 ：螺旋线与平面相交的状况多种多样,依据螺旋线与平面方程的不同可以相交,也可以不相交.在相交的状况下,可以交于一点,也可以交于好多点.对于各种相交的状况,要求其交点的坐标并不是一件简洁的事.本次试验就以此为背景争论下面的具体问题：螺旋线的参数方程为x4cos,y4sin,z,08.xy0.5z20求该螺旋线与平面的交点.要求：求出全部交点的坐标；在同一图形窗口画出螺旋线、平面和交点.试验过程：问题分析fsolve等.先对方程化简,削减变量个数,使用图解方法求方程的根.再分别画出螺旋线,平面,及其交点.算法描述与分析fsolve,选择适当的初值,求其数值解;再分别会出图形；最终对图形作出必要的修饰.源程序及注释cos4sin0.520cos4sin20.5建立下面M文件intersect_point.m%使用图解法求交点，并且三维图%画图确定解的个数和或许位置theta=0:0.01:8*pi;y1=4*(cos(theta)+sin(theta));y2=2-0.5*theta;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) holdon%透亮的画平面%x1=-5:0.1:5;[-4,4]有关.y1=x1;[X1Y1]=meshgrid(x1,y1);%网格化，画曲面Z1=4-2*X1-2*Y1;surf(X1,Y1,Z1) mesh(X1,Y1,Z1)shadingflatalpha(0.5) %设置透亮度alpha(”z”) %设置透亮度方向ttheta%i=1forn=[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”)测试结果〔写清输入输出状况〕交点坐标为：

08

内与三角曲线有4个交点.theta的数值解为：t=[2.1961 5.3759 9.1078 11.1023]四个交点的近似坐标为：x0=[-2.3413 2.4635 -3.8007 0.4261]y0=[3.2432-3.15141.2468-3.9772]z0=[2.19615.37599.107811.1023]调试和运行程序过程中产生的问题及实行的措施求交点的时候会消灭重根和漏根的情形,通过选择适当的初值避开了上述状况.对算法和程序的争论、分析,改进设想及其它阅历教训solve函数只能求解一个数值解,不能全部求出;fsolve函数好;为了满足更好的视觉1112效果,可以对图形进展进一步的修饰.习题1.f(x)3x5x42x3x23sin(xyyex0(1)x2y2 (2)求解方程:cos(x)xex求解多项式方程x9810求以下代数方程〔组〕的解：(1) x5x10(2) 2x3y0 ①4x23y1 ②〔1一般迭代法〔2与之相应的松弛迭代法和x33x10在1.4四周的根，准确到4变化．Altken迭代法、牛顿切法线等5方程txsin(x)的正的近似根，0t1．(建议取t0.5．时间许可的话，可进一步考虑t0.25的状况．)五、附录为供近似求根的算法1：对分法程序〔fulu1.m〕symsxfx;a=0;b=1;fx=x^3-3*x+1;x=(a+b)/2;k=0;ffx=subs(fx,”x”,x);ifffx==0;disp([”therootis:”,num2str(x)])elsedisp(”kakbkf(xk)”)whileabs(ffx)>0.0001&a<b;num2str(a),””,num2str(b),””,num2str(ffx)])fa=subs(fx,”x”,a);ffx=subs(fx,”x”,x);iffa*ffx<0b=x;elsea=x;endk=k+1;x=(a+b)/2;enddisp([num2str(k),” ”,num2str(a),” ”,num2str(b),” ”,num2str(ffx)])end注：试验时，可将第2行的a、b改为其它区间端点进展其它试验．2：一般迭代法〔fulu2.m〕symsxfxgx;gx=(x^3+1)/3;fx=x^3-3*x+1;disp(”k x x=0.5;k=0;ffx=subs(fx,”x”,x);whileabs(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〕symsxg1g2g3dg1dg2dg3;x1=0.347;x2=1.53;x3=-1.88;g1=(x^3+1)/3;dg1=diff(g1,”x”);g2=1/(3-x^2);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 ”,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,”x”,x3)))])4：松弛迭代法〔fulu4.m〕symsfxgxxdgx;gx=(x^3+1)/3;fx=x^3-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”)whileabs(ffx)>0.0001;w=1/(1-dgxx); disp([num2str(k),” ”,num2str(x),” ”,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);enddisp([num2str(k),” ”,num2str(x),” ”,nu