当前位置:文档之家› 常微分方程数值解实验

常微分方程数值解实验

X=dsolve(‘f1’,’f2’,…) 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求 出通解,如果有初始条件,则求出特解。
有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无 法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程 数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般 格式为:
Image
Image
如果微 分方程 由一个 或多个 高阶常微分方程给出,要得到该方程的数值解,可以将方程转换成一阶 常微分方程组。假设高阶常微分方程的一般形式为y( n) = f ( t, y, yʹ, ⋯,y( n - 1) ),而且函数y(t)的各阶导数初值为y(0),yʹ(0) ,…, y( n - 1) (0)可以选 择一组变量令: x1= y, x2 = yʹ,…, xn = y( n - 1) ,我们就可以把原高阶常微 分方程转换成下面的一阶常微分方程组形式: 而且初值x1(0)=y(0),x2(0)=yʹ(0),…,xn(0)=(0)。 转换以后就可以求原 高阶常微分方程的数值解了。 例2 求微分方程,,的数值解。 对方程进行变换,选择变量 (1) 建立自定义函数 用edit命令建立自定义函数名为f.m,内容为: function y =f(t,x) y=[x(2);x(3);-t^2*x(2)*x(1)^2-t*x(1)*x(3)+exp(t*x(1))]; (2)调用对微分方程数值解ode45函数求解 用edit命令建立一个命令文件f2. m,内容为: >>x0=[2;0;0]; >>[t,y] =ode45(’f’,[0,10],x0);plot(t,y); >>figure; >>plot3(y(:,1),y(:,2), y(:,3))得
四阶Runge-Kutta法
以Euler法为基础,继续推广和处理,可得一、二、三、四阶RungeKutta格式,最常用的一种经典Runge-Kutta格式的具体形式如下: (n=0,1,2,…)
利用matlab求解一阶常微分方程的初值解
在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其 具体格式如下:
a==b
处的近似值yn(n=1,2,…),求解过程是顺着节点的次序一步步的向前推 进,即按递推公式由已知的y1,y2,y3… yi求出yi+1。
建立数值解法的一些途径
用差商代替导数
若步长较小,则有 固有公式: 此即Euler法。
使用数值积分
对微分方程两边积分得
故有公式: 和 此即改进的Euler法。
y0初始条件。 MATLAB 中有几个专门用于求解常微分方程的函数,它们的设计思
想基于Runge - Kutta方法,基本设计思想为:从改进的欧拉方法比欧拉 方法精度高的缘由着手,如果在区间[ x1 , xi+1 ] 多取几个点的斜率 值,然后求取平均值,则可以构造出精度更高的计算方法。 这些函数主 要包括:ode45、ode23、ode15s、ode113、ode23s、ode23t、ode23tb. 其中最常用的是函数ode45,该函数采用变步长四阶五阶Runge - Kutta 法求数值解,并采用自适应变步长的求解方法。ode23采用二阶三阶 Runge - Kutta法求数值解,与ode45类似,只是精度低一些。ode15s用来 求刚性方程组。 ode45函数的调用格式为: [ tout, yout ] = ode45 ( ’yprime’, [ t0, tf ] , y0) 其中yprime是表示f ( t, y)的M文件名, t0表示自变量的初始值, tf表示自 变量的终值, y0表示初始向量值. 输出向量tout表示节点( t0 ; t1 ; ⋯; tn) ,输出矩阵yout表示数值。
1411 120
186

