数学实验报告四

数学实验报告四非线性方程的数值解2013012472生医三汪一凡2015/4/5一实验目的掌握用Matlab软件求解非线性方程额方程组的基本用法，并对结果作初步分析。练习用非线性方程和方程组建立实际问题的模型并进行求解。二．实验内容：1.1小张夫妇以按揭方式贷款买了一套价值20万的房子，首付5万元，每月还款1000元，15年还清。问贷款利率是多少。1.2某人欲贷款50万元购房，他咨询了两家银行。第一家开出的条件是每月还款4500元。第二家银行开出的条件是每年还45000元，20年还清。从利率方面哪一家更优惠？（简单假设年利率等于月利率*12）问题分析：假定在购买发生后的第n个月（有时是年，不过在本题中可以不用考虑年利率到月利率的转换问题）贷款者还清贷款。若购买时商品价值为P0，月利率为r，首付x0，每月还款x则在n个月之后，商品价值为P0（1+r）nx0（1+r）n+x+x2+……+x由还款基本知识可知，x0（1+r）n+x*1-（1+r）nr=P01+rn（1为了解答本题，首先写出函数monthrate：functiony=monthrate(p0,x0,x,n,r)y=(p0-x0)*(1+r)^n+x*((1+r)^n-1)/r;end下面运用monthrate函数求解小张夫妇的月利率：>>p0=200000;>>x0=50000;>>x=1000;>>n=180;>>r0=0.001;>>[r,fv,ef,out]=fzero(@monthrate,r0,[],p0,x0,x,n);但是这样总是会显示Errorusingfzero(line309)Functionvalueatstartingguessmustbefiniteandreal.具体原因是什么我也不清楚，减小倍数后给出的答案竟然不相同。用inline语句检查方程没有问题。难道是程序自己的故障？为此我还是采取运用inline语句。>>f=inline('150000*r*(1+r)^180/((1+r)^180-1)-1000');>>r=fzero(f,0.001)解得r=0.00208116388946即约为0.21%。对于1.2，同样可以采用inline语句>>f1=inline('500000*r*(1+r)^180/((1+r)^180-1)-4500');>>r1=fzero(f1,0.001)r1=0.005850792582845第二家银行>>f2=inline('500000*r*(1+r)^20/((1+r)^20-1)-45000');>>r2=fzero(f2,0.01);>>r2=0.063948777092386月利率=r2/12=0.51%两者相比显然第二家月利率较小，所以选择第二家银行。小结：需要注意的是，对于这样的银行贷款问题，实际生活中，年利率和月利率绝对不能简单的乘以12倍了事。但是一般以年为单位的利率为相对低一些。同时虽然第二家利率低，但因为年限较长，考虑到CPI等一些通胀因素，不一定会比第一种方案更为优惠。如果物价变化比较大，银行也一般不会改变一份合同中的利率，但不同的经济情况有不同的方案，这些都是要综合比较的。在一些情况需要考虑复利问题，那么式子就会变得稍微更加复杂。如果1.2考虑首付，由计算方法看出也不会出现方案一更优惠的结论。由汽缸控制关闭的门，关闭状态的示意图如图。门宽a，门枢在H处，与H相距b出有一门销，通过活塞与圆柱形的汽缸相连，活塞半径r，汽缸长l0汽缸内气体的压强p0。当用力F推门，使门打开一个角度α时（示意图），活塞下降的距离为c，门销与H的水平距离b保持不变，于是汽缸内的气体被压缩，对活塞的压强增加。已知在绝热条件下，气体的压强p和体积V满足pVγ=C，其中γ是绝热系数，C是常数。试利用开门力矩和作用在活塞上的力矩相平衡的关系（对门枢而言），求在一定的力F作用下，门打开的角度α。设a=0.8m，b=0.25m，r=0.04m，l0=0.5m，如图，在门打开一定的角度α时，由题目各个条件可得btanα=c，pVγFa由此方程写出函数alpfunctiony=alp(x)a=0.8;b=0.25;r=0.04;l0=0.5;p0=10^4;gamma=1.4;F=25;y=(l0/(l0-b*tan(x)))^gamma*p0*pi*r^2*b-F*a*cos(x)end>>x=fzero（@alp，0.01）>>x2=x*180/pi;>>x2解之得x2=24.81°，即门打开了24.81°。小结：这样数学物理结合的问题最主要要找到其对应的数学模型与其原理。力矩这一块我不是很熟悉了已经，还要重新找资料。对于物理量之间的变化关系，“牵一发而动全身”也要特别注意。3用迭代公式xk+1=axkexp⁡(-bxk)计算序列xk分析其收敛性，其中a分别取5，11，15；分析：非线性方程在求解过程中，可能出现迭代中得到的解相差甚远的情况，这是由于其自身属性造成的。一旦参数不在收敛域的范围之中，得到的解释一片混沌。在这里xk+1=axkexp⁡(-bxk)的收敛情况由方程x=axexp⁡(-bx)决定。X1=0，x2=lna/b,若x2<1,则有，平衡点收敛。所以首先按照书上的程序编写chaos.mfunctionchaos(iter_fun,x0,r,n)kr=0;forrr=r(1):r(3):r(2)kr=kr+1y(kr,1)=feval(iter_fun,x0,rr);fori=2:n(2)y(kr,i)=fetal(iter_fun,y(kr,i-1),rr);endendplot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');然后编写函数iter.mfunctiony=iter(x,a)y=a*x*exp(-5*x);end之后写主程序（以a为5为例）chaos(@iter,5,[1,30,0.02],[100,200]);a=5;n=40;x(1)=1;fork=1:(n-1)x(k+1)=a*x(k)*exp(-3*x(k));endk=1:n;plot(k,x);由程序得到以下图：由此图看出，a=5，b=5时，函数能够渐渐收敛趋于稳定，收敛于0.3218867。a=11时，迭代值成周期性变化，收敛到上下两个函数值；可以计算出它们分别收敛到0.794826和0.1643317。此时可以发现函数已经没有明显的收敛子列，不再收敛。所以混沌现象在a=15时已经开始出现下面来求分叉点：还是可以回到这张图，利用光标功能我们会发现，在a在1到7.26时，只有一个极限，没有分叉点，在x在7.26到12.42时有两个收敛极限，开始有分叉点。在12.42到14.22左右有4个极限，有更多的分叉点，在14.22到14.66时有8个收敛极限。此后开始一片混沌。在22.26开始回归2个极限，到23.52处有4个极限，24.24处有8个极限，然后又开始一片混沌。观察7.26,12.42,14.22,14.66和22.26，23.52，24.24两组数：12.42-7.2614.22-12.4214.22-12.4214.66-14.2223.52-22．应该说这和Feigenbaum常数还是有接近的地方，但可能因为读数不精确，取值没有足够大，所以再加上本身数字精度要求很高，所以误差比较大，但是在做到足够多的数据之后应该和Feigenbaum常数的4.6692会越来越接近。小结：这道题主要关注的是非线性方程中的混沌现象，具体从“感觉上”其实很难想象到“蝴蝶效应”的影响之大，这更加体验出运用matlab分析的重要意义。随着a的不断增大，函数图像的分叉也越来越多。但是当到某一个点之后可能又会一下子只有两个分叉，但是无论如何，整体的混沌趋向性是一致的。在运用程序计算时，因为舍入误差，其实显示的混沌图像也许不能精确的反映实际情况，这也是我们要牢记的一点。实验