常微分方程数值解实验
数值计算方法实验报告
数值计算方法实验报告一、实验介绍本次实验是关于数值计算方法的实验，旨在通过计算机模拟的方法，实现对于数值计算方法的掌握。
本次实验主要涉及到的内容包括数值微积分、线性方程组的求解、插值与拟合、常微分方程的数值解等。
二、实验内容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. 化学反应动力学研究：化学反应动力学研究需要求解涉及反应速率、物质浓度等的微分方程。
常微分方程数值解
dt (c x)2 (at y)2
由方程无法得到x(t), y(t)的解析解
需要用数值解法求解
0
Q(c,at)
P(x,y)
b
R(c,y )
cx
3
龙格—库塔方法的 MATLAB 实现
x(t) f (t, x), x(t0 ) x0 , x (x1, xn )T , f ( f1, fn )T
6
• clear,clf,shg • %Set the definied time • %ts=0:0.05:0.5; • %ts=0:0.1:1.6; • n=length(ts); • x0=[0 0]; • a=35;b=40;c=15; • opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); • [t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); • %a=35;b=40;c=15; • %opt1=odeset('RelTol',1e-6,'AbsTol',1e-9); • %[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c); • %exact solution x1=c • y1=a*t; • %output t,x(t),y(t) and draw x(t),y(t) • [t,x,y1] • plot(t,x),grid,gtext('x(t)','FontSize',16), • gtext('y(t)','FontSize',16),pause
常微分方程数值解
1
实例1 海上缉私
海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船 正以速度a向正北方向行驶，缉私艇立即以最大速度b(>a)前往 拦截。如果用雷达进行跟踪时，可保持缉私艇的速度方向始终 指向走私船。
【清华】实验4-常微分方程数值解(011813)
实验4 常微分方程数值解化学工程系分9班焦阳2009011813 【实验目的】1. 掌握用MATLAB软件求微分方程初值问题数值解的方法；2. 通过实例学习用微分方程模型解决简化的实际问题；3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】1.题目3：小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃烧用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
建立模型并进行分析：假设火箭在上升过程中，重力加速度g不随高度而变化，即固定g = 9.8m/s^2。
、(1)从火箭开始上升到引擎关闭：设火箭质量为m，高度为h，速度为v，加速度为a，阻力为f：，，ﭸ由牛顿第二定律可得：总ﭸ综上可得：；ﭸ；初值条件为：，；定义域为：。
根据常微分方程组的初值问题，在MATLAB中计算数值解，记，，， 。
通过解出微分方程的数值解，并进行绘图得到高度-时间曲线，速度-时间曲线，加速度-时间曲线如下：由MA度、速度t(s MATLAB 计算速度、加速度如下t(s) 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 计算得到的火箭度如下表：h(m)06.57326.4459.76106.5166.7240.2326.7425.7536.9659.8793.6937.81091.1254.1425.1604.1790.1983.2181.2384.的火箭从开始上h(m) 0 6.5737 26.444 59.762 106.57 166.79 240.27 326.72 425.79 536.99 659.8 793.63 937.85 1091.8 1254.7 1425.9 1604.8 1790.8 1983.1 2181.2 2384.5 开始上升到关闭引v(m/s)0 13.18926.57740.06253.53566.89 80.02192.829105.22117.11128.43139.14149.18158.55167.23175.22182.55189.22195.27200.75205.7 关闭引擎这段时间/s) 189 577 062 535 021 829 .22 .11 .43 .14 .18 .55 .23 .22 .55 .22 .27 .75 段时间内各时刻a(m/s^2)13.0571413.304513.4532813.4971913.4331313.2613112.9853412.6121912.1519511.6169311.0212710.38 9.7083329.0209048.3309057.6502496.9900526.3593395.7646125.2094884.694626各时刻的高^2)71404532871931313153421919569312733290490524905233961248862621 2592.4 210.18 4.22201222 2804.5 214.19 3.79432623 3020.6 217.79 3.41201724 3240.1 221.01 3.07303925 3462.7 223.92 2.7726326 3687.9 226.56 2.50440927 3915.6 228.97 2.26774828 4145.6 231.14 2.06332529 4377.8 233.11 1.88975930 4611.9 234.91 1.74334931 4847.7 236.57 1.6177932 5085 238.14 1.50617933 5323.8 239.61 1.40954434 5564.1 240.99 1.32933135 5805.8 242.28 1.2650336 6048.7 243.5 1.21393737 6292.9 244.68 1.17083938 6538.1 245.83 1.1302639 6784.5 246.96 1.09469840 7032 248.05 1.06634841 7280.5 249.1 1.04558942 7530.2 250.12 1.03075943 7780.9 251.14 1.01781844 8032.5 252.15 1.00237545 8285.1 253.16 0.98756846 8538.8 254.15 0.97632347 8793.4 255.12 0.96963948 9049 256.07 0.96629149 9305.6 257.03 0.96239250 9563.1 257.99 0.95272351 9821.5 258.95 0.9411852 10081 259.9 0.93372753 10341 260.83 0.93277654 10603 261.75 0.93633455 10865 262.67 0.9369556 11128 263.61 0.92584757 11392 264.54 0.91379358 11657 265.46 0.91062859 11923 266.35 0.9161260 12190 267.26 0.917011根据表格可以很容易得到：关闭引擎的瞬间，h=12190m，v=267.26m/s，a=0.917011m/s^2。
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解y′=f(x,y),x∈[x0,b]y(x0)=y0.定理1（利普希茨条件）若存在正数L，使得对任意，y1，y2，有|f(x,y1)−f(x,y2)|≤L|(y1−y2)|定理2（解存在性）①若函数f在方区域x∈[a,b],y∈R连续，②函数f关于y 满足利普希茨条件，则对任意x∈[a,b]，常微分方程存在唯一的连续可微数值解.两类问题：①单步法---计算下一个点的值yn+1只需要用到前面一个点的值yn②多步法---计算下一个点的值yn+1需要用到前面l个点的值yl1、欧拉法---下一个点的计算值等于前一个点的计算值加上步长乘以前一个点的函数值•具体过程一些批注：显式欧拉方程指下一步要计算的值，不在迭代方程中；隐式欧拉方程指下一步要计算的值，在迭代方程中。
怎么计算隐式欧拉方程----要借助显示欧拉迭代计算---一般用迭代法-----迭代---将微分方程在区间[xn,xn+1]进行积分，然后函数f进行近似，即可得到迭代方程-----迭代方程收敛性？由函数关于y满足利普希茨条件，可以推出迭代公式收敛。
•局部截断误差：假设前n步误差为0，我们计算第n+1步的误差，将次误差称为局部截断误差，且局部误差为O(hp+1)•p阶精度：由理论证明：若局部误差阶的时间复杂度为O(hp+1)，则整体误差阶为O(hp)我们称公式精度为p。
•显示欧拉法与隐式欧拉法•梯形方法----将显式欧拉迭代方程与隐式欧拉迭代方程做一下加权平均，构造的计算公式.•改进的欧拉方法---思想：因为梯形公式是隐式公式，将显式欧拉公式对下一步的计算值进行预估，用梯形公式对下一步的计算值进行校正.2、龙格-库塔方法思想：根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以前一个点的斜率；而这个斜率用该区间上的多个点的斜率的算数平均来逼近。
注意：怎么计算任意斜率Ki？第i个点的斜率Ki有微分方程可以算出f′=f(xn,yn)所以要算的f(xn,yn)值，由欧拉法即可算出, yn+1=yn+hf′•2阶-龙格-库塔方法----类似改进的欧拉法根据Lagrange中值定理，下一次的计算值可以用前一次的计算值加上h乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
一阶常微分方程初值问题的数值解fortron
一阶常微分方程初值问题的数值解是微分方程数值求解中的基础问题，对于工程、物理、生物等领域的科学计算和数值模拟具有重要意义。
本文将从常微分方程初值问题的数值解的基本原理和数值方法入手，详细介绍使用Fortran语言进行一阶常微分方程初值问题的数值解的实现过程和注意事项。
1. 常微分方程初值问题的数值解基本原理常微分方程初值问题的数值解是通过数值方法来逼近微分方程的解。
对于一阶常微分方程初值问题：dy/dx = f(x, y)y(x0) = y0其中，f(x, y)是给定的函数，y(x0) = y0是给定的初值条件。
求解该初值问题即是要找到一个函数y(x)近似地满足该微分方程，并且在点x = x0处与给定的初值条件相符。
2. 常微分方程初值问题的数值解的数值方法常见的数值方法包括欧拉方法、改进的欧拉方法、四阶龙格-库塔方法等。
其中，四阶龙格-库塔方法是最常用和最经典的数值方法之一。
该方法通过取若干个函数值点上的斜率的加权平均值来逼近微分方程的解，具有较高的数值精度和稳定性。
3. 使用Fortran语言实现一阶常微分方程初值问题的数值解Fortran是一种古老但经典的科学计算语言，以其高效的数值计算和科学工程计算而闻名。
下面将结合Fortran语言的特点，介绍如何使用Fortran语言实现一阶常微分方程初值问题的数值解。
（1）定义常微分方程的函数f(x, y)在Fortran程序中，首先需要定义常微分方程的函数f(x, y)，并将其定义为一个子程序或函数。
这里以一个简单的一阶线性常微分方程为例：f(x, y) = x + y则在Fortran程序中可以这样定义：```function f(x, y)real :: x, y, ff = x + yend function f```（2）实现四阶龙格-库塔方法在Fortran程序中，可以实现四阶龙格-库塔方法来数值解常微分方程初值问题。
具体做法是按照龙格-库塔方法的算法，在程序中编写相应的代码实现。
数值解常微分方程的方法和技巧
数值解常微分方程的方法和技巧在科学和工程领域，我们经常遇到一些复杂的常微分方程（Ordinary Differential Equations, ODEs），这些方程往往很难用解析方法得到精确解。
而数值解常微分方程的方法和技巧提供了一种有效的途径来近似求解这些方程。
本文将介绍一些常用的数值解ODEs的方法和技巧。
一、欧拉方法（Euler Method）欧拉方法是最简单的数值解ODEs的方法，它利用初始条件和微分方程的导数来计算下一个点的近似值。
具体来说，假设我们要求解的ODE为dy/dx = f(x, y)，其中f(x, y)是已知函数，初始条件为x0 = x(0)，y0 = y(0)。
欧拉方法的迭代公式为：y[i+1] = y[i] + h * f(x[i], y[i])其中，h是步长，x[i]表示第i个点的x坐标，y[i]表示对应的y坐标。
二、龙格-库塔方法（Runge-Kutta Method）龙格-库塔方法是一族常用的数值解ODEs方法，其基本思想是通过计算不同阶数的导数来提高求解的精度。
最常用的龙格-库塔方法是四阶龙格-库塔方法，也称为RK4方法。
它的迭代公式如下：k1 = h * f(x[i], y[i])k2 = h * f(x[i] + h/2, y[i] + k1/2)k3 = h * f(x[i] + h/2, y[i] + k2/2)k4 = h * f(x[i] + h, y[i] + k3)y[i+1] = y[i] + 1/6 * (k1 + 2*k2 + 2*k3 + k4)其中，k1、k2、k3、k4是中间变量，h是步长。
三、改进的欧拉方法（Improved Euler Method）改进的欧拉方法是对欧拉方法的改进，它通过使用导数的平均值来提高求解的精度。
其迭代公式为：k1 = h * f(x[i], y[i])k2 = h * f(x[i] + h, y[i] + k1)y[i+1] = y[i] + 1/2 * (k1 + k2)其中，k1、k2是中间变量，h是步长。
常微分方程数值解法
第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy微分方程的数值解：设微分方程的解y (x )的存在区间是[a,b ]，在[a,b ]内取一系列节点a= x 0< x 1<…< x n =b ,其中h k =x k+1-x k ；(一般采用等距节点，h=(b-a)/n 称为步长)。
在每个节点x k 求解函数y(x)的近似值:y k ≈y(x k )，这样y 0 , y 1 ,...,y n 称为微分方程的数值解。
用数值方法，求得f(x k )的近似值y k ，再用插值或拟合方法就求得y(x)的近似函数。
(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:⎪⎩⎪⎨⎧==00)(),(y x y y x f dx dy如果：(1) 在B y y A x x 00≤-≤≤,的矩形内),(y x f 是一个二元连续函数。
(2) ),(y x f 对于y 满足利普希茨条件，即2121y y L y x f y x f -≤-),(),(则在C x x 0≤≤上方程⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy的解存在且唯一，这里C=min((A-x 0),x 0+B/L),L 是利普希茨常数。
定义：任何一个一步方法可以写为),,(h y x h y y k k k 1k Φ+=+，其中),,(h y x k k Φ称为算法的增量函数。
收敛性定理：若一步方法满足： (1)是p 解的.(2) 增量函数),,(h y x k k Φ对于y 满足利普希茨条件.(3) 初始值y 0是精确的。
则),()()(p h O x y kh y =-kh =x -x 0,也就是有0x y y lim k x x kh 0h 0=--=→)((一)、主要算法 1．局部截断误差局部截断误差：当y(x k )是精确解时，由y(x k )按照数值方法计算出来的1~+k y 的误差y (x k+1)- 1~+k y 称为局部截断误差。
常微分方程数值解实验报告
常微分方程数值解实验报告学院：数学与信息科学专业：信息与计算科学：思义学号：201216524 课程：常微分方程数值解实验一：常微分方程的数值解法1、分别用Euler 法、改进的Euler 法（预报校正格式）和S —K 法求解初值问题。
（h=0.1）并与真解作比较。
⎩⎨⎧=++-=10(1y'）y x y 1.1实验代码：%欧拉法function [x,y]=naeuler(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值围，y0是初值，h 是步长 x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end%改进的欧拉法function [x,m,y]=naeuler2(dyfun,xspan,y0,h)%dyfun 是常微分方程，xspan 是x 的取值围，y0是初值，h 是步长。
%返回值x 为x 取值，m 为预报解，y 为校正解 x=xspan(1):h:xspan(2); y(1)=y0;m=zeros(length(x)-1,1); for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;end%四阶S—K法function [x,y]=rk(dyfun,xspan,y0,h)%dyfun是常微分方程，xspan是x的取值围，y0是初值，h是步长。
x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2);k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2);k4=feval(dyfun,x(n)+h,y(n)+h*k3);y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4);end%主程序x=[0:0.1:1];y=exp(-x)+x;dyfun=inline('-y+x+1');[x1,y1]=naeuler(dyfun,[0,1],1,0.1);[x2,m,y2]=naeuler2(dyfun,[0,1],1,0.1);[x3,y3]=rk(dyfun,[0,1],1,0.1);plot(x,y,'r',x1,y1,'+',x2,y2,'*',x3,y3,'o');xlabel('x');ylabel('y');legend('y为真解','y1为欧拉解','y2为改进欧拉解','y3为S—K解','Location','NorthWest');1.2实验结果：x 真解y 欧拉解y1 预报值m 校正值y2 S—K解y30.0 1.0000 1.0000 1.0000 1.00000.1 1.0048 1.0000 1.0000 1.0050 1.00480.2 1.0187 1.0100 1.0145 1.0190 1.01870.3 1.0408 1.0290 1.0371 1.0412 1.04080.4 1.0703 1.0561 1.0671 1.0708 1.07030.5 1.1065 1.0905 1.1037 1.1071 1.10650.6 1.1488 1.1314 1.1464 1.1494 1.14880.7 1.1966 1.1783 1.1945 1.1972 1.19660.8 1.2493 1.2305 1.2475 1.2500 1.24930.9 1.3066 1.2874 1.3050 1.3072 1.30661.0 1.3679 1.3487 1.3665 1.3685 1.36792、选取一种理论上收敛但是不稳定的算法对问题1进行计算，并与真解作比较。
常微分方程的数值解法实验报告
常微分方程的数值解法实验报告实验报告：常微分方程的数值解法摘要：常微分方程（ODE）是描述动力学系统中物理量随时间变化的数学方程，广泛应用于自然科学和工程领域。
然而，对于一些复杂的非线性ODE，很难找到解析解。
因此，我们需要数值解法来求解这些方程。
本实验报告将介绍四种常见的常微分方程数值解法，分别是欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法，并通过数值实验比较它们的精度和效率。
1.引言在实际问题中，许多物理量的变化规律可以由常微分方程描述。
然而，对于复杂的非线性ODE，很难找到解析解。
因此，为了解决这类问题，我们需要借助数值方法来求解。
2.方法本实验采用四种常见的常微分方程数值解法：欧拉法、改进的欧拉法、四阶龙格-库塔法和自适应步长的龙格-库塔法。
（1）欧拉法是最简单的数值解法，通过将微分方程转化为差分方程，使用离散的步长来近似微分方程。
（2）改进的欧拉法在欧拉法的基础上进行了改进，使用预估-校正的方法来提高精度。
（3）四阶龙格-库塔法是一种经典的数值解法，通过利用不同步长处的斜率来近似微分方程，具有较高的精度。
（4）自适应步长的龙格-库塔法是在四阶龙格-库塔法的基础上改进而来的，根据步长的大小自适应地选择不同的步长，同时保证精度和效率。
3.实验设计为了比较这四种数值解法的精度和效率，我们设计了两个实验。
实验一是求解一阶常微分方程：dy/dx = -2x，初始条件y(0) = 1，解析解为y = 1 - x^2、实验二是求解二阶常微分方程：d^2y/dx^2 + y = 0，初始条件y(0) = 0，dy/dx(0) = 1，解析解为y = sin(x)。
4.结果与分析实验一中，比较四种数值解法在不同步长下的近似解和解析解，计算其误差。
实验结果表明，四阶龙格-库塔法和自适应步长的龙格-库塔法具有较高的精度，而欧拉法和改进的欧拉法的精度较低。
实验二中，我们比较四种数值解法在不同步长下的近似解和解析解，并计算其误差。
微分方程数值解实验报告
微分方程数值解法课程设计报告班级：姓名：学号：成绩：2017年 6月 21 日摘要自然界与工程技术中的很多现象，可以归结为微分方程定解问题。
其中，常微分方程求解是微分方程的重要基础内容。
但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。
，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge—Kutta方法、Adams法以及椭圆型方程、抛物型方程的有限差分方法等，通过具体的算例，结合MATLAB求解画图，初步给出了一般常微分方程数值解法的求解过程。
同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。
关键词：微分方程数值解、MATLAB目录摘要 (2)目录 (3)第一章常微分方程数值解法的基本思想与原理 (4)1.1常微分方程数值解法的基本思路 (4)1.2用matlab编写源程序 (4)1.3常微分方程数值解法应用举例及结果 (5)第二章常系数扩散方程的经典差分格式的基本思想与原理 (6)2.1常系数扩散方程的经典差分格式的基本思路 (6)2.2 用matlab编写源程序 (7)2.3常系数扩散方程的经典差分格式的应用举例及结果 (8)第三章椭圆型方程的五点差分格式的基本思想与原理 (10)3.1椭圆型方程的五点差分格式的基本思路 (10)3.2 用matlab编写源程序 (10)3.3椭圆型方程的五点差分格式的应用举例及结果 (12)第四章总结 (12)参考文献 (12)第一章常微分方程数值解法的基本思想与原理1.1常微分方程数值解法的基本思路常微分方程数值解法(numerical methods forordinary differential equations)计算数学的一个分支.是解常微分方程各类定解问题的数值方法.现有的解析方法只能用于求解一些特殊类型的定解问题，实用上许多很有价值的常微分方程的解不能用初等函数来表示，常常需要求其数值解.所谓数值解，是指在求解区间内一系列离散点处给出真解的近似值.这就促成了数值方法的产生与发展.1.2用matlab编写源程序龙格库塔法：M文件:function dx=Lorenz(t,x)%r=28,sigma=10,b=8/3dx=[-10*(x(1)-x(2));-x(1)*x(3)+28*x(1)-x(2);x(1)*x(2)-8*x(3)/3];运行程序：x0=[1,1,1];[t,y]=ode45('Lorenz',[0,100],x0);subplot(2,1,1) %两行一列的图第一个plot(t,y(:,3))xlabel('time');ylabel('z');%画z-t图像subplot(2,2,3) %两行两列的图第三个plot(y(:,1),y(:,2))xlabel('x');ylabel('y'); %画x-y图像subplot(2,2,4)plot3(y(:,1),y(:,2),y(:,3))xlabel('x');ylabel('y');zlabel('z');%画xyz图像欧拉法：h=0.010;a=16;b=4;c=49.52;x=5;y=10;z=10;Y=[];for i=1:800x1=x+h*a*(y-x);y1=y+h*(c*x-x*z-y);z1=z+h*(x*y-b*z);x=x1;y=y1;z=z1;Y(i,:)=[x y z];endplot3(Y(:,1),Y(:,2),Y(:,3));1.3常微分方程数值解法的应用举例及结果应用举例：⎪⎪⎪⎩⎪⎪⎪⎨⎧-=--=--=)()()()()()()()()())()(()(t bz t y t x dt t dz t z t x t y t rx dt t dy t y t x a dt t dx a=10,b=8/3,0<r<+∞,当1<r<24.74时，Lorenz 方程有两个稳定的不动点c()1(-r b ，)1(-r b ，r-1)和c '（-)1(-r b ，-)1(-r b ，r-1），一个稳定的不动点0=（0，0，0），当r>24.74时，c 和c '都变成不稳定的，此时存在混沌和奇怪吸引子。
常微分方程初值问题的数值解法
常微分方程初值问题的数值解法在实际应用中，对于某些微分方程，我们并不能直接给出其解析解，需要通过数值方法来求得其近似解，以便更好地理解和掌握现象的本质。
常微分方程初值问题（IVP）即为一种最常见的微分方程求解问题，其求解方法有多种，本文将对常微分方程初值问题的数值解法进行较为详细的介绍。
一、欧拉法欧拉法是最基本的一种数值解法，它采用泰勒级数展开并截断低阶项，从而获得一个差分方程近似求解。
具体来讲，设 t 为独立变量，y(t) 为函数 y 关于 t 的函数，方程为：$$y'(t) = f(t, y(t)), \qquad y(t_0) = y_0$$其中 f(t,y(t)) 为已知的函数，y(t_0) 为已知的初值。
将函数 y(t) 进行泰勒级数展开：$$y(t+h) = y(t) + hf(t, y(t)) + O(h^2)$$其中 h 表示步长，O(h^2) 表示其他高阶项。
为了使误差较小，一般取步长 h 尽可能小，于是我们可以用欧拉公式表示数值解：$$y_{n+1} = y_n + hf(t_n, y_n), \qquad y_0 = y(t_0)$$欧拉法的优点是容易理解和实现，但是由于截取低阶项且使用的单步法，所以误差较大，精度较低，在具体应用时需要慎重考虑。
二、龙格-库塔法龙格-库塔法（Runge-Kutta method）是一种多步法，比欧拉法更加精确。
龙格-库塔法的主要思想是使用不同的插值多项式来计算近似解，并且将时间步长分解，每次计算需要多次求解。
以下简要介绍二阶和四阶龙格-库塔法。
二阶龙格-库塔法将时间步长 h 分解成两步 h/2，得到近似解表达式：$$\begin{aligned} k_1 &= hf(t_n, y_n)\\ k_2 &= hf(t_n+h/2,y_n+k_1/2)\\ y_{n+1} &= y_n+k_2+O(h^3)\\ \end{aligned}$$四阶龙格-库塔法四阶龙格-库塔法是龙格-库塔法中应用最为广泛的一种方法，其需要计算的中间值较多，但是具有更高的精度。
- 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).