数学实验非线性方程求解

1、大学数学实验作业非线性方程求解班级： 姓名： 学号： 日期： 目录目录2【实验目的】3【实验内容】3【题目1】（课本习题第六章第1题）31.1模型分析31.2模型求解31.3不同初值给出的结果51.4根收敛域的计算71.5构造迭代公式与牛顿法求解81.6结果分析与讨论10【题目2】（课本习题第六章第3题）102.1模型建立102.2第(1)问求解112.3第(2)问求解152.4本题小结20【题目3】（课本习题第六章第7题）203.1初步观测迭代序列变化的曲线203.2模型分析213.3模型求解223.4本题小结27【实验心得、体会】27注：本实验作业脚本文件均以ex6_1_1形式命名，其中e

2、x代表作业，6_1_1表示第六章第一小题第一个程序。本次试验编的脚本文件较多，word中每个脚本均有对应M文件。自编函数均以exf6_1_1形式命名，exf代表作业函数，6_1_1表示第六章第一题第一个自编函数。特殊函数，如混沌函数chaos 、分岔点函数fenchadian2、iter函数除外。【实验目的】1、掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作出初步的分析；2、练习用非线性方程组建立实际问题的模型并进行求解。【实验内容】【题目1】（课本习题第六章第1题）分别用fzero和fsolve程序求方程的所有根，准确到，取不同的初值计算，输出初值、根的近似值和迭代次数，

3、分析不同根的收敛域；自己构造某个迭代公式（如等）用迭代法求解，并自己编写牛顿法的程序进行求解和比较。1.1模型分析fzero命令主要用于单变量方程的求根，主要采用二分法、割线法和逆二次插值法等的混合方法。fzero至少需要两个输入参数：函数、迭代初值（或有根区间）。fsolve命令主要用于非线性方程组的求解，可以输出结果（如x点对应的雅可比矩阵等）。已知y=sinx是一个周期性有界函数，而二次函数在对称轴两边增长无界。我们可以直接观察出x=0是方程的解，再作出该方程两边所代表的函数的图像。从图上可以观察到另外一个根的范围。而由两函数性质，在第二根之外，二次函数增长，而三角函数波动，再也不会有交

4、点。从而大致可知此方程只有两解。1.2模型求解编写待求解函数M文件：% -作业题6_1待求解求解函数程序exf6_1-function y = exf6_1 (x)y=sin(x)-x2/2;end1.2.1fzero求解因另外一个解范围未知，故先编一段程序估计解的范围。%-作业题6_1估计解的范围源程序ex6_1_1-clear;clc;clf;x=-2:0.1:2; %估计根的范围y1=sin(x);y2=x.2/2; %分别先作出两函数图像plot(x,y1,x,y2); %图像交点处即为根，观察出解的数目与分布xlabel(x);ylabel(y1/y2);title(y1、y2函数图

5、像); % 加入X轴、Y轴标记和标题legend( y1=sin(x), y2=x.2/2); % 加入图例得到如下图像，已知一个解为0，由图像可得另一解在1,2之间。图1：y1=sin(x), y2=x.2/2交点图根据两解分布，x1-1,1，x21,2，编写matlab：%-作业题6_1用fzero求方程的所有根源程序ex6_1_2-clear;clc;opt=optimset(fzero);opt=optimset(opt,tolx,1e-10); %opt设定求解精度format long g;x,fv,ef,out=fzero(exf6_1,1,2,opt)x,fv,ef,out=f

6、zero(exf6_1,-1,1,opt) %根据解分布，用fzero函数求解输出结果如下:x = 1.40441482402454fv = 8.41122727024413e-011ef = 1out = intervaliterations: 0 iterations: 7 funcCount: 9 algorithm: bisection, interpolation message: Zero found in the interval 1, 2另一根为：x = 1.74713912083679e-011fv = 1.74713912082153e-011ef = 1out = int

