数学实验报告-利用MALTAB计算非线性方程近似解

1、实验四 非线性方程近似解一、按揭还贷问题描述（1）小张夫妇以按揭方式贷款买了一套价值20万元的房子，首付5万元，每月还款1000元，15年还清。问贷款利率是多少？（2）某人想贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。从利率方面看，哪家银行较优惠？（简单假设年利率=月利率*12）简要分析初看本题，一个简单的思路是每次测试一个利率值，以这个值为基础计算15年后所剩还款数量，通过结果判断应将利率值增大或减小，从而实现迭代。这其实是一个二重迭代的过程，之所以这样是因为不容易一眼看出本题的非线性方程。事实

2、上，转换思路后，可以利用一个简单的方程描述整个迭代过程。这样就将二重迭代转化为了一层迭代。使得处理更加简便。方法与公式1、解题方法(1)二次迭代给定总的本金，从每一次还款中扣去这段时间中增加的利息，再将其还到本金，使本金总量逐渐减少。代码： for i = 1:time*12 less = (repay-left*interest); left = left - less;(2)方程描述虽然并不是所有本金都在还款的整个期间中产生了相应的利息，但是可以设想成这样，与此同时，还款从在相应的还款时间开始产生利息，这样可以得出，两者最终的“本息和”相等，即A(1+q)n=i=1nP(1+q)n-i其中

3、A为总还款金额，q为了利率，P为每次还款金额。2、解方程方法(1)牛顿法xk+1=xk-f(xk)f'(xk)(2)直接使用公式fzero()结果与分析1、第一问：(1)二次迭代i,q=iterate(150000,1000,15,2,0,1,100,10-6);公式表意为：总贷款量=200000-50000=150000； 每月还款100元； 还款期限15年；还款方式为按月还款；迭代区间设定为0,1；最大迭代次数为100次；精度要求为10-6；最终结果为：迭代次数：45；使用时间0.5705s；利率为0.9457。(2)公式A(1+q)n=i=1nP(1+q)n-i=P(1+q)n-

4、1q由于函数中出现q在分母位置的情况，故不能将迭代区间设为0,1，为了能够更好地比较各种方法，因此，这里设置一端为1，一端尽量小。经测试，当设置值比10-15还小时，matlab自动认为其为0，故取10-15。所得结果如下:迭代次数：20使用时间不稳定，最长使用0.0066s，最短使用0.0014s。利率为0.9460(3)使用自己编写的牛顿迭代公式按照书中给出的牛顿迭代公式，进行迭代。迭代次数随初值的设定而不断变化，初值的设置要求更加复杂，不能达到10-15这样小。这里仅实验了10-6的情况，事实上，此时的迭代次数已经相当大，用时也不短。具体数据如下：迭代次数：33703使用时间：1.299

5、6s利率为0.8815 (4)简要小结从以上三种计算方法可以看出，各种计算方法的结果均相差不多，可见各种方法的正确性都是有保障的。与此同时不同的方法有着自己不同的特性：a.由于涉及公式的使用，第二、三种方法都出现了初值设置的问题。当然这个问题视具体情况而定，不能一概而论；b.本次实验中发现，前两种办法在给定区域较广的情况下能够较快地收敛到所求的位置，而牛顿法使用的时间反而较多。为了进一步确认，我将牛顿法的迭代初值置为1，此时共迭代395次，使用时间为0.1168。可以发现仍旧较大。可见自己编写的迭代方法距matlab所用的混合方法相差较远。2、第二问 (1)月付款仍旧使用三种方法并进行比较二次

6、迭代fzero自编牛顿法用时0.13 0.884 0.263 迭代次数47 23 170结果0.2844 0.2845 0.286 （2）年付款用时0.920 0.590 0.018 迭代次数442353结果0.241 0.239 0.245 从实验结果可以看出，各种方法所得的最终数据基本一致。由题目假设：年利率=月利率*12得到第一中按揭方式的年利率约为0.070209511，比第二种方法大。故第二家银行开出的条件较优惠。3、总结a.二分迭代和牛顿迭代法通常都能通过迭代收敛得到所要的结果；b.使用Matlab自带的fzero通常能够更快地得到结果，相比之下牛顿迭代法所用时间较长；程序清单1、

