数学实验非线性方程求解.docx
- 文档编号:8885946
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:19
- 大小:872.42KB
数学实验非线性方程求解.docx
《数学实验非线性方程求解.docx》由会员分享,可在线阅读,更多相关《数学实验非线性方程求解.docx(19页珍藏版)》请在冰豆网上搜索。
数学实验非线性方程求解
实验6非线性方程求解
分1黄浩43
一、实验目的
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)所提供的“负本金”:
则还清的标志为:
本金额=负本金的总额
公式即:
整理后得:
即:
为了寻找比较合适的初始值,我们先编写上述方程的matlab函数(程序见四.2),然后对
进行作图(程序见四.1),得:
上图中,横坐标为r,纵坐标为方程(*)左侧的值。
由上图可见,当i在区间[0.002,0.0025]内,方程有零点。
因此,令A=150000,N=180,R=1000,有根区间为[0.002,0.0025],使用fzero函数(程序见四.3)解得:
r=0.0089459
fv=4.0978177e-012
可见,fv已经十分接近于0,结合之前的图像,可以判定是一个零点而非近似间断点。
因此,贷款利率r≈0.21%。
得出结论:
贷款利率为0.21%
(2)对于这两种情况,仍然使用上一问的基本模型,代入两种情况下的参数值,构建两个方程:
其中r1为第一家银行的月利率,r2为第二家银行的年利率。
为了寻找比较合适的初始值,对
进行作图(程序见四.4):
由上述两图可见,当
时,方程有零点。
分别对上述两种情况,使用fzero函数求解(程序见四.5),得:
r1=0.0082846≈0.59%fv1=1.2482100e-010
r2=0.0092386≈6.39%fv2=-6.2765083e-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。
模型转换及实验过程:
根据力矩平衡,可得:
因为外压没有说明,而设外压为1个大气压显然不合理(因为
),所以可以简化模型,忽略外压的影响。
即:
又根据绝热过程的公式:
得:
代入第一式,得:
为了寻找比较合适的初始值,先编写上述方程的matlab函数(程序见四.6),然后对函数
在
上进行作图(程序见四.7),结果如下:
由上图可见,当
时,方程有零点。
然后,使用matlab解该方程的零点(程序见四.8),得:
由fv=0,及上面的图像,可以判定所得解为零点而非近似间断点。
结论:
在力F=25N的作用下,门打开的角度是24.81°
关于此题的小讨论:
个人感觉这道题有很多漏洞,在审题的时候,存在较大的障碍,最后的答案实际上是不断修改模型,不断与书后答案对照而“凑”出来的,现将其中的漏洞一一列出:
第一,原图的标注有问题。
此题的被控物是一道门,而根据经验,门的宽度一般是固定的,不会发生伸缩,但在原图(第一版,P135,图6.12)中,图(b)门的长度竟成了
,门的水平投影没有变化。
而我认为,真实的情况应该是:
门的长度不变,而门的水平投影变为
,这样的改动也与答案相符。
第二,关于汽缸外的压强。
如果要推出答案的结果,必须假设外压为0,而这与实际情况是完全不相符合的。
我认为,此题的问题在于将p0设的过小,才0.1atm,这根本不可能使门在没有外力的条件下保持水平,即由于汽缸活塞内压远小于正常大气压,门会自动“挤”开,直至内压=1atm。
我使用计算器计算了此时的角度,为58.22°,已经远远大于本题的结果:
24.81°。
显然,如果考虑外压,则真正角度会远远大于书后答案。
第三,此题的绝热气体公式
=C的适用条件为绝热可逆过程,即我们需要无限缓慢地推门,使气体始终处于力学平衡和热平衡,才能完全适用于这个模型。
而实际情况显然不是这样的,我们往往在一瞬间就完成了推门过程,可逆的要求是不能满足的。
这也是本题的一点漏洞。
3.《数学实验》第一版(问题7)
问题描述:
用迭代公式
计算序列
,分析其收敛性,其中a分别取5,11,15;b(>0)任意,初值x0=1。
观察是否有混沌现象出现,并找出前几个分岔点,观察分岔点的极限趋势是否符合Feigenbaum常数揭示的规律。
模型转换及实验过程:
直接使用matlab编程,固定b=1,分别调整a=5,11,15,迭代值
与k之间的图像如下:
当a=5时(程序见四.9):
由上图可见,最终序列
收敛到唯一值,即有唯一的收敛子列。
当a=11时(程序略):
a=11时,似乎有两个个收敛子列,但从该图像上无法观察仔细,我们将奇数子列与偶数子列分开制图(程序见四.10),得:
此图显然明了了许多,a=11时的确有两个收敛子列。
而当a=15时(程序见四.11),迭代了100次,但很难找到收敛的痕迹:
初步判定,当a=15时出现混沌状态。
因为处于混沌状态的系统,对于初值十分敏感,因此,将初值改为1.001,与之前的图像进行对比(程序见四.12,取迭代次数为100~200的数据作图,蓝、绿线分别为更改初值前后的迭代值):
如图,当a=5、11时,对初值的变化并不敏感,当a=15时,初值会对迭代产生明显影响。
a=15时,新旧迭代值的差值(程序见四.13)如下图所示:
由上图可见,在混沌状态下,初值不同而造成的迭代值改变关于迭代次数成近似的周期性关系,且差值时大时小,最大时有0.8的绝对误差。
因初值的相对误差为:
而迭代值的最大相对误差(图中取点而得)近似为:
,即1个单位的初值改变会造成最大143倍的迭代值改变!
为了观察当参数a变化时,迭代值的分岔与混沌情况,通过编制chaos函数(程序见四.14)与迭代函数(程序见四.15),取a=0:
0.01:
25,利用迭代次数为200到500次的数据为锚点,输出了收敛、分岔和混沌状态图(程序见四.16.第1行):
由上图也可以佐证a=15时出现混沌态的判断。
为了观察分岔点,我们将上图放大后(放大图略),标记分岔点,得到:
分别计算
可得到三组数据:
n=1时,
n=2时,
n=3时,
因此,可以认为
,即分叉点的极限趋势符合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)贷款模型的函数
functiony=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.第二题——汽缸模型的函数
functiony=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;
fori=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;
fori=2:
50
x(i)=11*x(i-1)*exp(-x(i-1));
end
k=1:
50;
k1=1:
2:
50;
k2=2:
2:
50;
fori=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;
fori=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];
forj=1:
3
fori=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;
forj=1:
3
subplot(3,1,j);
plot(k(100:
200),x(j,100:
200),k(100:
200),xx(j,100:
200));
grid;axis([10020006])
end
13.第三题——混沌状态下,初值改变对于迭代值的改变(差值)
x=1;
xx=1.001;
a=15;
fori=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
functionchaos(iterfun,x0,r,n)
kr=0;
forrr=r
(1):
r(3):
r
(2)
kr=kr+1;
y(kr,1)=feval(iterfun,x0,rr);
fori=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.第三题——迭代函数
functiony=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;
fori=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;
fori=2:
50
x(i)=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;
fori=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;
fori=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
仅供个人用于学习、研究;不得用于商业用途。
Forpersonaluseonlyinstudyandresearch;notforcommercialuse.
NurfürdenpersönlichenfürStudien,Forschung,zukommerziellenZweckenverwendetwerden.
Pourl'étudeetlarechercheuniquementàdesfinspersonnelles;pasàdesfinscommerciales.
толькодлялюдей,которыеиспользуютсядляобучения,исследованийинедолжныиспользоватьсявкоммерческихцелях.
以下无正文
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数学 实验 非线性 方程 求解