机械 2592 夏俊超 实验二.docx
《机械 2592 夏俊超 实验二.docx》由会员分享，可在线阅读，更多相关《机械 2592 夏俊超 实验二.docx（30页珍藏版）》请在冰点文库上搜索。
机械2592夏俊超实验二
重庆大学
学生实验报告
实验课程名称数学实验
开课实验室DS1401
学院机械工程年级2009级专业班机自14班
学生姓名夏俊超学号20092592
开课时间2010至2011学年第一学期
总成绩
教师签名
数理学院制
开课学院、实验室：
数学学院DS1401实验时间：
2011年3月30日
课程
名称
数学实验
实验项目
名称
实验项目类型
验证
演示
综合
设计
其他
指导
教师
成绩
实验目的
[1]复习求解方程及方程组的基本原理和方法；
[2]掌握迭代算法；
[3]熟悉MATLAB软件编程环境；掌握MATLAB编程语句（特别是循环、条件、控制等语句）；
[4]通过范例展现求解实际问题的初步建模过程；
通过该实验的学习，复习和归纳方程求解或方程组求解的各种数值解法（简单迭代法、二分法、牛顿法、割线法等），初步了解数学建模过程。
这对于学生深入理解数学概念，掌握数学的思维方法，熟悉处理大量的工程计算问题的方法具有十分重要的意义。
基础实验
一、实验内容
1．方程求解和方程组的各种数值解法练习
2．直接使用MATLAB命令对方程和方程组进行求解练习
3．针对实际问题，试建立数学模型，并求解。
二、实验过程（一般应包括实验原理或问题分析，算法设计、程序、计算、图表等，实验结果及分析）
第一题：
用图形放大法求解方程xsin（x）=1.
x=1:
0.01:
5;
y=x.*sin（x）;
>>y1=zeros（size（x））;
>>plot（x,y,x,y1）
x=3:
0.001:
3.5;
y=x.*sin（x）;
y1=zeros（size（x））;
plot（x,y,x,y1）
x=3.1:
0.00001:
3.15;
y=x.*sin（x）;
y1=zeros（size（x））;
plot（x,y,x,y1）
x=3.141:
0.00001:
3.143;
y=x.*sin（x）;
y1=zeros（size（x））;
plot（x,y,x,y1）
图中只求了一个根，因为sin（x）是以2π为周期的周期函数所以共有无数个根。
第二题：
．将方程x5+5x3-2x+1=0改写成各种等价的形式进行迭代，观察迭代是否收敛，并给出解释。
解：
第一步构造迭代函数：
x=j（x）
f=（x.^5+5*x.^3+1）/2f（x）
y=（（2*x-1-x.^5）/5）^（1/3）y（x）
z=（2*x-1-5*x^3）^（1/5）z（x）
第二步迭代：
设定初值x0=1
f
（1）=1;y
（1）=1;z
（1）=1;
fork=1:
20
f（k+1）=（f（k）^5+5*（f（k）^3）+1）/2;%f（x）
y（k+1）=（（2*y（k）-1-y（k）^5）/5）^（1/3）;%y（x）
z（k+1）=（2*z（k）-1-5*（z（k）^3））^（1/5）;%z（x）
end
x
y
z
第三步计算结果
x=
Columns1through18
-6.0000-5.9900-5.9800-5.9700-5.9600-5.9500-5.9400-5.9300-5.9200-5.9100-5.9000-5.8900-5.8800-5.8700-5.8600-5.8500-5.8400-5.8300……
y=
1.0e+003*
Columns1through9
0.001000.0003+0.0005i0.0005+0.0004i0.0005+0.0003i0.0004+0.0002i0.0004+0.0003i0.0004+0.0003i0.0004+0.0003i……
z=
Columns1through10
1.00001.0675+0.7756i1.5607-0.3514i1.5818+0.8381i1.8323-0.6099i1.8377+0.8353i1.9443-0.7354i1.9458+0.8285i1.9884-0.7881i1.9889+0.8249i……
第四步结论
f（x）,y（x）,z（x）均收敛！
近似解为x1=6.0000，x2=7.8520,x3=2.0162-0.8222i。
第三题：
用solve（）和fsolve（）对方程组求解
（1）解
[x1,x2]=solve（'2*x1-x2-exp（-x1）','-x1+2*x2-exp（-x2）','x1','x2'）
x1=
0.56714329040978387299996866221036
x2=
0.56714329040978387299996866221036
建立M函数文件（ex.m）
functioneq=ex（x）
eq
（1）=2*x
（1）-x
（2）-exp（-x
（1））;
eq
（2）=-x
（1）+2*x
（2）-exp（-x
（2））;
运行文件（ex.m）
y=fsolve（'ex',[2,2],1）
运行结果：
y=
0.56710.5671
（2）解A
[x1,x2,x3]=solve（'x1^2-5*x2^2+7*x3^2+12','3*x1*x2+x1*x3-11*x1','2*x2*x3+40*x1','x1','x2','x3'）
x1=
0
0
0
0
1
1/20*（1/9*（2052199+60*94819590129^（1/2））^（1/3）-69599/9/（2052199+60*94819590129^（1/2））^（1/3）+7/9）*（1/3*（2052199+60*94819590129^（1/2））^（1/3）-69599/3/（2052199+60*94819590129^（1/2））^（1/3）-26/3）
1/20*（-1/18*（2052199+60*94819590129^（1/2））^（1/3）+69599/18/（2052199+60*94819590129^（1/2））^（1/3）+7/9+1/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3）））*（-1/6*（2052199+60*94819590129^（1/2））^（1/3）+69599/6/（2052199+60*94819590129^（1/2））^（1/3）-26/3+3/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3）））
1/20*（-1/18*（2052199+60*94819590129^（1/2））^（1/3）+69599/18/（2052199+60*94819590129^（1/2））^（1/3）+7/9-1/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3）））*（-1/6*（2052199+60*94819590129^（1/2））^（1/3）+69599/6/（2052199+60*94819590129^（1/2））^（1/3）-26/3-3/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3）））
x2=
-2/5*15^（1/2）
2/5*15^（1/2）
0
0
5
1/9*（2052199+60*94819590129^（1/2））^（1/3）-69599/9/（2052199+60*94819590129^（1/2））^（1/3）+7/9
-1/18*（2052199+60*94819590129^（1/2））^（1/3）+69599/18/（2052199+60*94819590129^（1/2））^（1/3）+7/9+1/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3））
-1/18*（2052199+60*94819590129^（1/2））^（1/3）+69599/18/（2052199+60*94819590129^（1/2））^（1/3）+7/9-1/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3））
x3=
0
0
2/7*i*21^（1/2）
-2/7*i*21^（1/2）
-4
-1/3*（2052199+60*94819590129^（1/2））^（1/3）+69599/3/（2052199+60*94819590129^（1/2））^（1/3）+26/3
1/6*（2052199+60*94819590129^（1/2））^（1/3）-69599/6/（2052199+60*94819590129^（1/2））^（1/3）+26/3-3/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3））
1/6*（2052199+60*94819590129^（1/2））^（1/3）-69599/6/（2052199+60*94819590129^（1/2））^（1/3）+26/3+3/2*i*3^（1/2）*（1/9*（2052199+60*94819590129^（1/2））^（1/3）+69599/9/（2052199+60*94819590129^（1/2））^（1/3））
B．建立方程组的M-函数文件：
functioneq=xf（x）
eq
（1）=（x
（1））^2-5*（x
（2））^2+7*（x（3）^2）+12;
eq
（2）=3*x
（1）*x
（2）+x
（1）*x（3）-11*x
（1）;
eq（3）=2*x（3）*x
（2）+40*x
（1）;
B．运行程序：
y=fsolve（'xf',[1,1,1],1）
C.运行结果：
y=
0.00001.54920.0000
第四题：
迭代以下函数，分析其收敛性。
1）解
（1）线性连接图
%当a等于1时
a=1;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,1）,plot（n,x1）,title（'a=1,x0=0.5'）
%当a等于1.1时
a=1.1;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,2）,plot（n,x1）,title（'a=1.1,x0=0.5'）
%当a等于1.2时
a=1.2;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,3）,plot（n,x1）,title（'a=1.2,x0=0.5'）
%当a等于1.3时
a=1.3;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,4）,plot（n,x1）,title（'a=1.3,x0=0.5'）
%当a等于1.4时
a=1.4;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,5）,plot（n,x1）,title（'a=1.4,x0=0.5'）
%当a等于2时
a=2;
x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a-（x1（i-1）-a^（1/2））^.2;end;
n=1:
20;
subplot（3,3,6）,plot（n,x1）,title（'a=2,x0=0.5'）
（2）蛛网图法
x1=[];a=0.5;
x=-1.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,1）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=0.5'）
end
>>x1=[];a=1;
x=-0.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,2）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=1'）
end
>>x1=[];a=1.5;
x=-0.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,3）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=1.5'）
end
>>x1=[];a=2;
x=-0.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,4）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=2'）
end
>>x1=[];a=2.5;
x=-0.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,5）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=2.5'）
end
>>x1=[];a=3;
x=-0.2:
0.01:
1.2;
y=a.*x.*（1-x）;subplot（3,3,6）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*x1（i-1）*（1-x1（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=3'）
end
3）混沌图法
建立M文件（ft.m）
functionroot=ft（x,a）
x=[];x
（1）=0.2;
fori=2:
100
x（i）=a-（x（i-1）-sqrt（a））^2;
end
root=x;
程序
>>clf
>>x=[];
holdon;
fora=2:
0.01:
4
root=ft（x,a）;
plot（a.*ones（size（root（51:
100）））,root（51:
100）,'.'）
end
xlabel（'parametera'）;
ylabel（'迭代序列（51-100）'）;
混沌图
2，
折线图
程序：
a=0.5;x1=[];
x1
（1）=0.5;
fori=2:
20
x1（i）=a*sin（x1（i-1））;
end
n=1:
20;
subplot（2,2,1）,plot（n,x1）,title（'a=0.5,x0=0.5'）
a=2.0;x1=[];
x1
（1）=0.6;
fori=2:
20
x1（i）=a*sin（x1（i-1））;
end
n=1:
20;
subplot（2,2,2）,plot（n,x1）,title（'a=2.0,x0=0.6'）
a=3.0;x1=[];
x1
（1）=0.6;
fori=2:
20
x1（i）=a*sin（x1（i-1））;
end
n=1:
20;
subplot（2,2,3）,plot（n,x1）,title（'a=3.0,x0=0.6'）
a=3.5;x1=[];
x1
（1）=0.6;
fori=2:
20
x1（i）=a*sin（x1（i-1））;
end
n=1:
20;
subplot（2,2,4）,plot（n,x1）,title（'a=3.5,x0=0.6'）
（2）蛛网图法
x1=[];a=0.5;
x=-1.2:
0.01:
1.2;
y=a*sin（x）;subplot（2,2,1）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a.*sin（x（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=0.5',’x=0.2’）
end
>>x1=[];a=1.5;
x=-0.2:
0.01:
1.2;
y=a.*sin（x）;subplot（2,2,2）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.2;
fori=2:
50
x1（i）=a*sin（x（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=1.5',’x=0.2’）
end
>>x1=[];a=2.5;
x=-0.2:
0.01:
1.2;
y=a.*sin（x（i-1））;subplot（2,2,1）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.1;
fori=2:
50
x1（i）=a*sin（x（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=2.5',’x=0.1’）
end
>>x1=[];a=3;
x=-0.2:
0.01:
1.2;
y=a.*sin（x（i-1））;subplot（2,2,2）
plot（x,y）,holdon
ezplot（'x',[-0.2,1.2]）
x1
（1）=0.1;
fori=2:
50
x1（i）=a*sin（x（i-1））;
plot（[x1（i-1）,x1（i-1）],[x1（i-1）,x1（i）]）;
plot（[x1（i-1）,x1（i）],[x1（i）,x1（i）]）;
title（'a=3',’x=0.1’）
end
3）混沌图法
建立m文件（fl.m）
functionroot=fl（x,a）
x=[];x
（1）=0.2;
fori=2:
100
x（i）=a*sin（x（i-1））;
end
root=x;
程序：
clf
x=[];
holdon;
fora=2:
0.01:
4
root=fl（x,a）;
plot（a.*ones（size（root（51:
100）））,root（51:
100）,'.'）
end
xlabel（'parametera'）;
ylabel（'迭代序列（51-100）'）;
混沌图：
应用实验（或综合实验）
一、实验内容
1．油价与船速的优化问题
油价的上涨，将影响大型海船确定合理的航行速度，以优化航行收入。
直观地，油耗的多少直接影响船速的快慢，因而直接影响航行时间的长短，进而影响支付船员人工费用数量。
过去有一些经验表明：
（1）油耗正比于船速的立方；
（2）最省油航速的基础上改变20%的速度；则引起50%的油耗的变化。
作为一个例子：
某中型海船，每天油耗40吨，减少20%的航速，省油50%、即20吨。
每吨油价250美元，由此每天减少耗油费用5000美元，而航行时间的增加将增加对船员支付的费用的增加，如何最优化?
算例：
航程L=1536海里，标准最省油航速20节，油耗每天50吨，航行时间8天。
最低航速10节，本次航行总收入为84600美元。
油价250美元/吨，日固定开支1000美元。
试确定最佳航速。
解1）取航速为自变量设为：
x;
2）则油耗的话费为：
0.5/（0.8^3）*50*（x/20）^3*250*（1536/（x*24））=97.656*x^2
3）船员的日付费用即日固定开销：
升级会员 