7、二重迭代(1)计算剩余本金的函数function left = Return(All,repay,interest)All = 200000-50000;repay = 1000;left = All;for i = 1:15*12 less = (repay-left*interest); left = left - less;end(2)主迭代函数function i,q = iterate (All,repay,time,YearOrMonth,a,b,n,tol)interest = (a+b)/2;i=0;diff = Inf;before = Inf;while(i<n &a

8、mp;& diff>tol) left = Return(All,repay,time,YearOrMonth,interest); if(left>0) b=interest; interest = (a+interest)/2; else a=interest; interest = (b+interest)/2; end i = i+1; diff = abs(before-left); before = left;endq=interest; 2、fzero(1)函数function y = Function(x)A=500000;P=45000;n=20;y =

9、A*(1+x)n-P*(1+x)n-1)/x;(2)调用方法tic;q,fv,ef,out = fzero(Function,1*10-15,1)time = toc;3、自编牛顿迭代法 (1)函数的导数function y = DeltaFunction(x)A=500000;P=45000;n=20;y = A*n*(1+x)(n-1)-P*(1+x)(n-2)*(n-2)*x-1)/x2;(2)迭代主方法x0 =1; n = 500000; tol = 1e-6;x = x0;i=1;diff = inf;tic;while(abs(diff)>tol && i&l

10、t;n) Fxk = Function(x); FDxk = DeltaFunction(x); x = x - Fxk/FDxk; diff = Function(x)-Fxk; i = i+1;endtime = toc;二、气缸阀门问题描述由气缸控制的门关闭状态如图所示。门宽为a，门枢在H处，与H相距b处有一门销，通过活塞与圆柱形气缸相连。气缸长l0，缸内气体压强为p0。用力F推开门时，气缸绝热压缩，压强增大。已知绝热方程为pV=c,试通过力矩平衡求在一定力作用下开门角度。方法与公式1、MATLAB自带fzero()函数x,fv,ef,out = fzero(valve,0);2、牛顿迭

11、代法xk+1=xk-f(xk)f'(xk)3、二分迭代法4、建模F*a=P*S*bS=pi*r2P*(l0-x)=P0*l0综上，有：y=F*a/b-pi*r2*(l0/(l0-x)Gamma*p0;结果与分析1、fzero()取初值为0x,fv,ef,out = fzero(valve,0);得到结果：c = 0.141232770027359alpha = 0.514234297016288 = 29.463454899910120用时0.6654s2、牛顿迭代法得到结果：c = 0.141232754515091alpha = 0.514234249979021 = 29.463

12、452204873240用时0.8827s3、二分迭代法得到结果：c = 0.141232490539551alpha = 0.514233449535557 = 29.463406342841026用时0.1362s4、总结 这道题比较简单，计算过程中没有涉及过多的问题。不过有关建模，有一个问题值得讨论，即阀门收到的大气压强。事实上，大气压强大小约为105，远大于题目中阀门的内压强，这样算来题目所给参数有误。所以只好认为题目所给参数已经考虑了大气压。有关各种方法的比较如下：a.三种方法所得结果相差无几，可见结果的正确性；b.时间使用上与第一题的实验结果基本一致使用Matlab自带fzero和

