实验报告七 常微分方程初值问题的数值解法.docx
《实验报告七 常微分方程初值问题的数值解法.docx》由会员分享,可在线阅读,更多相关《实验报告七 常微分方程初值问题的数值解法.docx(12页珍藏版)》请在冰豆网上搜索。
实验报告七常微分方程初值问题的数值解法
资料范本
本资料为word版本,可以直接编辑和打印,感谢您的下载
实验报告七常微分方程初值问题的数值解法
地点:
__________________
时间:
__________________
说明:
本资料适用于约定双方经过谈判,协商而共同承认,共同遵守的责任与义务,仅供参考,文档可直接下载或修改,不需要的部分可直接删除,使用时请详细阅读内容
浙江大学城市学院实验报告
课程名称数值计算方法
实验项目名称常微分方程初值问题的数值解法
实验成绩指导老师(签名)日期2015/12/16
一.实验目的和要求
用Matlab软件掌握求微分方程数值解的欧拉方法和龙格-库塔方法;
通过实例学习用微分方程模型解决简化的实际问题。
二.实验内容和原理
编程题2-1要求写出Matlab源程序(m文件),并有适当的注释语句;分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab源程序和运行结果和结果的解释、算法的分析写在实验报告上。
编程
编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序,问题如下:
在区间内个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句。
Euler法y=euler(a,b,n,y0,f,f1,b1)
改进Euler法y=eulerpro(a,b,n,y0,f,f1,b1)
分析应用题
假设等分区间数,用欧拉法和改进欧拉法在区间内求解初值问题
并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度。
分析应用题
用以下三种不同的方法求下述微分方程的数值解,取
画出解的图形,与精确值比较并进行分析。
1)欧拉法;
2)改进欧拉法;
3)龙格-库塔方法;
分析应用题
考虑一个涉及到社会上与众不同的人的繁衍问题模型。
假设在时刻(单位为年),社会上有人口人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。
而固定比例为的所有其他的后代也是与众不同的人。
如果对所有人来说出生率假定为常数,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:
其中变量表示在时刻社会上与众不同的人的比例,表示在时刻人口中与众不同的人的数量。
1)假定和,当步长为年时,求从到解的近似值,并作出近似解的曲线图形。
2)精确求出微分方程的解,并将你当时在分题(b)中得到的结果与此时的精确值进行比较。
【MATLAB相关函数】
求微分方程的解析解及其数值的代入
dsolve(‘egn1’,‘egn2’,‘’)
subs(expr,{x,y,…},{x1,y1,…})
其中‘egn’表示第个方程,‘’表示微分方程中的自变量,默认时自变量为。
subs命令中的expr、x、y为符合型表达式,x、y分别用数值x1、x2代入。
>>symsxyz
>>subs('x+y+z',{x,y,z},{1,2,3})
ans=
6
>>symsx
>>subs('x^2',x,2)
ans=
4
s=dsolve(‘’,‘’,‘’)
ans=
>>symsx
>>subs(s,x,2)
ans=
-0.3721
右端函数的自动生成
f=inline(‘expr’,’var1’,‘var2’,……)
其中’expr’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为f(var1,var2,……)。
>>f=inline('x+3*y','x','y')
f=
Inlinefunction:
f(x,y)=x+3*y
>>f(2,3)
ans=
11
4,5阶龙格-库塔方法求解微分方程数值解
[t,x]=ode45(f,ts,x0,options)
其中f是由待解方程写成的m文件名;x0为函数的初值;t,x分别为输出的自变量和函数值(列向量),t的步长是程序根据误差限自动选定的。
若ts=[t0,t1,t2,…,tf],则输出在自变量指定值,等步长时用ts=t0:
k:
tf,输出在等分点;options用于设定误差限(可以缺省,缺省时设定为相对误差,绝对误差),程序为:
options=odeset(‘reltol’,rt,’abstol’,at),这里rt,at分别为设定的相对误差和绝对误差。
常用选项见下表。
例:
解微分方程
在命令窗口执行
=(‘’,‘’,‘’);
;
ans=
01.0000
0.0502 1.0490
0.10051.0959
0.15071.1408
3.85072.9503
3.90052.9672
3.95022.9839
4.00003.0006
plot(,,‘o-’,)%解函数图形表示
%不用输出变量,则直接输出图形
;
ans=
01.0000
1.0000 1.7321
2.00002.2361
3.00002.6458
4.00003.0006
三.操作方法与实验步骤(包括实验数据记录和处理)
2-1编程
编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab程序,问题如下:
在区间内个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句。
Euler法y=euler(a,b,n,y0,f,f1,b1)
改进Euler法y=eulerpro(a,b,n,y0,f,f1,b1)
Euler法
y=euler(a,b,n,y0,f,f1,b1)
y=zeros(1,n+1);
y
(1)=y0;
h=(b-a)/n;
x=a:
h:
b;
fori=1:
n;
y(i+1)=y(i)+h*f(x(i),y(i));
end
plot(x,y)
holdon
%求微分方程的精确解
x1=linspace(a,b,100);
'精确解为'
s=dsolve(f1,b1,'x')
symsx
y1=zeros(1,100);
for
i=1:
100
y1(i)=subs(s,x,x1(i));
end
plot(x1,y1,'r')
title('红色代表精确解')
改进Euler法
y=eulerpro(a,b,n,y0,f,f1,b1)
%求微分方程的数值解
y=zeros(1,n+1);
y
(1)=y0;
h=(b-a)/n;
x=a:
h:
b;
for
i=1:
n;
T1=f(x(i),y(i));
T2=f(x(i+1),y(i)+h*T1);
y(i+1)=y(i)+(h/2)*(T1+T2);
end
plot(x,y)
holdon
%求微分方程的精确解
x1=linspace(a,b,100);
'精确解为'
s=dsolve(f1,b1,'x')
symsx
y1=zeros(1,100);
fori=1:
100
y1(i)=subs(s,x,x1(i));
end
plot(x1,y1,'r')
title('红色代表精确解')
2-2分析应用题
假设等分区间数,用欧拉法和改进欧拉法在区间内求解初值问题
并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度。
(1)向前欧拉法
>>euler(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10')
ans=
精确解为
s=
20-10*exp(x)
ans=
1.0e+005*
Columns1through8
0.00010.00010.00010.00010.00010.00000.0000
0.0000
Columns9through16
-0.0000-0.0000-0.0001-0.0001-0.0001-0.0001-0.0002
-0.0002
Columns17through24
-0.0003-0.0003-0.0004-0.0004-0.0005-0.0005-0.0006
-0.0007
Columns25through32
-0.0008-0.0009-0.0010-0.0011-0.0012-0.0014-0.0015
-0.0017
Columns33through40
-0.0019-0.0021-0.0024-0.0026-0.0029-0.0032-0.0035
-0.0039
Columns41through48
-0.0043-0.0048-0.0053-0.0058-0.0064-0.0071-0.0078
-0.0086
Columns49through56
-0.0095-0.0105-0.0115-0.0127-0.0140-0.0154-0.0170
-0.0187
Columns57through64
-0.0206-0.0227-0.0250-0.0275-0.0302-0.0333-0.0366
-0.0403
Columns65through72
-0.0444-0.0488-0.0537-0.0591-0.0651-0.0716-0.0788
-0.0867
Columns73through80
-0.0954-0.1049-0.1154-0.1270-0.1397-0.1537-0.1691
-0.1860
Columns81through88
-0.2046-0.2251-0.2477-0.2724-0.2997-0.3297-0.3627
-0.3990
Columns89through96
-0.4389-0.4828-0.5311-0.5842-0.6427-0.7070-0.7777
-0.8555
Columns97through101
-0.9410-1.0352-1.1387-1.2526-1.3779
改进欧拉法
>>eulerpro(0,10,100,10,inline('y-20','x','y'),'Dy=y-20','y(0)=10')
ans=
精确解为
s=
20-10*exp(x)
ans=
1.0e+005*
Columns1through8
0.00010.00010.00010.00010.00010.00000.0000-0.0000
Columns9through16
-0.0000-0.0000-0.0001-0.0001-0.0001-0.0002-0.0002-0.0002
Columns17through24
-0.0003-0.0003-0.0004-0.0005-0.0005-0.0006-0.0007-0.0008
Columns25through32
-0.0009-0.0010-0.0011-0.0013-0.0014-0.0016-0.0018-0.0020
Columns33through40
-0.0022-0.0025-0.0028-0.0031-0.0034-0.0038-0.0042-0.0047
Columns41through48
-0.0052-0.0058-0.0064-0.0071-0.0079-0.0087-0.0097-0.0107
Columns49through56
-0.0119-0.0131-0.0145-0.0161-0.0178-0.0197-0.0218-0.0241
Columns57through64
-0.0266-0.0294-0.0325-0.0360-0.0398-0.0440-0.0486-0.0537
Columns65through72
-0.0594-0.0656-0.0726-0.0802-0.0886-0.0980-0.1083-0.1197
Columns73through80
-0.1323-0.1462-0.1615-0.1785-0.1973-0.2180-0.2409-0.2663
Columns81through88
-0.2942-0.3251-0.3593-0.3971-0.4388-0.4849-0.5358-0.5921
Columns89through96
-0.6543-0.7230-0.7989-0.8828-0.9755-1.0780-1.1912-1.3163
Columns97through101
-1.4545-1.6073-1.7760-1.9626-2.1686
改进欧拉法的精度比向前欧拉法更高。
2-3分析应用题
用以下三种不同的方法求下述微分方程的数值解,取
画出解的图形,与精确值比较并进行分析。
1)欧拉法;
2)改进欧拉法;
2-4分析应用题
考虑一个涉及到社会上与众不同的人的繁衍问题模型。
假设在时刻(单位为年),社会上有人口人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。
而固定比例为的所有其他的后代也是与众不同的人。
如果对所有人来说出生率假定为常数,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:
其中变量表示在时刻社会上与众不同的人的比例,表示在时刻人口中与众不同的人的数量。
1)假定和,当步长为年时,求从到解的近似值,并作出近似解的曲线图形。
2)精确求出微分方程的解,并将你当时在分题(b)中得到的结果与此时的精确值进行比较。
1)
>>
euler(0,50,50,0.01,inline('0.002-0.002*p','t','p'),'Dp=0.002-0.002*p','p(0)=0.0
1')
ans=
精确解为
s=
1-99/(100*exp(x/500))
ans=
Columns1through8
0.01000.01200.01400.01590.01790.01990.02180.0238
Columns9through16
0.02570.02770.02960.03160.03350.03540.03740.0393
Columns17through24
0.04120.04310.04500.04700.04890.05080.05270.0546
Columns25through32
0.05640.05830.06020.06210.06400.06580.06770.0696
Columns33through40
0.07140.07330.07510.07700.07880.08070.08250.0844
Columns41through48
0.08620.08800.08980.09170.09350.09530.09710.0989
Columns49through51
0.10070.10250.1043
(2)
>>dsolve('Dp=0.002-0.002*p','p(0)=0.01','t')
ans=1-99/(100*exp(t/500))
>>1-99/(100*exp(0.1))
ans=0.1042
与欧拉法求得的精确值0.1043差0,0001
四.实验结果与分析