数学实验——非线性方程求解.docx

数学实验——非线性方程求解

———————————————————————————————— 作者：
———————————————————————————————— 日期：

2

数学实验报告
实验6 非线性方程求解

实验6 非线性方程求解
分1 黄浩 2011011743
一、 实验目的
1. 掌握用MATLAB软件求解非线性方程和方程组的基本用法，并对结果作初步分析。
2. 练习用非线性方程和方程组建立实际问题的模型并进行求解。
二、 实验内容
1.《数学实验》第一版（问题3）
问题叙述：
（1）小张夫妇以按揭方式贷款买了一套价值20万的房子，首付了5万元，每月还款1000元，15年还清。问贷款利率是多少？
（2）某人欲贷款50万元购房，他咨询了两家银行，第一家银行开出的条件是每月还4500元，15年还清；第二家银行开出的条件是每年还45000元，20年还清。从利率方面看，哪家银行较优惠（简单的假设年利率=月利率*12）？
模型转换及实验过程：
（1）本题的按揭贷款属于等额还款类型，即每月的还款额度相同，还款额先抵消当前本金的月利息，然后剩余还款额用来偿还本金。再进一步简化等效，可以设想本金额不被抵消，始终按照一定的利率呈指数增长，同时所交款成为“负本金”，也按照利率成指数增长，二者的差值即为待交的款额。
设本金为A=15万，从第一个月(n=1）开始，贷款月利率为i，每月交款额为R=1000元。
则第n月所交款在还清时(n=N)所提供的“负本金”：
an=R(1+r)N-n
则还清的标志为：本金额=负本金的总额
公式即：
A(1+r)N=1Naj=R(1+r)N-1r
整理后得：
R=Ar1+rN1+rN-1
即：
Ar-R1+rN+R=0…………(*)
为了寻找比较合适的初始值，我们先编写上述方程的matlab函数（程序见四.2），然后对r∈[0,0.005]进行作图（程序见四.1），得：

上图中，横坐标为r，纵坐标为方程(*)左侧的值。由上图可见，当i在区间[0.002,0.0025]内，方程有零点。
因此，令A=150000，N=180，R=1000，有根区间为[0.002,0.0025]，使用fzero函数（程序见四.3）解得：
r = 0.002081163889459
fv = 4.092726157978177e-012
可见，fv已经十分接近于0，结合之前的图像，可以判定是一个零点而非近似间断点。因此，贷款利率r≈0.21%。
得出结论：贷款利率为0.21%
（2）对于这两种情况，仍然使用上一问的基本模型，代入两种情况下的参数值，构建两个方程：
(500000*r1-4500)*(1+r1)^180+4500=0 (500000*r2-45000)*(1+r2)^20+45000=0
其中r1为第一家银行的月利率，r2为第二家银行的年利率。
为了寻找比较合适的初始值，对r1∈0,0.01,r2∈0.01,0.1进行作图（程序见四.4）:

由上述两图可见，当r1∈0005,0.006,r2∈0.06,0.07时，方程有零点。
分别对上述两种情况，使用fzero函数求解（程序见四.5），得：
r1 = 0.005850792582846≈0.59% fv1 =1.273292582482100e-010
r2 = 0.063948777092386≈6.39% fv2 =-6.548361852765083e-011
因为fv1、fv2都很接近0，结合上面的两幅图像，可以判定r1、r2分别为两方程的零点而非近似间断点
因此，第一家银行的年利率R1 = 0.59%*12 = 7.08% > R2 = r2 = 6.39%
得出结论：
第二家银行更优惠。
2.《数学实验》第一版（问题5）
问题叙述：
由汽缸控制关闭的门，门宽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=10000N/m2，γ=1.4，F=25N。
模型转换及实验过程：
根据力矩平衡，可得：
Facosα=p-p外S0∙cosα∙bcosα
因为外压没有说明，而设外压为1个大气压显然不合理（因为p0=104pa≪1atm=105pa）,所以可以简化模型，忽略外压的影响。即：
Facosα=pS0∙b
又根据绝热过程的公式：
p0V0γ=pVγ
得：
p=(l0l0-c)γp0
代入第一式，得：
Facosα=(l0l0-btanα)γp0πr2b
为了寻找比较合适的初始值，先编写上述方程的matlab函数（程序见四.6），然后对函数y=Facosα-(l0l0-btanα)γp0πr2b在α∈[0,π4]上进行作图（程序见四.7），结果如下：

由上图可见，当α∈[0.4,0.5]时，方程有零点。
然后，使用matlab解该方程的零点（程序见四.8），得：
α= 0.432945646735310≈24.81° fv=0
由fv=0，及上面的图像，可以判定所得解为零点而非近似间断点。
结论：
在力F=25N的作用下，门打开的角度是24.81°
关于此题的小讨论：
个人感觉这道题有很多漏洞，在审题的时候，存在较大的障碍，最后的答案实际上是不断修改模型，不断与书后答案对照而“凑”出来的，现将其中的漏洞一一列出：
第一，原图的标注有问题。此题的被控物是一道门，而根据经验，门的宽度一般是固定的，不会发生伸缩，但在原图（第一版，P135，图6.12）中，图(b)门的长度竟成了acosα，门的水平投影没有变化。而我认为，真实的情况应该是：门的长度不变，而门的水平投影变为acosα，这样的改动也与答案相符。
第二，关于汽缸外的压强。如果要推出答案的结果，必须假设外压为0，而这与实际情况是完全不相符合的。我认为，此题的问题在于将p0设的过小，才0.1atm，这根本不可能使门在没有外力的条件下保持水平，即由于汽缸活塞内压远小于正常大气压，门会自动“挤”开，直至内压=1atm。我使用计算器计算了此时的角度，为58.22°,已经远远大于本题的结果：24.81°。显然，如果考虑外压，则真正角度会远远大于书后答案。
第三，此题的绝热气体公式pVγ=C的适用条件为绝热可逆过程，即我们需要无限缓慢地推门，使气体始终处于力学平衡和热平衡，才能完全适用于这个模型。而实际情况显然不是这样的，我们往往在一瞬间就完成了推门过程，可逆的要求是不能满足的。这也是本题的一点漏洞。
3.《数学实验》第一版（问题7）
问题描述：
用迭代公式xk+1=axkexp⁡(-bxk)计算序列{xk}，分析其收敛性，其中a分别取5,11,15；b(>0)任意，初值x0=1。观察是否有混沌现象出现，并找出前几个分岔点，观察分岔点的极限趋势是否符合Feigenbaum常数揭示的规律。
模型转换及实验过程：
直接使用matlab编程，固定b=1，分别调整a=5,11,15，迭代值xk与k之间的图像如下：
当a=5时（程序见四.9）：

由上图可见，最终序列{xk}收敛到唯一值，即有唯一的收敛子列。
当a=11时（程序略）：

a=11时，似乎有两个个收敛子列，但从该图像上无法观察仔细，我们将奇数子列与偶数子列分开制图（程序见四.10），得：

此图显然明了了许多，a=11时的确有两个收敛子列。
而当a=15时（程序见四.11），迭代了100次，但很难找到收敛的痕迹：

初步判定，当a=15时出现混沌状态。因为处于混沌状态的系统，对于初值十分敏感，因此，将初值改为1.001，与之前的图像进行对比（程序见四.12，取迭代次数为100~200的数据作图，蓝、绿线分别为更改初值前后的迭代值）：
a=5
a=11
a=15

如图，当a=5、11时，对初值的变化并不敏感，当a=15时，初值会对迭代产生明显影响。a=15时，新旧迭代值的差值（程序见四.13）如下图所示：

由上图可见，在混沌状态下，初值不同而造成的迭代值改变关于迭代次数成近似的周期性关系，且差值时大时小，最大时有0.8的绝对误差。因初值的相对误差为：1.001-11=0.001,而迭代值的最大相对误差（图中取点而得）近似为：0.85.6=0.143=143*0.001，即1个单位的初值改变会造成最大143倍的迭代值改变！
为了观察当参数a变化时，迭代值的分岔与混沌情况，通过编制chaos函数（程序见四.14）与迭代函数（程序见四.15），取a=0:0.01:25，利用迭代次数为200到500次的数据为锚点，输出了收敛、分岔和混沌状态图（程序见四.16.第1行）：

由上图也可以佐证a=15时出现混沌态的判断。
为了观察分岔点，我们将上图放大后（放大图略），标记分岔点，得到：
a0=7.32 a1=12.45 a2=14.23 a3=14.65 a4=14.74
分别计算an-an-1an+1-an,可得到三组数据：
n=1时，an-an-1an+1-an=2.882
n=2时，an-an-1an+1-an=4.238
n=3时，an-an-1an+1-an=4.667
因此，可以认为limn→∞an-an-1an+1-an=4.6692，即分叉点的极限趋势符合Feigenbaum常数解释的规律。
此外，通过观察上图，我们还不难注意到，在a=16和18.5时，似乎出现了一条缝隙，而在a=22~25的区间内，仿佛混沌现象已经消失了，迭代重新归于有序的收敛。将a的取值扩大为40，如下图所示（程序见四.16.第2行）：

由上图可见，当a=36~37的区间内，同样出现了非混沌区。为了观察、验证这些非混沌区，我取a=16、18.5、23、37，将他们的迭代值绘制成曲线（程序见四.17）：

由上图可见，当a=16、18.5、23、37时，均出现了有序的收敛现象，而且，通过观察放大图，可以辨明，这四幅图分别是周期6、周期5、周期3、周期8收敛，这说明在一段混沌之后，系统又重新归于有序的周期性的收敛，而且随着收敛子列的增加，又会再次出现混沌的情况。
将a的取值扩大到100，得混沌图如下（程序见四.16.第3行）：

由上图可以观察到，当a增加时，迭代序列的混沌态与收敛态会交替出现，且随a的增加，收敛态的宽度越来越小，收敛态之间的间隔越来越大，也即混沌态的宽度越来越大，混沌态之间的间隔也越来越小。
得出结论：
1) 当a=5时，有唯一的收敛序列，当a=11时，有两个收敛子列，当a=15时，出现了混沌现象，混沌态的系统对于初值变化十分敏感。
2) 观察到的分叉点的极限趋势符合Feigenbaum常数解释的规律。
3) 随着参数a的增加，混沌态与收敛态交替出现，且随a增加，收敛态的宽度越来越小，收敛态之间的间隔越来越大。
三、 实验总结
本次实验是求非线性方程的数值解，通过求解两道实际应用题与一道理论数学题目，初步达到了实验要求的目标，掌握了使用matlab求解非线性方程的基本方法。此外，在第二题中，针对题目的漏洞，提出了自己的看法，而在最后一题中，针对实验过程中出现的一些新奇的现象，还进行了自主探究，并得到了相应的结论。
当然，本次实验也经历了很多困难，主要是前两道题，难点在于建立模型，即构建贷款的公式（第一题）和分析力矩平衡（第二题），由于缺乏相应的理论知识，我在这两块上花费了很多时间。这也说明，只有专业素养，而缺乏全面丰富的社科知识，是很难适应当今学科交叉的大环境的。这也是我今后要努力的方向。
四、 程序清单
1. 第一题——（1）小题作图寻找有根区间
fplot('(150000*x-1000)*(1+x)^180+1000',[0,0.005])
grid