5月13 日
5月14 日
5月15 日
5月16 日
5月17 日
5月18 日
5月19 日
5月20 日
5月21 日
5月22 日5 日
5月26 日
5月27 日
5月28 日
5月29
2304 2347 2370 2388 2405 2420 2434 2437 2444 2444 2456 2465 2490 2499 2504 2512 2514 2517
活带来了很大影响,附表给出了从2003年4月20日开始到6月23
日北京市按天公布的SARS传播数据,请建立微分方程模型对
SARS的传播过程进行预测,并对当时采取的隔离等措施进行评
估。另外,参考文献部分提供了一些和该问题相关的一些文献
资料可供大家参考。
参考文献
[1] 李贝, 徐海譞, 郭佳佳, 考虑自愈的SARS的传播模型, 工程 数学学报, (2003). [2] 谭欣欣, 冯恩民, 徐恭贤, 王宗涛, 修志龙, SARS流行病动力 学建模及其参数控制系统的研究, 工程数学学报, (2003). [3] 王议锋, 田一, 杨倩, 尚寿亭, 非典数学模型的建立与分析, 工程数学学报, (2003). [4] 肖红江, 吴彤, 李名科, 贺祖国, SARS传播的研究, 工程数学 学报, (2003). [5] 周义仓, 唐云, SARS传播预测的数学模型, 工程数学学报, (2003). [6] 邹宇庭, 郑晓练, 缪旭晖, 谭忠, SARS传播的数学原理及预 测与控制, 工程数学学报, (2003). [7] 李海龙, 任筱钰, 刘双, 用数学模型分析非典型肺炎预防和隔 离措施的有效性, 生物数学学报, (2004). [8] 刘云忠, 宣慧玉, 林国玺, SARS传染病数学建模及预防、控 制机理研究, 中国管理科学, (2004). [9] 刘云忠, 宣慧玉, 林国玺, SARS传染病数学建模及预测预防 控制机理研究, 中国工程科学, (2004).
一步算法,4,5阶
ode45
非刚性 Runge-Kutta 方法累积截断误差
大部分场合的 首选算法
一步算法,2,3阶
ode23
非刚性 Runge-Kutta 方法累积截断误差
使用于精度较 低的情形
多步法,Adams算 ode113非刚性 法,高低精度均可
达到
计算时间比 ode45短
适度刚 ode23t 性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法,Gear’s反 向
数值积分,精度中 等
若ode45失效 时,
可尝试使用
一步法,2阶 ode23s 刚性 Rosebrock算法,
低精度。
当精度较低 时,
计算时间比 ode15s短
odefx为显式常微分方程
中的 ,t为求解区间,要获得
问题在其他指定点
上的解,则令t=[t0,t1,t2,…](要求 单调),
附表:北京市疫情的数据
( 据:/Resource/Detail.asp?ResourceID=66070 )
已确诊
日 期 病例累

4月20 日
339
4月21 日
482
4月22 日
588
4月23 日
693
4月24 日
774
现有疑 死亡 治愈出 似病例 累计 院累计
402
18
33
610
25
43
666
28
46
782
35
55
863
39
64
4月25 日
4月26 日
4月27 日
4月28 日
4月29 日
4月30 日
5月01 日
5月02 日
5月03 日
5月04 日
5月05 日
5月06 日
5月07 日
5月08 日
5月09 日
5月10 日
5月11 日
5月12
877
988 1114 1199 1347 1440 1553 1636 1741 1803 1897 1960 2049 2136 2177 2227 2265
954
42
73
1093
48
76
1255
56
78
1275
59
78
1358
66
83
1408
75
90
1415
82
100
1468
91
109
1493
96
115
1537 100
118
1510 103
121
1523 107
134
1514 110
141
1486 112
152
1425 114
168
1397 116
175
x2 ( t)
用刚性方程求解函数可以快速求出该方程的数值解,并且画出两个状态 变量的时间曲线。x1 ( t) 曲线变化比较平滑, x2 ( t) 曲线变化在某些点上 较快.
6.3 综合实验任务
SARS的传播(根据CUMCM 2003 A改编)
SARS(Severe Acute Respiratory Syndrome,严重急性呼吸 道综合症, 俗称:非典型肺炎)是21世纪第一个在世界范围内 传播的传染病。SARS的爆发和蔓延给我国的经济发展和人民生
[10] 吕巍, 李海龙, 北京市SARS数学模型的建立与拟合, 辽宁 师范大学学报(自然科学版), (2004). [11] 王鑫, 郭玉翠, 用常微分方程模型分析预防和隔离措施对 SARS发病率的影响, 数学的实践与认识, (2004). [12] 刘双, 李海龙, 用差分方程模型模拟北京2003年SARS疫 情, 生物数学学报, (2006). [13] 莫小平, 一类SARS流行病模型的改进及分析, 数学的实践 与认识, (2007). [14] 杨玉华, 传染病模型的研究及应用, 数学的实践与认识, (2007). [15] 刘双, 王天辉, 北京2003年SARS疫情的数值模拟, 生物数 学学报, (2010).
相关主题