数学实验报告四.doc

数学实验报告四

非线性方程的数值解

2013012472
生医三 汪一凡
2015/4/5

一 实验目的
1. 掌握用Matlab软件求解非线性方程额方程组的基本用法，并对结果作初步分析。
2. 练习用非线性方程和方程组建立实际问题的模型并进行求解。
二．实验内容：
1.1小张夫妇以按揭方式贷款买了一套价值20万的房子，首付5万元，每月还款1000元，15年还清。问贷款利率是多少。
1.2某人欲贷款50万元购房，他咨询了两家银行。第一家开出的条件是每月还款4500元。第二家银行开出的条件是每年还45000元，20年还清。从利率方面哪一家更优惠？（简单假设年利率等于月利率*12）
问题分析：
假定在购买发生后的第n个月（有时是年，不过在本题中可以不用考虑年利率到月利率的转换问题）贷款者还清贷款。若购买时商品价值为P0，月利率为r，首付x0，每月还款x则在n个月之后，商品价值为P0（1+r）n，在这n个月中，他所支付的钱的总价值为
x0（1+r）n+x+x2+……+xn=x0（1+r）n+x*1-（1+r）nr
由还款基本知识可知，x0（1+r）n+x*1-（1+r）nr= P01+rn（1）。在小张夫妇的那一问中，我们已知x0=5，P0=20，x=0.1，n=180，即可由（1）式求出贷款月利率。
为了解答本题，首先写出函数monthrate：
function y=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);
但是这样总是会显示Error using fzero (line 309)
Function value at starting guess must be finite and real.具体原因是什么我也不清楚，减小倍数后给出的答案竟然不相同。用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考虑首付，由计算方法看出也不会出现方案一更优惠的结论。

3. 由汽缸控制关闭的门，关闭状态的示意图如图。门宽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，p0 =104Nm2，γ=1.4，F=25N。

如图，在门打开一定的角度α时，由题目各个条件可得btanα=c，pVγ=C，S=πr2，V0=Sl0 ，V=S(l0 –c)，p=（l0 l0 –C）γp0 。假设推力始终与门所在平面垂直，假设此时气缸的气压为p，则根据受力力矩相等，可得：Facosα=pSb，结合以上各式，得到
Facosα=（l0 l0 –btanα）γp0 πr2b
由此方程写出函数alp
function y=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；b(>0)任意，初值x0=1，观察是否有混沌现象出现，并找出前几个分岔点，观察分岔点的极限趋势是否符合Feigenbaum常数提示的规律。
分析：
非线性方程在求解过程中，可能出现迭代中得到的解相差甚远的情况，这是由于其自身属性造成的。一旦参数不在收敛域的范围之中，得到的解释一片混沌。在这里xk+1=axkexp⁡(-bxk)的收敛情况由方程x=axexp⁡(-bx)决定。X1=0，x2=lna/b,若x2<1,则有，平衡点收敛。所以b>lna.因为a取值为5,11,15，可以取b=5.
首先按照书上的程序编写chaos.m
function chaos(iter_fun,x0,r,n)
kr=0;
for rr=r(1):r(3):r(2)
kr=kr+1
y(kr,1)=feval(iter_fun,x0,rr);
for i=2:n(2)
y(kr,i)=fetal(iter_fun,y(kr,i-1),rr);
end
end
plot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');
然后编写函数iter.m
function y=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;
for k=1:(n-1)
x(k+1)=a*x(k)*exp(-3*x(k));
end
k=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.42=2.867
14.22-12.4214.66-14.22=4.091
23.52-22．2624.24-23.52=1.78
应该说这和Feigenbaum常数还是有接近的地方，但可能因为读数不精确，取值没有足够大，所以再加上本身数字精度要求很高，所以误差比较大，但是在做到足够多的数据之后应该和Feigenbaum常数的4.6692会越来越接近。
小结：
这道题主要关注的是非线性方程中的混沌现象，具体从“感觉上”其实很难想象到“蝴蝶效应”的影响之大，这更加体验出运用matlab分析的重要意义。随着a的不断增大，函数图像的分叉也越来越多。但是当到某一个点之后可能又会一下子只有两个分叉，但是无论如何，整体的混沌趋向性是一致的。在运用程序计算时，因为舍入误差，其实显示的混沌图像也许不能精确的反映实际情况，这也是我们要牢记的一点。

实验总结：
本次试验因为书上程序比较完整，所以调试的时间相对减少了很多。通过对一系列实际生活中出现的非线性方程的求解，我相对清楚的掌握了matlab fzero函数的使用以及一些常见模型的建立。这些来自于生活的问题让我更加深刻的感受到了数学模型和思想的重要性。但是坦诚说对于混沌效应，其发展原因趋势等等，有很多地方我还是很不清楚的，希望今后能够有进一步的深入了解。