2. 第一题——（1）贷款模型的函数
function y=monthmoney(r,A,N,R)
y=(A*r-R)*(1+r)^N+R;
end

3. 第一题——（1）求解方程
A=150000;
N=180;
R=1000;
[r,fv]=fzero(@monthmoney,[0.002,0.0025],[],A,N,R)

4. 第一题——（2）作图寻找有根区间
fplot('(500000*r-4500)*(1+r)^180+4500',[0.001,0.01])
grid
xlabel('月利率');
pause
fplot('(500000*r-45000)*(1+r)^20+45000'[0.01,0.1])
grid
xlabel('年利率');

5. 第一题——（2）求解方程
A=500000;
N=180;
R=4500;
[r1,fv1]=fzero(@monthmoney,[0.005,0.006],[],A,N,R)
A=500000;
N=20;
R=45000;
[r2,fv2]=fzero(@monthmoney,[0.06,0.07],[],A,N,R)

6. 第二题——汽缸模型的函数
function y=angle(al)
k=2/(2-tan(al));
y=20*cos(al)-4*pi*k^1.4;
end

7. 第二题——作图寻找有根区间
fplot('20*cos(al)-4*pi*(2/(2-tan(al)))^1.4',[0,pi/4])
grid
xlabel('角度')

8. 第二题——求解方程
[al,fv]=fzero(@angle,[0.4,0.5])