13、二分迭代两种方法所用时间相差不大，两者都好于自己编写的牛顿迭代法。程序清单1、fzero()(1)函数function y=valve(x)a=0.8;b=0.25;r=0.04;l0=0.5;p0=104;Gamma = 1.4;F=25;y=F*a/b-pi*r2*(l0/(l0-x)Gamma*p0;(2)调用脚本tic;x,fv,ef,out = fzero(valve,0);time = toc;b=0.25;alpha = atan(x/b);angel = alpha/pi*180;2、牛顿迭代法(1)函数与fzero方法一致(2)函数的导数function y=Deltaval

14、ve(x)a=0.8;b=0.25;r=0.04;l0=0.5;p0=104;Gamma = 1.4;F=25;y=-9*pi*r2*(l0/(l0-x)Gamma*p0/(l0-x);(3)迭代主方法x0 =0; n = 500000; tol = 1e-6;x = x0;b=0.25;i=1;diff = inf;tic;while(abs(diff)>tol && i<n) Fxk = valve(x); FDxk = Deltavalve(x); x = x - Fxk/FDxk; diff = valve(x)-Fxk; i = i+1;endtime =

15、 toc;alpha = atan(x/b);angel = alpha/pi*180;3、二分迭代法a=0.8;b=0.25;r=0.04;l0=0.5;p0=104;Gamma = 1.4;F=25;n=500;tol = 1e-6;x0 = 0;xx = l0;x = (x0+xx)/2;i=0;diff = Inf;before = Inf;tic;while(i<n && diff>tol) deltax = valve(x); if(deltax<0) xx=x; x = (x0+x)/2; else x0=x; x = (xx+x)/2; end

16、 i = i+1; diff = abs(before-x); before = x;endtime = toc;alpha = atan(x/b);angel = alpha/pi*180;三、Feigenbaum问题描述给定一迭代公式，计算序列，分析其收敛性。给定不同的参数值，观察是否由混沌现象出现，并找到前几个分叉点，观察分叉点的极限趋势是否符合Feigenbaum常数解释的规律。方法与公式1、为判断相应a的收敛性，对不同的a值进行迭代（给定一个b值），并记录迭代结果，做出图像，从而得出结论。2、为说明b的任意性，借助书中所给的有关混沌的示例，将其进行改写，并重写迭代函数，得到如下迭代过

17、程for bb = b(1):b(3):b(2) kr = kr+1; y(kr,1)=feval(iter_fun,x0,a,bb); for i=2:n(2) y(kr,i)=feval(iter_fun,y(kr,i-1),a,bb); endend3、迭代函数xk+1=a*xk*exp(-b*xk)4、为判断是否有混沌现象以及确定分叉点，将上述二重迭代过程的对象由b转为a，即给定一个b值，观察a的变化对迭代终值的影响即可。结果与分析1、a固定(1)a=5有图可以看出，当迭代次数足够大时，值趋于稳定，收敛到一点，这一点的值约为0.3218876。从图中可以看出，a=5时，对任意b>

18、0，xk收敛，且由图中趋势以及函数的形式可以看出，b增大后，xk将仍旧收敛，而且收敛值逐渐减小，趋近于0。(2)a=11由图可以看出，当a=11时，迭代不再收敛到一个点，而是有两个收敛的子序列，分别趋于0.794826376和0.164331679。从图中可以看出，当a=11时，对任意b，xk均收敛到两个子列中，且两个子列的极限值随着b的增加不断减小，趋于0。(3)a=15 实验结果表示，我们已经已经无法从图中准确地判断迭代是否收敛。大体上看并不收敛。至于具体情况，需要进一步分析。实验结果表示，当a=15时，对任意b，已经不存在确定收敛子列。2、观察分叉与混沌 将a设定在一个范围内变动，固定b

19、的值。从图中可以看出，实验出现了明显的分岔与混沌的现象，通过图像以及数据结果，我们可以大体上判断分岔点的位置。改变b值继续观察。可以发现，b的更改对实验结果没有本质的影响，从坐标上看我们发现迭代值发生了变化，但是它们所构成的图像相差无几。3、计算分岔点由xk+1=a*xk*exp(-b*xk)得xk+2=a2*xk*exp(-b*xk*(1+a*exp(-b*xk)令xk+2=xk得到1=a2*exp(-b*xk*(1+a*exp(-b*xk)由稳定性有y'=xk+2'=1联立两个式子，进行化简，令t=b*x最终得到exp(-x(1)+1/(x(2)*(1-x(1)= 0exp

20、(x(1)/(1-x(1)-(x(1)-1)/x(2)= 0其中x(1)= t = b*x, x(2)= a。使用fsolve方法求解，得到：x*b = 1.999999993196870a1=7.3884再由y = -1，同理得到x*b=0.74621a2=12.570对于更高阶的分岔点，由于计算非常复杂，这里没有实际计算，仅从图中大致读出其值。a3=14.30a4=14.624、观察是否符合规律A=an-an-1an+1-ann=2时，A= 2.8593n=3时，A= 5.5969有向4.6692收敛的态势。5、小结 本道题应当说是这次作业中最复杂的一道，题目意思明确，但是实验过程中则需测试相当多的数据。另一方面，在计算分岔点时，由于题目所给迭代函数为超越函数，因此计算很复杂，但是这一项本身就又需要用到刚刚学习的非线性方程求解。当然，这也很好地锻炼了我们解决实际问题的能力。总体来说收获较大。程序清单1、迭代函数function y=iter(x,a,b)y=a*x*exp(-b*x);2、b为变量的迭代function chaos1(iter_fun,a,b,n)x0=1;kr=0;for b