常微分方程数值解实验
数值计算方法实验报告
数值计算方法实验报告一、实验介绍本次实验是关于数值计算方法的实验，旨在通过计算机模拟的方法，实现对于数值计算方法的掌握。
本次实验主要涉及到的内容包括数值微积分、线性方程组的求解、插值与拟合、常微分方程的数值解等。
二、实验内容1. 数值微积分数值微积分是通过计算机模拟的方法，实现对于微积分中的积分运算的近似求解。
本次实验中，我们将会使用梯形公式和辛普森公式对于一定区间上的函数进行积分求解，并比较不同公式的计算误差。
2. 线性方程组的求解线性方程组求解是数值计算领域中的重要内容。
本次实验中，我们将会使用高斯消元法、LU分解法等方法对于给定的线性方程组进行求解，并通过比较不同方法的计算效率和精度，进一步了解不同方法的优缺点。
3. 插值与拟合插值与拟合是数值计算中的另一个重要内容。
本次实验中，我们将会使用拉格朗日插值法和牛顿插值法对于给定的数据进行插值求解，并使用最小二乘法对于给定的函数进行拟合求解。
4. 常微分方程的数值解常微分方程的数值解是数值计算中的难点之一。
本次实验中，我们将会使用欧拉法和龙格-库塔法等方法对于给定的常微分方程进行数值解的求解，并比较不同方法的计算精度和效率。
三、实验结果通过本次实验，我们进一步加深了对于数值计算方法的理解和掌握。
在数值微积分方面，我们发现梯形公式和辛普森公式都能够有效地求解积分，但是辛普森公式的计算精度更高。
在线性方程组求解方面，我们发现LU分解法相对于高斯消元法具有更高的计算效率和更好的数值精度。
在插值与拟合方面，我们发现拉格朗日插值法和牛顿插值法都能够有效地进行插值求解，而最小二乘法则可以更好地进行函数拟合求解。
在常微分方程的数值解方面，我们发现欧拉法和龙格-库塔法都能够有效地进行数值解的求解，但是龙格-库塔法的数值精度更高。
四、实验总结本次实验通过对于数值计算方法的模拟实现，进一步加深了我们对于数值计算方法的理解和掌握。
在实验过程中，我们了解了数值微积分、线性方程组的求解、插值与拟合、常微分方程的数值解等多个方面的内容，在实践中进一步明确了不同方法的特点和优缺点，并可以通过比较不同方法的计算效率和数值精度来选择合适的数值计算方法。
数值计算实验报告-欧拉法常微分方程
数学与计算科学学院实验报告实验项目名称欧拉法解常微分方程所属课程名称数值计算实验类型验证型实验日期2012-6- 4班级隧道1002班学号201008020233姓名李彬彬成绩一、实验概述：【实验目的】 通过运用相关的数值计算软件，解决最基本的常微分方程的数值计算，并且能够熟练的运用这种方法。
【实验原理】 欧拉法1．对常微分方程初始问题(9.2))((9.1)),(00⎪⎩⎪⎨⎧==y x y y x f dxdy用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。
欧拉法是解初值问题的最简单的数值方法。
从（9.2）式由于y (x 0) = y 0已给定，因而可以算出),()('000y x f x y =设x 1 = h 充分小，则近似地有：),()(')()(00001y x f x y hx y x y =≈-(9.3)记 ,n ,,i x y y i i 10 )(== 从而我们可以取),(0001y x hf y y ==作为y (x 1)的近似值。
利用y 1及f (x 1, y 1)又可以算出y (x 2)的近似值：),(1112y x hf y y +=一般地，在任意点x n +1 = (n + 1)h 处y (x )的近似值由下式给出),(1n n n n y x hf y y +=+（9.4）这就是欧拉法的计算公式，h 称为步长。
不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。
还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。
【实验环境】Windows XP 环境下运行 NumericalAnalyse 软件二、实验内容：【实验方案】在区间[0,1]上以h=0.1为步长，分别用欧拉法与预估-校正法求初值问题y’=y-2x/y且 y|x=0 =1的数值解。
将上述方程输入到软件NumericalAnalyse中步骤如图选择常微分方程的数值解法。
微分方程数值解法实验报告
微分方程数值解法实验报告2篇微分方程数值解法实验报告（一）在实际科学与工程问题中，我们经常会遇到微分方程的求解。
然而，许多微分方程往往没有解析解，这就需要我们利用数值方法来获得近似解。
本篇实验报告将介绍两种常见的微分方程数值解法：欧拉方法和改进的欧拉方法。
一、欧拉方法欧拉方法是最简单的微分方程数值解法之一。
其基本原理为离散化微分方程，将微分方程中的导数用差商代替。
设要求解的微分方程为dy/dx = f(x, y)，步长为h，则可用以下公式进行递推计算：y_{n+1} = y_n +hf(x_n, y_n)二、算法实现为了对欧拉方法进行数值实验，我们以一阶线性常微分方程为例：dy/dx = x - y, y(0) = 1步骤如下：（1）选择合适的步长h和求解区间[a, b]，这里我们取h=0.1，[a, b] = [0, 1]；（2）初始化y_0 = 1；（3）利用欧拉方法递推计算y_{n+1} = y_n + 0.1(x_n - y_n)；（4）重复步骤（3），直到x_n = 1。
三、实验结果与讨论我们通过上述步骤得到了在[0, 1]上的近似解y(x)。
下图展示了欧拉方法求解的结果。
从图中可以看出，欧拉方法得到的近似解与精确解有一定的偏差。
这是因为欧拉方法只是通过递推计算得到的近似解，并没有考虑到更高阶的误差。
所以在需要高精度解时，欧拉方法并不理想。
四、改进的欧拉方法针对欧拉方法的不足，我们可以考虑使用改进的欧拉方法（也称为改进的欧拉-柯西方法）。
它是通过利用前后两个步长欧拉方法得到的结果作为差商的中间项，从而提高了求解精度。
一阶线性常微分方程的改进欧拉方法可以表示为：y_{n+1} = y_n + h(f(x_n, y_n) + f(x_{n+1}, y_n + hf(x_n,y_n)))/2五、算法实现与结果展示将改进的欧拉方法应用于前述的一阶线性常微分方程，我们同样选择h=0.1，[a, b] = [0, 1]。
实验报告七常微分方程初值问题的数值解法
浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1． 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格－库塔方法； 2． 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句；分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法； 3龙格－库塔方法；2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; 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‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格－库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名；x0为函数的初值；t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点；options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为：options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例：在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’；>>[],45(,[0,4],1)t y ode odefun =；ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =；[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下：在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8（2）改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法； 2改进欧拉法；2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为：其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
实验5常微分方程的数值解
实验5 常微分方程的数值解概要：将装满放射性废物的圆桶扔到水深300ft 的海底，圆桶体积55gal ，装满废料的桶重为527.436lbf ，在海中浮力为470.327lbf 。
此外，下沉时受到的阻力与速度成正比，比例系数为0.08lbf/s 。
实验发现当圆桶速度超过40ft/s 时，就会因与海底冲撞而破裂。
要求：（1）建立解决上述问题的微分方程模型（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司。
模型建立由牛顿第二定律可列出圆桶下沉速度的微分方程：dv G F f G F bv dt m m ----==其中G 为圆桶重量，F 为浮力，b 为下沉阻力与速度关系的比例系数。
换算到国际单位制，dept=300*0.3048=91.4400 海深（m ）ve=40*0.3048=12.1920 速度极限（超过就会使圆筒碰撞破裂）(m/s) G=527.436*0.4536*9.8＝2344.6 圆筒重量（N ） F=470.327*0.4536*9.8＝2090.7 浮力(N)m=527.436*0.4536＝239.24 圆筒质量（kg ） b=0.08*0.4536*9.8/0.3048＝1.1667 比例系数(Ns/m) 模型求解 一．求数值解Matlab 程序如下： sd.m:function dx=sd(t,x,G,F,m,b) dx=[(G-F-b*x)/m];%微分方程sddraw.m: clear;G=527.436*0.4536*9.8;%圆筒重量（N ） F=470.327*0.4536*9.8;%浮力(N)m=527.436*0.4536;%圆筒质量（kg ）b=0.08*0.4536*9.8/0.3048%比例系数(Ns/m) h=0.1;%所取时间点间隔ts=[0:h:2000];%粗略估计到时间2000 x0=0;%初始条件opt=odeset('reltol',1e-3,'abstol',1e-6);%相对误差1e-6，绝对误差1e-9 [t,x]=ode45(@sd,ts,x0,opt,G ,F,m,b);%使用5级4阶龙格—库塔公式计算 %[t,x]%输出t,x(t),y(t)plot(t,x,'-'),grid%输出v(t)的图形 xlabel('t'); ylabel('v(t)');%用辛普森公式对速度积分求出下沉深度 T=20;%估计20s 以内降到海底for i=0:2:10*T%作图时间间隔为0.2 y=x(1:(i+1)); k=length(y);a1=[y(2:2:k-1)];s1=sum(a1); a2=[y(3:2:k-1)];s2=sum(a2);z4((i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3;%辛普森公式求深度 endi=[0:2:10*T]; figure;de=300.*0.3048.*ones(5*T+1,1);%海深ve=40.*0.3048*[1 1];%速度极限值（超过就会使圆筒碰撞破裂）plot(x(i+1),z4',x(i+1),de,ve,[0 z4(5*T+1)]);%作出速度－深度图线，同时画出海深和速度要求grid;gtext('dept'),gtext('Vmax');xlabel('v');ylabel('dept(v)');figure;plot(i/10,z4');%作出时间－下降深度曲线grid;xlabel('t');ylabel('dept(t)');求解结果如下图：速度—时间曲线：可以看到经过足够长的时间后，如若桶没有落到海底，它的速度会趋于常值。
第9章 常微分方程初值问题数值解法
数值分析
第9章 常微分方程初值问题数值解法
《常微分方程》中介绍的微分方程主要有:
(1)变量可分离的方程 (2)一阶线性微分方程(贝努利方程) (3)可降阶的一类高阶方程 (4)二阶常系数齐次微分方程 (5)二阶常系数非齐次微分方程 (6)全微分方程 本章主要介绍一阶常微分方程初值问题的数值解法。
进一步: 令
y n1 y n
xn 1 xn
y n 1 y( x n 1 ) , y n y( x n )
f ( x , y( x ))dx h f ( x n , y n )
宽
9
高
实际上是矩形法
数值分析
第9章 常微分方程初值问题数值解法
(3)
用Taylor多项式近似并可估计误差
解决方法：有的可化为显格式，但有的不行 18
数值分析
第9章 常微分方程初值问题数值解法
与Euler法结合，形成迭代算法 ，对n 0，2， 1，
( yn0 )1 yn hf x n , yn ( k 1) h ( yn1 yn f x n , yn f x n1 , ynk )1 2
7
数值分析
第9章 常微分方程初值问题数值解法
建立数值解法的常用方法
建立微分方程数值解法,首先要将微分方程离散 化. 一般采用以下几种方法: (1) 用差商近似导数
dy yx yx x x dx x y
n 1 n n 1 n
n
,
n
进一步: 令
yn1 y( xn1 ) , yn y( xn )
由 x0 , y0 出发取解曲线 y y x 的切线（存在！），则斜率
实验4 常微分方程的数值解法
[内容] 1. 欧拉格式（或改进的欧拉格式），编写 相应的程序并能正确运行。 2. 经典四：先描述清楚问题。
实验4 常微分方程的数值解法 [要求]
1.程序的调试要耐心、细致；
2.语句应尽可能加注注释； 3.本次实验的各个程序(M文件)打包成压缩文件 （格式：学号姓名.RAR，如：200910119李娟.RAR）, 按时提交。
实验4 常微分方程的数值解法
[目的]
1.常微分方程差分算法的计算机实现；
2.进一步理解欧拉格式、改进的欧拉格式（预报—
校正系统）、龙格—库塔格式等算法，会运用这些方
法解决初步的常微分方程的求解问题； 3.进一步熟悉MATLAB数学软件的使用，锻炼程序 调试、排错的能力。
实验4 常微分方程的数值解法
常微分方程的数值解法与实际应用研究
常微分方程的数值解法与实际应用研究引言：常微分方程是数学中一种重要的数学工具，广泛应用于物理、经济、生物等领域的实际问题的数学建模。
在解析求解常微分方程存在困难或不可行的情况下，数值解法提供了一种有效的求解方法，并被广泛应用于实际问题的研究中。
本文将介绍常微分方程的数值解法以及一些实际应用的研究案例。
一、常微分方程的数值解法：1. 欧拉法：欧拉法是一种基础的数值解法，通过将微分方程离散化，近似得到方程的数值解。
欧拉法的基本思想是根据微分方程的导数信息进行近似计算，通过逐步迭代来逼近真实解。
但是欧拉法存在截断误差较大、收敛性较慢等问题。
2. 改进的欧拉法（改进欧拉法推导过程略）：为了解决欧拉法的问题，改进的欧拉法引入了更多的导数信息，改善了截断误差，并提高了算法的收敛速度。
改进欧拉法是一种相对简单而可靠的数值解法。
3. 四阶龙格-库塔法：四阶龙格-库塔法是常微分方程数值解法中最常用和最经典的一种方法。
通过多次迭代，四阶龙格-库塔法可以获得非常精确的数值解，具有较高的精度和稳定性。
二、常微分方程数值解法的实际应用研究：1. 建筑物的结构动力学分析：建筑物的结构动力学分析需要求解一些动力学常微分方程，例如考虑结构的振动和应力响应。
利用数值解法可以更好地模拟建筑物的振动情况，并对其结构进行安全性评估。
2. 生态系统模型分析：生态系统模型通常包含一系列描述物种数量和相互作用的微分方程。
数值解法可以提供对生态系统不同时间点上物种数量和相互作用的变化情况的模拟和预测。
这对于环境保护、物种保护以及生态系统可持续发展方面具有重要意义。
3. 电路模拟与分析：电路模拟与分析通常涉及电路中的电容、电感和电阻等元件，这些元件可以通过常微分方程进行建模。
数值解法可以提供电路中电压、电流等关键参数的模拟和分析，对电路设计和故障诊断具有重要帮助。
4. 化学反应动力学研究：化学反应动力学研究需要求解涉及反应速率、物质浓度等的微分方程。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无 法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程 数值解方面，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).