9. 第三题——a=5的情况
x(1)=1;
for i=2:50
x(i)=5*x(i-1)*exp(-x(i-1));
end
k=1:50;
plot(k,x);
xlabel('迭代次数');
ylabel('迭代值');
title('a=11');
grid

10. 第三题——a=11的情况
x(1)=1;
for i=2:50
x(i)=11*x(i-1)*exp(-x(i-1));
end
k=1:50;
k1=1:2:50;
k2=2:2:50;
for i=1:25
x1(i)=x(2*i-1);
x2(i)=x(2*i);
end
subplot(3,1,1);
plot(k,x);
grid;ylabel('迭代值');title('迭代子列全集');
subplot(3,1,2);
plot(k1,x1);
grid;ylabel('迭代值');title('奇数子列');
subplot(3,1,3);
plot(k2,x2);
grid;ylabel('迭代值');title('偶数子列');

11. 第三题——a=15的情况
x(1)=1;
for i=2:100
x(i)=15*x(i-1)*exp(-x(i-1));
end
k=1:100;
plot(k,x);
grid

12. 第三题——混沌状态对于初值的敏感特性
x=[1;1;1];
xx=[1.001;1.001;1.001];
a=[5,11,15];
for j=1:3
for i=2:200
x(j,i)=a(j)*x(j,i-1)*exp(-x(j,i-1));
xx(j,i)=a(j)*xx(j,i-1)*exp(-xx(j,i-1));
end
end
k=1:200;
for j=1:3
subplot(3,1,j);
plot(k(100:200),x(j,100:200),k(100:200),xx(j,100:200));
grid;axis([100 200 0 6])
end