7、ervaliterations: 0 iterations: 8 funcCount: 10 algorithm: bisection, interpolation message: Zero found in the interval -1, 1从而可知，方程一根为x = 1.40441482402454，另一根即为x = 0。1.2.2用fslove求解%-作业题6_1用fsolve求方程的所有根源程序ex6_1_3-clear;clc;opt=optimset(fzero);opt=optimset(opt,tolx,1e-10); %opt设定求解精度format long g;x,f

8、v,ef,out=fsolve(exf6_1,1,opt)x,fv,ef,out=fsolve(exf6_1,0,opt) %根据解分布，用fsolve函数求解输出结果如下：x = 1.40441482411066fv = -2.25777174733821e-011ef = 1out = iterations: 6 funcCount: 13 algorithm: trust-region-dogleg firstorderopt: 2.79692747209721e-011 message: 1x76 charjac = -1.23879992536652另一根为：x = 0fv = 0e

9、f = 1out = iterations: 0 funcCount: 2 algorithm: trust-region-dogleg firstorderopt: 0message: 1x76 char表1：fzero与fsolve求解结果根fzero求解fsolve求解x11.404414824024541.40441482411066x21.74713912083679e-011x = 01.3不同初值给出的结果1.3.1变化fzero的初值fzero主要是利用二分法、割线法等的混合来求解，因此x的初值对其计算出数值解是有影响的。下面的列表给出了不同x的初值，fzero函数输出的不同结

10、果。根据根的范围，x取值从-2变化到3.5。表2：不同x的初值，fzero函数输出的不同结果x初值根的近似值迭代次数函数调用次数-25.824116516731E-13831-1.58.278977421518E-12731-1-8.948271056829E-11630-0.5-1.904088574048E-1563000010.51.619924113706E-1673111.404414824090276251.51.4044148240919851321.404414824092435212.51.4044148240924352331.404414824092428283.51.4

11、0441482405939626由表2可得，x在0附近取值时，根的近似值都为0，迭代次数在7次左右，函数调用次数约为30次；x在1.4附近取值，根的近似值都为1.404414，迭代次数在6次左右。发现x取1.5时函数调用次数较少，说明离根很近；x=0时，直接取到根，因此迭代次数与函数调用次数均最少。1.3.2变化fsolve的初值根据根的范围，x取值从-2变化到3.5。表3：不同x的初值，fsolve函数输出的不同结果x初值根的近似值迭代次数函数调用次数-2-2.6071879498599e-010512-1.5-1.04384412156096e-012512-1-2.60718794985

12、99e-010410-0.5-1.0438441213492e-01241000020.5-8.02139615750056e-01361311.404414824110666131.51.404414824927073821.404414840973194102.51.4044148241100351231.404414840973195123.51.40441482411003614由表3可得，x在0附近取值时，根的近似值都为0，迭代次数在5次左右，函数调用次数约为12次；x在1.4附近取值，根的近似值都为1.404414，迭代次数在5次左右。发现x取1.5时函数调用次数较少，说明离根很近

