常微分方程数值解.docx
《常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《常微分方程数值解.docx（8页珍藏版）》请在冰点文库上搜索。
常微分方程数值解
淮海工学院
实验报告书
课程名称：
数学实验
实验名称：
常微分方程数值解
班级数学091
姓名：
耿萍学号：
090911107
日期：
2012.4.6地点数学实验室
指导教师：
曹卫平成绩：
数理科学系
1.实验目的：
1．用MATLAB软件掌握求微分方程数值解的方法。
2．通过实例学习用微分方程模型解决简化的实际问题。
2.实验内容：
（1）容器盛满水后，底端直径为
的小孔开启。
根据水力学知识，当水面高度为h时，水从小孔中流出的速度为v=0.6
（g为重力加速度，0.6为空口收缩系数）。
1）若容器为倒圆锥形，现测得容器高和上底面直径均为1.2m，小孔直径为3cm，问水从小孔中流完需要多少时间；2分钟时水面高度是多少。
2）若容器为倒葫芦形，现测得容器高1.2m，小孔直径3cm，由底端向上每隔0.1m测出容器的直径如下表所示，问水从小孔中流完需要多少时间；2分钟水面高度是多少。
x
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
1.1
1.2
D
0.03
0.05
0.08
0.14
0.19
0.33
0.45
0.68
0.98
1.10
1.20
1.13
1.00
（图略）
（2）一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。
已知河水流速为v1与船在静水中的速度v2之比为k。
1）建立小船航线的方程，求其解析解。
2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。
3.实验步骤：
（1）水面直径等于水深，设水深为h时，流量为
0.6（π/4）d^2*（gh）^0.5，0.6*（g*h）^（0.5）*π*（d0/2）^2*dt
=π/4*h^2*dh
则水深下降dh所需时间：
dt=-[（π/4）h^2*dh]/[0.6（π/4）d^2*（gh）^0.5]
=-[h^1.5*dh]/[0.6d^2*（g）^0.5]
水深由1.2m至0定积分得水从小孔流完的时间：
T（其中已知d=0.03m，g=9.8m*s（-2）
设两分钟（120S）后水深为Xm，由dt=-[（π/4）h^2*dh]/[0.6*（π/4）*d^2*（gh）^0.5]
=-[h^1.5*dh]/[0.6d^2*（g）^0.5]
则263.93-120=X^2.5/[1.5*d^2*（g）^0.5]
以d=0.03m，g=9.8m*s*-2代入上式得水深：
X
由知容器高1.2m，水深为h时，流量为0.6（π/4）d^2*（gh）^0.5，由于不同高度，倒葫芦形半径不同，用欧拉方程和龙格—库塔方法则水深下降dh所需时间：
dt=t（k+1）-t（k）
=-[（π/4）h^2*dh]/[0.6（π/4）d^2*（gh）^0.5]
=-[h^1.5*dh]/[0.6d^2*（g）^0.5]
然后利用循环
fork=1：
length（L），
t（k）=（（h（k+1）-h（k））*（π/4）*d（k）^2）/（0.6*（π/4）*d^2*（g（1.2-h（k）））^0.5），
T=sum（t）.可以求得水从小孔流完的总时间。
④设两分钟（120S）后水深为Xm，由S=0,利用条件，当120-s<0.0001时s=s+t（k），x（k）把d=0.03m，g=9.8m*s（-2）代入上式得水深：
X
程序
>>g=9.8;
d=0.03;
>>g=9.8;
d=0.03;
symsh
y=-h^（1.5）/（0.6*d^2*sqrt（2*g））;
T=int（y,h,1.2,0）;
t1=eval（T）
x=（（T-120）*（1.5*d^2*sqrt（2*g）））^（0.4）;
h1=eval（x）
t1=
263.9316
h1=
0.9416
clearall,clc;
g=9.8;
d=0.03;
k1=0;
k2=0;
h=0.4;
x
（1）=1.2;
forn=1:
1000k1=0.15*sqrt（g*（x（n）））*d^2/（-43.6359*x（n）^8+213.0457*x（n）^7-414.873*x（n）^6+410.2075*x（n）^5-218.8936*x（n）^4+62.553*x（n）^3-8.3215*x（n）^2+0.49619*x（n）+.014892）^2;k2=0.15*sqrt（g*（x（n））-h*k1）*d^2/（-43.6359*（x（n）-h*k1）^8+213.0457*（x（n）-h*k1）^7-414.873*（x（n）-h*k1）^6+410.2075*（x（n）-h*k1）^5-218.8936*（x（n）-h*k1）^4+62.553*（x（n）-h*k1）^3-8.3215*（x（n）-h*k1）^2+0.49619*（x（n）-h*k1）+.014892）^2;
x（n+1）=x（n）-h*（k1+k2）/2;
end
h2=x（300）
t=0:
h:
1000*h;
plot（t,x）;
axis（[0,400,0,1.21]）;
h2=
1.0278
（2）functiondx=key（t,x）
d=100;v1=1;v2=2;
s=sqrt（（d-x
（2））^2+（（v1*t）^2））;
dx=[v2*x
（1）/s;v2*（d-x
（2））/s];
ts=1:
1:
50;
x0=[0,0];
[t,x]=ode45（@key,ts,x0）;
[t,x]
plot（t,x）,grid,
gtext（'x（t）'）,gtext（'y（t）'）,pause
plot（x（:
1）,x（:
2））,grid,
gtext（'x'）,gtext（'y'）
ans=
2.000000
4.000003.9980
6.000007.9923
8.0000011.9801
10.0000015.9582
12.0000019.9225
14.0000023.8680
16.0000027.7883
18.0000031.6757
20.0000035.5209
22.0000039.3125
24.0000043.0372
26.0000046.6794
28.0000050.2218
30.0000053.6456
32.0000056.9316
34.0000060.0612
36.0000063.0182
38.0000065.7902
40.0000068.3690
42.0000070.7524
44.0000072.9424
46.0000074.9458
48.0000076.7726
50.0000078.4348
52.0000079.9454
54.0000081.3178
56.0000082.5652
58.0000083.6999
60.0000084.7333
62.0000085.6759
64.0000086.5371
66.0000087.3254
68.0000088.0482
70.0000088.7123
4.实验数据记录及分析（或程序及运行结果）：
通过这次试验，我学会了用MATLAB软件求微分方程数值解的方法，并且通过实例学习了微分方程模型解决简化的实际问题，我觉得学的知识只有运用到实际生活中才发挥了它的作用。
出师表
两汉：
诸葛亮
先帝创业未半而中道崩殂，今天下三分，益州疲弊，此诚危急存亡之秋也。
然侍卫之臣不懈于内，忠志之士忘身于外者，盖追先帝之殊遇，欲报之于陛下也。
诚宜开张圣听，以光先帝遗德，恢弘志士之气，不宜妄自菲薄，引喻失义，以塞忠谏之路也。
宫中府中，俱为一体；陟罚臧否，不宜异同。
若有作奸犯科及为忠善者，宜付有司论其刑赏，以昭陛下平明之理；不宜偏私，使内外异法也。
侍中、侍郎郭攸之、费祎、董允等，此皆良实，志虑忠纯，是以先帝简拔以遗陛下：
愚以为宫中之事，事无大小，悉以咨之，然后施行，必能裨补阙漏，有所广益。
将军向宠，性行淑均，晓畅军事，试用于昔日，先帝称之曰“能”，是以众议举宠为督：
愚以为营中之事，悉以咨之，必能使行阵和睦，优劣得所。
亲贤臣，远小人，此先汉所以兴隆也；亲小人，远贤臣，此后汉所以倾颓也。
先帝在时，每与臣论此事，未尝不叹息痛恨于桓、灵也。
侍中、尚书、长史、参军，此悉贞良死节之臣，愿陛下亲之、信之，则汉室之隆，可计日而待也
。
臣本布衣，躬耕于南阳，苟全性命于乱世，不求闻达于诸侯。
先帝不以臣卑鄙，猥自枉屈，三顾臣于草庐之中，咨臣以当世之事，由是感激，遂许先帝以驱驰。
后值倾覆，受任于败军之际，奉命于危难之间，尔来二十有一年矣。
先帝知臣谨慎，故临崩寄臣以大事也。
受命以来，夙夜忧叹，恐托付不效，以伤先帝之明；故五月渡泸，深入不毛。
今南方已定，兵甲已足，当奖率三军，北定中原，庶竭驽钝，攘除奸凶，兴复汉室，还于旧都。
此臣所以报先帝而忠陛下之职分也。
至于斟酌损益，进尽忠言，则攸之、祎、允之任也。
愿陛下托臣以讨贼兴复之效，不效，则治臣之罪，以告先帝之灵。
若无兴德之言，则责攸之、祎、允等之慢，以彰其咎；陛下亦宜自谋，以咨诹善道，察纳雅言，深追先帝遗诏。
臣不胜受恩感激。
今当远离，临表涕零，不知所言。
升级会员 