13. 第三题——混沌状态下，初值改变对于迭代值的改变（差值）
x=1;
xx=1.001;
a=15;
for i=2:200

x(i)=a*x(i-1)*exp(-x(i-1));
xx(i)=a*xx(i-1)*exp(-xx(i-1));
end
k=1:200;
plot(k(100:200),x(100:200)-xx(100:200));
grid;

14. 第三题——迭代序列随参数变化的分岔和混沌图函数chaos
function chaos(iterfun,x0,r,n)
kr=0;
for rr=r(1):r(3):r(2)
kr=kr+1;
y(kr,1)=feval(iterfun,x0,rr);
for i=2:n(2)
y(kr,i)=feval(iterfun,y(kr,i-1),rr);
end
end
plot([r(1):r(3):r(2)],y(:,n(1)+1:n(2)),'k.');
title('收敛、分岔和混沌现象图');
xlabel('参数a的取值');
ylabel('迭代值');
end

15. 第三题——迭代函数
function y=iter01(x,a)
y=a*x*exp(-x);
end

16. 第三题——输出分岔与混沌图（a取值从0分别到25、40、100）
chaos(@iter01,1,[0,25,0.01],[200,500])
chaos(@iter01,1,[0,40,0.01],[200,500])
chaos(@iter01,1,[0,100,0.01],[200,800])

17. 第三题——输出a=16、18.5、23、37时的迭代值曲线
x(1)=1;
for i=2:50
x(i)=16*x(i-1)*exp(-x(i-1));
end
k=1:50;
subplot(2,2,1);
plot(k,x);xlabel('迭代次数');ylabel('迭代值');title('a=16');grid
x(1)=1;
for i=2:50
xi)=18.5*x(i-1)*exp(-x(i-1));
end
k=1:50;
subplot(2,2,2);
plot(k,x);xlabel('迭代次数');ylabel('迭代值');title('a=18.5');grid
x(1)=1;
for i=2:50
x(i)=23*x(i-1)*exp(-x(i-1));
end
k=1:50;
subplot(2,2,3);
plot(k,x);xlabel('迭代次数');ylabel('迭代值');title('a=23');grid
x(1)=1;
for i=2:50
x(i)=37*x(i-1)*exp(-x(i-1));
end
k=1:50;
subplot(2,2,4);
plot(k,x);xlabel('迭代次数');ylabel('迭代值');title('a=37');grid

20