13、；x=0时，直接取到根，因此迭代次数与函数调用次数均最少。1.4根收敛域的计算有了以上不同x初值时根的近似值结果，进一步计算函数收敛域1.4.1fzero的收敛域首先进行简单的试探，大致得出如下结论：函数有两个根，根x=0的收敛域形式为(-,a)，大致确定a0,1根x=1.4044的收敛域形式为(a,b)，大致确定b14,15当x0b时，fzero解法不收敛，将会解出x=NaN为了确定ab两个数，编写两条循环语句如下:%-作业题6_1根收敛域的计算源程序ex6_1_4-clear all;clc;a=0.7; %设定a的初值opt=optimset(fzero);opt=optimset(op

14、t,tolx,1e-10); %opt设定求解精度while(x1) a=a+0.00001; x=fzero(exf6_1,a,opt);enda %输出函数值趋向无穷大时的ab=14; %设定b的初值x=1;while(xtol)&(sin(x(i)eps) x(i+1)=(2*sin(x(i)0.5; u=x(i+1)-x(i); i=i+1; if(in) error(n is full); endendxi-1 %输出迭代结果及迭代次数输出结果如下：ans = 1 1.29728253268738 1.38767986886849 1.40234159737526 1.4041688

15、093621 1.40438579137907 1.40441140012416 1.40441442031852 1.40441477647754 1.40441481847747 1.40441482343029 1.40441482401435 1.40441482408323ans = 12迭代次数为12。对该方法迭代次数进行统计：表6：用 x=(2sinx)1/2求解函数不同初值的迭代次数用 x=(2sinx)1/2求解初值根的近似值迭代次数11.404414824083231221.404414824088321231.404414824083531471.404414824083

16、531381.404414824083531091.4044148240907314131.404414824090814141.404414824090811当初值取其他值时，x=(2 sinx )(1/2)无意义。该迭代方法收敛次数多于fzero法和fsolve法。1.5.2牛顿法求解牛顿法程序代码：%-作业题6_1牛顿法求解源程序ex6_1_6-clear;clc;x0=1;x(1)=x0; %给定迭代初值tol=1e-6; %给定误差极限u=1;n=10; %迭代计数，i大于10时跳出i=1;while(abs(u)tol) x(i+1)=x(i)-(sin(x(i)-x(i)2/2)

17、/(cos(x(i)-x(i); u=x(i+1)-x(i); i=i+1; if(in) error(n is full); endendxi-1 %输出迭代结果及迭代次数迭代次数为6。对牛顿法迭代次数进行统计：表7：用牛顿法求解函数不同初值的迭代次数用牛顿法求解初值根的近似值迭代次数-206-10600111.40441482409243621.40441482409243631.40441482409243741.40441482409243751.40441482409243761.40441482409243771.40441482409243881.4044148240924389

18、1.404414824092438101.404414824092438111.404414824092438121.404414824092438131.404414824092438141.404414824092439牛顿法迭代次数随初值渐渐远离根的真实值而逐渐增大，增大速度递减。【公式迭代与牛顿法比较】公式法：由于选取的迭代公式在某些初值设定下，会使表达式无意义，故该公式迭代只给出部分值的迭代结果。另外，在尝试中，我发现当表达式在实数范围内无意义，给出复数解时，给出的实部总是1.4044这个根。如果给定的初值不为0，就永远不可能迭代出0这个根。这是由于的图形在1.4044附近斜率是小于

19、1的，根据局部收敛性定理，它在1.4044附近连续可微，导数绝对值小于1，因此会收敛到根。但是在0附近不满足不动点收敛的条件，是一个蛛网型的迭代。故不论以其他任何值迭代，都不会迭代到0。牛顿法：比较可知，等精度的情况下牛顿法的迭代次数较少，是一种更有效的算法。而且牛顿法适当改变初值就可以得到两个根，而迭代法由于局部收敛性的限制，则只能得到一个根。总体来说不如牛顿法。1.6结果分析与讨论以上四种迭代方法（fzero、fsolve、x=(2 sinx )(1/2)）、牛顿），均可得到根x2的近似值。其中，自己构造的迭代公式无法得到根x1的值。四种迭代方法中fsolve和牛顿法迭代次数较少，收敛较快

20、。从本质上讲，牛顿法只是一种特殊的迭代法，它在迭代函数的选择上有讲究。相对于其他迭代方法，也许牛顿法不是最有效的，但是它是一种有具体表达形式的迭代法。但是它利用导数构造公式，可能让运算稍慢。【题目2】（课本习题第六章第3题）(1)小张夫妇以按揭方式贷款买了一套价值20万的房子，首付了5万元，每月还款1000元，15年还清。问贷款利率是多少？(2)某人欲贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。从利率方面看，那家银行较优惠（简单的假设年利率=月利率*12）？2.1模型建立设yk为第k个月的欠款数，

21、设月利率为r，房价为c万元，首付为d万元，每月还款b万元，则有 整理得根据递推关系知即根据题意，c=20，d=5，b=0.1，y180=0，求解月利率r。上式即为按揭方式贷款购房（按月计算）的数学模型。2.2第(1)问求解按照已经建立好的按揭方式贷款购房模型，将第一问中的条件数据（房价c=20，首付d=5，每月还款数b=0.1，还款月数k=180）代入，可得2.2.1求根的大致范围：令利用MATLAB输出f(r)的图像。%-作业题6_3第一问图解根的范围ex6_3_1-clear all; clc;n=30; for i=1:n r(i)=0.0001*i;b=r(i); f(i)=150*b

22、*(1+b).180; g(i)=(1+b).180-1; c=f(i);d=g(i); h(i)=c-d; end; plot(r,h),hold on,grid on;xlabel(r);ylabel(f(r）);title(f(r)图像); % 加入X轴、Y轴标记和标题图2：函数零点图得到上图，从上图可作出初步估计，r在0.00200和0.00225之间。2.2.2下面利用二分法、牛顿法、fzero函数三种不同的方法求解：利用二分法求解下面的程序实现的是用二分法求解方程，需要说明几点： 1、r1和r2分别表示之前估计的r的范围的上下限。 2、由于在这个实际问题中，不可能出现某个r3=0.

23、5(r1+ r2)正好是非线性方程（1）的解，故在编写程序的时候未将这种情况考虑进去，而是假设f(r1)f(r3)和f(r2)f(r3)总有一个大于零另一个小于零。 3、二分的次数达到100次时，认为已经求出足够精确的解。%-作业题6_3第一问二分法求解源程序ex6_3_2-clear all; clc;m=200000;a=50000;n=15*12;x=1000;p=m-a;i=0; r1=0.002;r2=0.00225; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; end; r=(r1+r2)/

24、2;format long g;r输出的结果如下：r = 0.00208116388945975 0.2081%利用牛顿法求解 已知 那么对r求导得 故有 根据迭代公式 可求出非线性方程（1）的解。先利用MATLAB输出 的曲线，观察rk是否收敛。%-作业题6_3第一问牛顿法观察图解观察根ex6_3_3-clear all; clc;n=15; for i=1:n r(i)=0.0001*(i+14);b=r(i); c=(1+b).180;d=(1+b).179; c1=b*c;d1=b*d; f(i)=b-(150*c1-c+1)/(27000*d1+150*c-180*d); g(i)=

25、b; end; plot(r,f,r,g),hold on,grid on;xlabel(r);ylabel(y=r/y=(r);title( y=r/y=(r)图像); % 加入X轴、Y轴标记和标题legend(y=(r) , y=r)图3：牛顿法求解收敛性判断图从上图可以看出，蓝色曲线在交点处的斜率的绝对值显然小于1，说明序列rk收敛。下面的程序实现的是用牛顿法求解方程,需要说明几点： 1、取初值r0=0.0015。 2、根据上图可知，当取上条中指定的初值时，序列rk收敛，可得到方程的解。 3、取迭代次数为100次，认为此时得到的r已经足够精确。%-作业题6_3第一问牛顿法求解源程序ex6

26、_3_4-clear;clc;r=0.0015;m=0; while (m=100) b=r; c=(1+b).180; d=(1+b).179; c1=b*c;d1=b*d; a=b-(150*c1-c+1)/(27000*d1+150*c-180*d); r=a;m=m+1; end format long g;r输出的结果为: r = 0.00208116388945918 0.2081%利用fzero求解fzero是MATLAB自带的求解单变量非线性方程的命令，所采用的算法主要是二分法、割线法和逆二次插值法等的混合方法。 需要说明的是：fzero求解可以不事先知道解的范围。在Matla

27、b中编写程序如下：%-作业题6_3按揭贷款购房模型函数M文件源程序exf6_3-% b:每月还款数，c:房价，d:首付数，n:还款年数，x:贷款利率function y=exf6_3(x,c,d,k,b) y=(c-d)-b/x)*(1+x)k+b/x; % 按照差分方程的解构造函数 end%-作业题6_3第一问求解脚本M文件源程序ex6_3_5-clear all;clc;c=20;d=5;b=0.1;k=180; % 代入条件：房价c=20，首付d=5，每月还款数b=0.1，还款月数k=180x0=1; % 设定迭代初值x,fv,ef,out=fzero(exf6_3,x0,c,d,k,b

28、) % 利用fzero命令求解贷款利率得到的运行结果如下：x = 0.00208116388945936 % x点对应的函数值fv = -2.62900812231237e-013ef = 1 % 程序停止时的状态，正数（1）表示找到异号点out = intervaliterations: 12iterations: 22 % 迭代了22次funcCount: 46 % 函数被调用了46次algorithm: bisection, interpolation % 算法为二分法和插值法message: Zero found in the interval -0.28, 1.9051输出的结果为:

29、 r = 0.00208116388945936 0.2081%表8：二分法、牛顿法、fzero求解结果二分法求出的解 牛顿法求出的解 利用fzero求出的解 0.00208116388945975 0.00208116388945918 0.00208116388945936故小张夫妇的贷款的银行的月利率是0.2081%。2.3第(2)问求解2.3.1解取值范围判断根据模型，第一家银行的月利率r1满足下述方程：第二家银行的年利率r2满足下述方程：设利用MATLAB输出f1(r)和f2(r)的图像。%-作业题6_3第二问求解脚本M文件源程序ex6_3_6-clear all;clc; n=90

30、;m=800; for i=1:n r1(i)=0.0001*i;b=r1(i); f1(i)=1000*b*(1+b).180;g1(i)=9*(1+b).180-9; c1=f1(i);d1=g1(i);h1(i)=c1-d1; end; for i=1:m r2(i)=0.0001*i;b=r2(i); f2(i)=100*b*(1+b).20;g2(i)=9*(1+b).20-9; c2=f2(i);d2=g2(i);h2(i)=c2-d2; end; subplot(1,2,1),plot(r1,h1), grid on; xlabel(r);ylabel(y=f1(r);title

31、(y=f1(r)图像); % 加入X轴、Y轴标记和标题subplot(1,2,2),plot(r2,h2), grid on;xlabel(r);ylabel(y=f2(r);title( y=f2(r)图像); % 加入X轴、Y轴标记和标题图4：解取值范围判断图观察上两图，可知r1的范围在0.0055到0.0060之间，r2的范围在0.060到0.065之间。2.3.2下面分别用二分法、牛顿法、fzero函数三种方法求解 二分法求解下面的两个程序实现的分别是利用二分法的原理求解方程。之前的说明仍然有效。第一家银行：%-作业题6_3第二问二分法求解第一家银行利率ex6_3_7-clear al

32、l; clc;i=0;r1=0.0055;r2=0.0060; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; end; r1=(r1+r2)*6 输出的结果为年利率 :r1= 0.07020951099414317.021%第二家银行%-作业题6_3第二问二分法求解第二家银行利率ex6_3_8-clear all; clc;i=0;r1=0.060;r2=0.065; while (i0 & f2*f30) r1=r3; else if (f1*f30) r2=r3; end; end; i=i+1; e

33、nd; r2=(r1+r2)/2 输出的结果为 r2=0.06394877709238636.395% 牛顿法求解已知 那么对r求导得 故有 根据迭代公式 rk+1=1,2(rk)可求出非线性方程的解。 先利用MATLAB输出 的曲线，观察rk是否收敛%-作业题6_3第二问牛顿收敛性判断程序ex6_3_9-clear;clc;n=15; m=15; for i=1:n r1(i)=0.0001*(i+50);b=r1(i);c=(1+b).180;d=(1+b).179;c1=b*c;d1=b*d; f1(i)=b-(1000*c1-9*c+9)/(180000*d1+1000*c-1620*

34、d);g1(i)=b; end; subplot(1,2,1),plot(r1,f1,r1,g1),grid on;xlabel(r);ylabel(y=1(r)/y=r);title( y=1(r)/y=r图像); % 加入X轴、Y轴标记和标题legend(y=1(r) ,y=r);for i=1:m r2(i)=0.001*(i+54);b=r2(i);c=(1+b).20;d=(1+b).19;c1=b*c;d1=b*d; f2(i)=b-(100*c1-9*c+9)/(2000*d1+100*c-180*d);g2(i)=b; end; subplot(1,2,2),plot(r2,f

35、2,r2,g2), grid on;xlabel(r);ylabel(y=2(r)/y=r);title( y=2(r)/y=r图像); % 加入X轴、Y轴标记和标题legend(y=2(r) ,y=r);得到以下图像：图5：牛顿法求解收敛性判断图从上面两幅图可以看出，两条蓝色曲线在交点处的斜率的绝对值显然都小于1，说明对r1和r2序列rk都收敛。下面的程序实现的是用牛顿法求解第一家银行利率。 取初值为0.0052，其他说明同之前的类似。%-作业题6_3第二问牛顿法求解第一家银行年利率ex6_3_10-clear;clc;r=0.0052;m=0; while (m=100) b=r; c=(

36、1+b).180; d=(1+b).179; c1=b*c;d1=b*d; a=b-(1000*c1-9*c+9)/(180000*d1+1000*c-1620*d); r=a;m=m+1; end r1=12*r输出的结果为年利率 ：r1= 0.07020951099414477.021%下面的程序实现的是用牛顿法求解方程第二家银行利率。 取初值为0.056，其他说明同之前的类似。%-作业题6_3第二问牛顿法求解第二家银行利率ex6_3_11-clear;clc;r=0.056;m=0; while (m0)任意，初值观察是否有混沌现象出现，并找出前几个分岔点，观察分岔点的极限趋势是否符合F

37、eigenbaum常数揭示的规律。3.1初步观测迭代序列变化的曲线%-作业题6_7输出xk随着迭代次数变化的曲线源程序ex6_7_1-clear;clc;a=5,11,15;b=2; i=1;c=1; j=1;x=0:1:25; for j=1:3for i=1:26 x(i)=c;p=c*a(j)*exp(-b*c); c=p; k(i)=i-1; end subplot(1,3,j), plot(k,x), grid ;xlabel(k);ylabel(Xk);title( 序列xk迭代图像); % 加入X轴、Y轴标记和标题enda=5迭代图像 a=11迭代图像 a=15迭代图像图6：a分

38、别取5、11、15，b=2时函数迭代图像观察上面三幅图片，可知： a=5时，xk收敛于1个值，这个值是0.8047； a=11时，xk有两个收敛的子序列，呈周期2收敛，分别趋向0.4108和1.9871； a=15时，xk没有收敛的子序列，变化没有规律。3.2模型分析非线性差分方程参数不同时，迭代结果会呈现不同的走向。有可能收敛，有可能存在多个收敛的子列，甚至可能完全不收敛。本题通过计算求解一个迭代公式的趋势，以研究混沌、分岔现象。选择函数时，应考虑其平衡点问题。本题由方程确定。此方程有两个解x=0与，令x1，也就是lnaln15，在以下的计算中，取b=5。3.3模型求解3.3.1混沌分岔图：

39、%-作业题6_7混沌分岔图函数M文件源程序chaos-function chaos(iter_fun,a,b,n) %该函数没有返回值；iter_fun是迭代函数（句柄）x0=1; %迭代初值kr=0;for aa = a(1):a(3):a(2) %输入中a(1),a(2)是参数变化的范围，a(3)是步长 kr = kr+1; y(kr,1)=feval(iter_fun,x0,aa,b); for i=2:n(2) %输入中n(2)是迭代序列的长度，但画图时前n(1)个迭代值被舍弃 y(kr,i)=feval(iter_fun,y(kr,i-1),aa,b); endendplot(a(1