重庆大学数学实验3微分方程.docx
《重庆大学数学实验3微分方程.docx》由会员分享，可在线阅读，更多相关《重庆大学数学实验3微分方程.docx（22页珍藏版）》请在冰豆网上搜索。
重庆大学数学实验3微分方程
重庆大学
学生实验报告
实验课程名称数学实验
开课实验室DS1422
学院
学生姓名
开课时间
总成绩
教师签名
数学与统计学院制
开课学院、实验室：
数统学院DS1422实验时间：
课程
名称
数学实验
实验项目
名称
实验项目类型
验证
演示
综合
设计
其他
指导
教师
肖剑
成绩
实验目的
[1]归纳和学习求解常微分方程（组）的基本原理和方法；
[2]掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析；
[3]熟悉MATLAB软件关于微分方程求解的各种命令；
[4]通过范例学习建立微分方程方面的数学模型以及求解全过程；
通过该实验的学习，使学生掌握微分方程（组）求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
这对于学生深入理解微分、积分的数学概念，掌握数学的分析思维方法，熟悉处理大量的工程计算问题的方法是十分必要的。
基础实验
一、实验内容
1．微分方程及方程组的解析求解法；
2．微分方程及方程组的数值求解法——欧拉、欧拉改进算法；
3．直接使用MATLAB命令对微分方程（组）进行求解（包括解析解、数值解）；
4．利用图形对解的特征作定性分析；
5．建立微分方程方面的数学模型，并了解建立数学模型的全过程。
二、实验过程（一般应包括实验原理或问题分析，算法设计、程序、计算、图表等，实验结果及分析）
1．求微分方程的解析解,并画出它们的图形，
y’=y+2x,y（0）=1,0y’’+ycos（x）=0,y（0）=1,y’（0）=0;
（1）通解
y=dsolve（'Dy=y+2*x','x'）
y=
-2*x-2+exp（x）*C1
特解
y1=dsolve（'Dy=y+2*x','y（0）=1','x'）
y1=
-2*x-2+3*exp（x）
作图（y=3*exp（x）-2*x–2）
x=0:
0.02:
1;
y1=3*exp（x）-2*x-2;
plot（x,y1,'b-'）
（2）算法设计
通解y=dsolve（'D2y+y*cos（x）=0'）
y=
C19*exp（t*（-cos（x））^（1/2））+C20/exp（t*（-cos（x））^（1/2））
特解
y=dsolve（'D2y+y*cos（x）=0','y（0）=1,Dy（0）=0','x'）
结果y=
[emptysym]
分析：
此微分方程没有解析解.
该方程数值解的求解过程：
设y1=y,y2=y’则可得y1’=y2，y2’+y1cos（x）=0
首先建立M-文件：
functionf=fun（x,y）
f=[y
（2）;-y
（1）*cos（x）];
输入命令：
[x,y]=ode23（'fun',[0,1],[1,0]）
结果：
x=
0
0.0001
0.0005
0.0025
0.0125
0.0625
0.1541
0.2541
0.3541
0.4541
0.5541
0.6541
0.7541
0.8541
0.9541
1.0000
y=
1.00000
1.0000-0.0001
1.0000-0.0005
1.0000-0.0025
0.9999-0.0125
0.9980-0.0624
0.9882-0.1529
0.9680-0.2487
0.9386-0.3397
0.9003-0.4243
0.8540-0.5011
0.8004-0.5693
0.7404-0.6280
0.6751-0.6772
0.6053-0.7168
0.5721-0.7319
作图：
y1=y（:
1）;
>>y2=y（:
2）;
>>plot（x,y1,'b*',x,y2,'r'）,gtext（'y1'）,gtext（'y2'）
2．用向前欧拉公式和改进的欧拉公式求方程y’=y-2x/y,y（0）=1（0≤x≤1,h=0.1）的数值解，要求编写程序，并比较两种方法的计算结果，说明了什么问题？
程序：
向前欧拉和向后欧拉法:
x1
（1）=0;y1
（1）=1;y2
（1）=1;h=0.1;
fork=1;10
x1（k+1）=x1（k）+h;
y1（k+1）=y1（k）+h*（y1（k）-2*x1（k）/y1（k））
y2（k+1）=y2（k）+h*（y2（k+1）-2*x1（k+1）/y2（k+1））;
end
x1,y1,y2
ans=
10
y1=
1.0000
1.1000
1.0000
1.0000
0.9999
0.9980
0.9882
0.9680
0.9386
0.9003
0.8540
0.8004
0.7404
0.6751
0.6053
0.5721
x1=
00.1000
y1=
1.0000
1.1000
1.0000
1.0000
0.9999
0.9980
0.9882
0.9680
0.9386
0.9003
0.8540
0.8004
0.7404
0.6751
0.6053
0.5721
y2=
1.0000
251.0000
-0.0005
-0.0025
-0.0125
-0.0624
-0.1529
-0.2487
-0.3397
-0.4243
-0.5011
-0.5693
-0.6280
-0.6772
-0.7168
-0.7319
改进欧拉公式法：
算法设计：
x1
（1）=0;y1
（1）=1;h=0.1;
fork=1:
10
x1（k+1）=x1（k）+h;
y1（k+1）=y1（k）+1/2*h*（y1（k）-2*x1（k）/y1（k）+y1（k+1）-2*x1（k+1）/y1（k+1））;
end
x1,y1
结果：
x1=
Columns1through6
00.10000.20000.30000.40000.5000
Columns7through11
0.60000.70000.80000.90001.0000
y1=
1.0000
1.0957
1.1828
1.2624
1.3345
1.3989
1.4544
1.4994
1.5313
1.5462
1.5380
0.8004
0.7404
0.6751
0.6053
0.5721
3．Rossler微分方程组;
当固定参数b=2,c=4时，试讨论随参数a由小到大变化（如a∈（0,0.65））而方程解的变化情况，并且画出空间曲线图形，观察空间曲线是否形成混沌状？
建立M-file文件
functionf=fun（t,x）
f=[0,-1,-1;1,0.1,0;x（3）,0,-4]*x+[0,0,2]';
输入命令：
x0=[000.1];
[t,x]=ode45（'fun',[0,10],x0）;
plot（t,x（:
1）,'r',t,x（:
2）,'b-',t,x（:
3）,'--'）
pause
plot3（x（:
1）,x（:
2）,x（:
3））,gridon
当a=0.2时
当a=0.3时
当a=0.4时
当a=0.5
当a=0.6时
结论：
空间曲线不会形成混沌状
4.Apollo卫星的运动轨迹的绘制
则编写M文件函数：
functionf=fun（x,y）
u=1/82.45,u1=1-u,r1=sqrt（（y
（1）+u）^2+y（3）^2）,r2=sqrt（（y
（1）-u1）^2+y（3）^2）
f=[y
（2）;
2*y（4）+y
（1）-u1*（y
（1）+u）/（r1）^3-u*（y
（1）-u1）/（r2）^3;
y（4）;
-2*y
（2）+y（3）-u1*y（3）/（r1）^3-u*y（3）/（r2）^3]
输入：
[x,y]=ode45（'fun',[0,10],[1.2,0,0,-1.04935751]）
plot（y（:
1）,y（:
3）'）,gridon
作图：
应用实验（或综合实验）
一、实验内容
5．盐水的混合问题
一个圆柱形的容器，内装350升的均匀混合的盐水溶液。
如果纯水以每秒14升的速度从容器顶部流入，同时，容器内的混合的盐水以每秒10.5升的速度从容器底部流出。
开始时，容器内盐的含量为7千克。
求经过时间t后容器内盐的含量。
二、问题分析
已知：
水的密度为1kg/L,盐溶解度为36g（在假设
（1）中）
由题意可以知道，在流入到流出的过程中，由于混合在水中的盐含量是不同的，所以溶解于水中的盐的量每一时刻都是不同的。
而所求的这个瞬时量即为微分方程的解。
可计算出7kg盐所需要的溶剂为194L水。
因此，由混合液体积即可知开始时刻的7kg盐是完全溶于水中的,并且没有饱和。
所以，整个过程为食盐水被再次稀释的过程，则不会出现有盐析出现象。
三、数学模型的建立与求解（一般应包括模型、求解步骤或思路，程序放在后面的附录中）
假设：
1）温度对盐在水中的溶解度变化影响不大，
2）任意时刻容器内混合的、流出的盐水都均匀
3）水流入及盐水流出的速度均为匀速的。
4）假设注水时间为t、t时刻时容器内含盐量为P（t）、纯水流入容器的速度为v1,混合液流出的速度为v2。
计算过程如下：
1）t时刻容器内含盐量可表示为P（t）
2）纯水以v1=14L/s的速度流入容器，混合液以v2=10.5L/s的速度流出容器，则根据溶解度公式：
，可推出7kg食盐是完全溶解在水中的
3）根据假设1）、2），则可列出关于盐在水中的浓度的比例式，即为
（假设其中
足够小）
4）根据上述比例式即可得到相应的方程为：
而使
可得微分方程：
。
5）求解该微分方程当
时的曲线图像
四、实验结果及分析
编写M文件，并根据得出的结果作图
编写M文件：
y=dsolve（'Dy=-10.5*y/（350+3.5*t）','y（0）=7','t'）
得y=7000000/（t+100）^3
作图：
ezplot（'7000000/（t+100）^3',[0,100]）
总结与体会
通过本实验，我学会掌握微分方程（组）求解方法（解析法、欧拉法、梯度法、改进欧拉法等），对常微分方程的数值解法有一个初步了解，同时学会使用MATLAB软件求解微分方程的基本命令，学会建立微分方程方面的数学模型。
但对建立数学模型及matlab应用还不够灵活。
还需继续努力。
教师签名
年月日
升级会员 