实验一 常微分方程的数值解法
常微分方程初值问题数值解法
常微分方程初值问题的数值解法在自然科学、工程技术、经济和医学等领域中，常常会遇到一阶常微分方程初值问题：(,),,(),y f x y a x b y a y '=≤≤⎧⎨=⎩ (1) 此处f 为,x y 的已知函数，0y 是给定的初始值。
本章讨论该问题的数值解法，要求f 在区域{(,)|,}G x y a x b y =≤≤<∞内连续，并对y 满足Lipschitz 条件，从而初值问题(1)有唯一的连续可微解()y y x =，且它是适定的。
1 几个简单的数值积分法1.1 Euler 方法（1）向前Euler 公式（显式Euler 公式）10(,),0,1,2,,(),n n n n y y hf x y n y y a +=+=⎧⎨=⎩(2) 其中h 为步长。
由此便可由初值0y 逐步算出一阶常微分方程初值问题(1)的解()y y x =在节点12,,x x 处的近似值12,,y y 。
该公式的局部截断误差为2()O h ，是一阶方法。
（2）向后Euler 公式（隐式Euler 公式）1110(,),0,1,2,,(),n n n n y y hf x y n y y a +++=+=⎧⎨=⎩(3) 这是一个隐格式，也是一阶方法。
这类隐格式的计算比显格式困难，一般采用迭代法求解。
首先用向前Euler 公式提供迭代初值，然后迭代计算：(0)1(1)()111(,),(,),0,1,2,n n n n k k n n n n y y hf x y y y hf x y k +++++⎧=+⎨=+=⎩ (4)1.2 梯形方法1110[(,)(,)],2(),(0,1,2,)n n n n n n h y y f x y f x y y y a n +++⎧=++⎪⎨⎪=⎩= (5) 这也是一个隐格式，是二阶方法。
一般也采用迭代法求解。
迭代公式如下：(0)1(1)()111(,),[(,)(,)],0,1,2,2n n n n k k n n n n n n y y hf x y h y y f x y f x y k +++++⎧=+⎪⎨=++=⎪⎩ (6)1.3 改进的Euler 方法11110(,),[(,)(,)],0,1,2,,2(),n n n n n n n n n n y y hf x y h y y f x y f x y n y y a ++++⎧=+⎪⎪=++=⎨⎪=⎪⎩(7) 为了便于上机编程计算，(7)可改写为110(,),(,),0,1,2,,1(),2(),p n n n cn n p n p c y y hf x y y y hf x y n y y y y y a ++=+⎧⎪=+⎪⎪=⎨=+⎪⎪=⎪⎩(8) 该格式是显式，也是二阶方法。
常微分方程的数值解法
常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。
然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。
在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。
这些方法可以给出解的近似表达式，通常称为近似解析方法。
还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。
利用计算机解微分方程主要使用数值方法。
我们考虑一阶常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy在区间[a, b]上的解，其中f (x, y )为x, y 的已知函数，y 0为给定的初始值，将上述问题的精确解记为y (x )。
数值方法的基本思想是：在解的存在区间上取n + 1个节点b x x x x a n =<<<<= 210这里差i i i x x h -=+1，i = 0,1, …, n 称为由x i 到x i +1的步长。
这些h i 可以不相等，但一般取成相等的，这时na b h -=。
在这些节点上采用离散化方法，（通常用数值积分、微分。
泰勒展开等）将上述初值问题化成关于离散变量的相应问题。
把这个相应问题的解y n 作为y (x n )的近似值。
这样求得的y n 就是上述初值问题在节点x n 上的数值解。
一般说来，不同的离散化导致不同的方法。
§1 欧拉法与改进欧拉法 1.欧拉法1．对常微分方程初始问题(9.2))((9.1) ),(00⎪⎩⎪⎨⎧==y x y y x f dx dy用数值方法求解时，我们总是认为(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)的近似值。
常微分方程的数值解
f ( x, y1 ) f ( x, y2 ) L y1 y2
（其中 L 为 Lipschitz 常数）则初值问题（ 1 ）存 在唯一的连续解。
求问题（1）的数值解，就是要寻找解函数在一 系列离散节点x1 < x2 <……< xn < xn+1 上的近似 值y1, y 2,…,yn 。 为了计算方便，可取 xn=x0+nh，(n=0,1,2,…)， h称为步长。
(1)，(2)式称为初值问题，(3)式称为边值问题。 在实际应用中还经常需要求解常微分方程组：
f1 ( x, y1 , y2 ) y1 ( x0 ) y10 y1 (4) f 2 ( x, y1 , y2 ) y2 ( x0 ) y20 y2
本章主要研究问题(1)的数值解法，对(2)~(4)只 作简单介绍。
得 yn1 yn hf ( xn1 , yn1 )
上式称后退的Euler方法，又称隐式Euler方法。 可用迭代法求解
二、梯形方法 由
y( xn1 ) y( xn )
xn1 xn
f ( x, y( x))dx
利用梯形求积公式： x h x f ( x, y( x))dx 2 f ( xn , y( xn )) f ( xn1 , y( xn1 ))
常微分方程的数言 简单的数值方法 Runge-Kutta方法 一阶常微分方程组和高阶方程
引言
在高等数学中我们见过以下常微分方程：
y f ( x, y, y) a x b y f ( x, y ) a x b (2) (1) (1) y ( x ) y , y ( x ) y 0 0 0 0 y ( x0 ) y0 y f ( x, y, y) a x b (3) y(a) y0 , y(b) yn
常微分方程(组)的数值解法
刚性常微分方程组求解
function demo1 figure ode23s(@fun,[0,100],[0;1]) hold on, pause ode45(@fun,[0,100],[0;1]) %-------------------------------------------------------------------------function f=fun(x,y) dy1dx = 0.04*(1-y(1))-(1-y(2)).*y(1)+0.0001*(1-y(2)).^2; dy2dx = -1e4*dy1dx + 3000*(1-y(2)).^2; f = [dy1dx; dy2dx];
解算指令的使用方法
[T,Y]=ode45(@fun, TSPAN,Y0) 输出变量T为返回时间列向量；解矩阵Y的每一行对应于T的 一个元素，列数与求解变量数相等。
@fun为函数句柄，为根据待求解的ODE方程所编写的ode文
件（odefile）； TSPAN＝[T0 TFINAL]是微分系统y'＝F(t,y)的积分区间；Y0 为初始条件
2.3 常微分方程（组）的数值解法
知识要点
常微分方程初值问题---ode45,0de23
微分方程在化工模型中的应用
•间歇反应器的计算 •活塞流反应器的计算
•全混流反应器的动态模拟
•定态一维热传导问题
•逆流壁冷式固定床反应器一维模型
•固定床反应器的分散模型
Matlab常微分方程求解问题分类
边值问题：
ode解算指令的选择(2)
2.根据常微分方程组是否为刚性方程
y ' Ay b( x) y ( x0 ) y0
微分方程数值解法实验报告
微分方程数值解法实验报告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]。
常微分方程初值问题数值解法
0.4 1.3582 1.3416 0.9 1.7178 1.6733
0.5 1.4351 1.4142 1.0 1.7848 1.7321
7
初值问题（2.2）有解 y ，1按2这x 个解析式子
算出的准确值 y(x同n )近似值 一y起n 列在表9-1中，两者 相比较可以看出欧拉方法的精度很差.
17
所以，局部截断误差可理解为用方法（2.10）计算一步的 误差，也即公式（2.10）中用准确解y(x代) 替数值解产生
的公式误差.
根据定义，显然欧拉法的局部截断误差
Tn1 y( xn1) y( xn ) hf ( xn , y( xn ))
y(xn h) y(xn ) hy(xn )
y(2) n1
yn
hf
( xn1,
y (1) n1
).
11
如此反复进行，得
y (k 1) n1
yn
hf
( xn1,
y(k) n1
),
(k 0,1, ).
（2.6）
由于 f (x,对y) 满足y 利普希茨条件（1.3）. 由（2.6）减 （2.5）得
y (k 1) n 1
yn1
h
f
( xn1,
y(k) n 1
积分曲线上一点 (x的, y切)线斜率等于函数 值.
的f (x, y)
如果按函数 f (在x, y) 平x面y上建立一个方向场，那 么，积分曲线上每一点的切线方向均与方向场在该点的方 向相一致.
基于上述几何解释，从初始点 P0 (x出0 ,发y0，) 先依 方向场在该点的方向推进到 x 上x1一点 ，P然1 后再从 P1 依方向场的方向推进到 x 上x2一点 ，循P2此前进做出
常微分方程数值解法
欧拉方法
总结词
欧拉方法是常微分方程数值解法中最基础的方法之一，其基本思想是通过离散化时间点上的函数值来 逼近微分方程的解。
详细描述
欧拉方法基于微分方程的局部线性化，通过在时间点上逐步逼近微分方程的解，得到一系列离散点上 的近似值。该方法简单易行，但精度较低，适用于求解初值问题。
龙格-库塔方法
总结词
影响
数值解法的稳定性对计算结果的精度和可靠 性有重要影响。
判断方法
通过分析数值解法的迭代公式或离散化方法， 判断其是否具有稳定性和收敛性。
数值解法的收敛性
定义
数值解法的收敛性是指随着迭代次数的增加， 数值解逐渐接近于真实解的性质。
影响
数值解法的收敛性决定了计算结果的精度和 计算效率。
分类
根据收敛速度的快慢，可以分为线性收敛和 超线性收敛等。
判断方法
通过分析数值解法的迭代公式或离散化方法， 判断其是否具有收敛性。
误差分析
定义
误差分析是指对数值解法计算过程中 产生的误差进行定量分析和估计的过 程。
分类
误差可以分为舍入误差、截断误差和 初始误差等。
影响
误差分析对于提高计算精度和改进数 值解法具有重要意义。
分析方法
通过建立误差传递公式或误差估计公 式，对误差进行定量分析和估计。
生物学
生态学、生物种群动态和流行病传播 等问题可以通过常微分方程进行建模
和求解。
化学工程
化学反应动力学、化学工程流程模拟 等领域的问题可以通过常微分方程进 行描述和求解。
经济学
经济系统动态、金融市场模拟和预测 等问题可以通过常微分方程进行建模 和求解。
02 常微分方程的基本概念
常微分方程的定义
常微分方程模型及其数值解
Q(c,at)
P(x,y)
R(c,y )
0
y
x
c
例2 弱肉强食
问题 自然界中在同一环境下的两个种群之间存在着几种不同的生存方式，比如相互竞争，即争夺同样的食物资源，造成一个种群趋于灭绝，而另一个趋向环境资源容许的最大容量；或者相互依存，即彼此提供部分食物资源，二者和平共处，趋于一种平衡状态；再有一种关系可称之为弱肉强食，即某个种群甲靠丰富的自然资源生存，而另一种群乙靠捕食种群甲为生，种群甲称为食饵（Prey），种群乙为捕食者（Predator），二者组成食饵-捕食者系统。海洋中的食用鱼和软骨鱼（鲨鱼等）、美洲兔和山猫、落叶松和蚜虫等都是这种生存方式的典型。这样两个种群的数量是如何演变的呢？近百年来许多数学家和生态学家对这一系统进行了深入的研究，建立了一系列数学模型，本节介绍的是最初的、最简单的一个模型，它是意大利数学家Volterra在上个世纪20年代建立的。
0.00 0.40 0.80 1.20 1.60 2.00
0.00000 0.36085 0.51371 0.50961 0.45872 0.40419
0.00000 0.34483 0.48780 0.49180 0.44944 0.40000
0.00000 -0.01603 -0.02590 -0.01781 -0.00928 -0.00419
而
从而有: y(xn+1)-yn+1=O(h3)
2.4 Taylor展开方法
设y(x)是初值问题（4)的精确解, 利用Taylor展开式可得
称之为p阶Taylor展开方法. …… …… …… 因此，可建立节点处近似值yn满足的差分公式 其中
所以,此差分公式是p阶方法.
02
常微分方程初值问题的数值解法
常微分方程初值问题数值解法初值问题：即满足初值条件的常微分方程的解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乘以斜率；而这个斜率用区间上的端点和中点的斜率的算数平均来逼近。
常微分方程初值问题的数值积分法
y( xn1) y( xn ) hfh ( xn , y( xn )) Rh ( xn ),
并且当 h 0时，
yh,0
y(x0 )
h1 max xIh
Rh (x)
0,
(7.4.2)
则称(7.4.1) 式为初值问题(7.1.1)的一个相容近似 ，
或称此格式满足相容条件即(7.4.2)式。
若
Rh (xn ) O(h p1), yh,0 y(x0) O(h p )(h 0)
a x0 x1 L xN1 xN b,
令 hn xn1 xn，称为积分网格的步长。
用y0 , y1, y2 , , yN 表示精确解 y(x)在节点 x0 , x1, x2 , , xN 上函数值 y( x0 ), y( x1), y( x2 ), , y( xN )的近似值。
对给定的数值积分法，各个 yn 是按某一递推算法确 定的。若在计算 yn1 时只用到已求出的 y0, y1,L , yn中的 yn ，而无须使用其余值 y0, y1,L , yn1 中的任何一个， 则称此法为单步法，否则，称之为多步法。
xn
的右端积分中用梯形公式，则得
yn1
yn
h 2
[
f
(
xn
,
yn
)
f (xn1, yn1)],
n 0,1,L
称该递推公式为梯形方法。
梯形公式
b f (x)dx (b a) ( f (a) f (b))
a
2
梯形方法
yn1
yn
h 2
[
f
( xn
,
yn
)
f (xn1, yn1)],
n 0,1,L
7.2 几个简单的数值积分法
常微分方程数值解法
第八章 常微分方程的数值解法一．内容要点考虑一阶常微分方程初值问题:⎪⎩⎪⎨⎧==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进行计算，并与真解作比较。
第一章常微分方程初值问题数值解法
（1.2.3）
其中rn,k(t)为插值余项。 代到（1.2.2）式中得
u ( tn +1 ) = u ( tn ) +
舍去余项 并用uj代替u(tj)即得
∫
tn+1 tn
Ln , k ( t ) dt + ∫ t rn , k ( t ) dt
tn+1
n
(1.2.4 (1.2.5)
Rn , k = ∫
⎡ ∑ ⎣α u
j =0
j =0
αk ≠ 0
（1.2.1）
j n+ j
⎤ − hβ j f n + j ⎦ = 0（数值解满足的差分方程）
因此称（1.2.1）为多步法 或 k-步法。 又因为(1.2.1)关于 u n + j , f n + j 是线性的，所以称为线性多步法。 为使多步法的计算能够进行，除给定的初值u0 外，还要 知道附加初值u1,u2,…,uk-1 ，这可用其它方法计算。 若 β k = 0 则称(1.2.1)是显式的； 若 β k ≠ 0 则方法（1.2.1）是隐式的。 例如，一般线性二步法可写成：
f ( t , u ( t ) ) = Ln , k +1 ( t ) + rn , k +1 ( t )
其中rn,k+1(t)为插值余项。 同理即
un +1 = un + h ∑ bk +1i f ( tn −i +1 , un −i +1 )
i =0
k +1
其中
bk +1i
=∫ ∏
−1
j =0 j ≠i
0
k +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.结果与分析实验一中，比较四种数值解法在不同步长下的近似解和解析解，计算其误差。
实验结果表明，四阶龙格-库塔法和自适应步长的龙格-库塔法具有较高的精度，而欧拉法和改进的欧拉法的精度较低。
实验二中，我们比较四种数值解法在不同步长下的近似解和解析解，并计算其误差。
常微分方程的数值解法2010
（11）
要使近似公式(8)的局部截断误差为O(h3)，则应要求(10)和(11) 式前三项相同：
c1 c 2 1 c2a2 1 / 2 c 2 b 21 1 / 2
以上方程组有无穷多组解，如取c1=c2=1/2，a2=b21=1，近似公 式(8)即为改进的Euler公式：
常微分方程的数值解法
对一阶常微分方程的初值问题，其一般形式是
y f ( x, y ) y (a ) y0 a x b
(1)
在下面的讨论中，假定f(x,y)连续，且关于y满足李
普希兹(Lipschitz)条件，即存在常数L，使得
f ( x, y ) f ( x, y ) L y y
Y(I+1)=Y(I)+(F(X(I),Y(I))+F(X(I+1),Y(I)+H*F(X(I),Y(I))))*H/2
设 x i x 0 ih ( i 1, 2 ,3 , ); y i , z i 为节点上的近似解， 则有改进的Euler格式为
y i 1 y i hf ( x i , y i ) h f ( x i , y i ) f ( x i 1 , y i 1 ) y i 1 y i 2
则初值问题(1)的解必定存在且唯一。
常微分方程的数值解法
所谓数值解法，就是要求问题(1)的解在若干点：
a x 0 x 1 x n 1 x n b
处的近似值yi(i=0,1,2…n)的方法，yi称为问题(1)的数
值解。相邻两个节点的间距
h n x i 1 x i
(5 )
计算方法 第6章 常微分方程数值解
已知Euler格式 yn1 yn hf ( xn , yn )
h2 y( xn1 ) yn1 2 y''( xn )
即Euler格式具有一阶精度
如果令
y( xn1 ) y( xn1 ) 2h
y'( xn )
f ( xn , yn )
并假定 y( xn1 ) yn1, y( xn ) yn
常微分方程数值解
常微分方程的数值解法
§1 引 言 §2 欧拉方法 §3 龙格-库塔方法
2
§1 引 言
在工程和科学技术的实际问题中，常需要解常微 分方程。但常微分方程组中往往只有少数较简单和典 型的常微分方程（例如线性常系数常微分方程等）可 求出其解析解。对于变系数常微分方程的解析求解就 比较困难，而一般的非线性常微分方程就更不用说了。 在大多数情况下，常微分方程只能用近似法求解。这 种近似解法可分为两大类：一类是近似解析法，如级 数解法、逐次逼近法等；另一类则是数值解法，它给 出方程在一些离散点上的近似解。
yn
2
xn yn
令 h 0.1 将 x0 0, y0 1 代入Euler格式
步进计算结果见P106表5.1
第五章：常微分方程数值解
Euler值
y 1 2x
第五章：常微分方程数值解
Euler格式的误差分析
pn1
事实上Euler格式的每一步都存在误差，为了方便讨论y算( x)
d2x
dt 2 x(t
0
c
m )
x x
0
0 (t
t) 0
x(t ) x
0
0
5
常微分方程初值问题的的数值解法
本章讨论常微分方程初值问题的数值解法
2
考虑一阶常微分方程的初值问题
⎧ dy ⎪ = f ( x, y ) ⎨ dx ⎪ ⎩ y (a ) = y0
x ∈ [a, b]
只要 f (x, y) 在[a, b] × R1 上连续，且关于 y 满足 Lipschitz 条 件，即存在与 x, y 无关的常数 L 使对任意x∈[a, b] ，和y1， y2 ∈ R1 都有 | f ( x, y1) − f ( x, y2 ) | ≤ L| y1 − y2 | 在唯一解。 成立， 则上述问题存
⎧ ⎪ ⎨ ⎪ ⎩ y n +1 = yn + hf ( xn , yn ), h yn +1 = yn + [ f ( xn , yn ) + f ( xn +1 , y n +1 )] 2
改进的Euler方法：y0=1,
y1=y0+hf (x0, y0) =1.1, y1=1+01./2 ×[(1−2 ×0/1)+(1.1−2 ×0.1/1.1)] =1.095909, …… y11=…… y11=1.737869.
1 yn +1 = yn + h[ f ( xn , yn ) + f ( xn +1 , yn +1 )] 2
12
称之为梯形公式。这是一个隐式的计算公式，欲求的yn+1需 解一个方程。
3．截断误差
定义 在假设 yn = y(xn)，即第 n 步计算是精确的前提下，考 虑的截断误差 εn+1 = y(xn+1) − yn+1 称为局部截断误差
⎧ y n +1 = y n + k1 ⎨ ⎩k1 = hf ( xn ,y n )
第六章常微分方程的数值解法
第六章常微分方程的数值解法第六章常微分方程的数值解法在自然科学研究和工程技术领域中，常常会遇到常微分方程的求解问题。
传统的数学分析方法仅能给出一些简单的、常系数的、经典的线性方程的解析表达式，不能处理复杂的、变系数的、非线性方程，对于这些方面的问题，只能求诸于近似解法和数值解法。
而且在许多实际问题中，确确实实并不总是需要精确的解析解，往往只需获得近似的解或者解在若干个点上的数值即可。
在高等数学课程中介绍过的级数解法和逐步逼近法，能够给出解的近似表达式，这一类方法称为近似解法。
还有一类方法是通过计算机来求解微分方程的数值解，给出解在一些离散点上的近似值，这一类方法称作为数值方法。
本章主要介绍常微分方程初值问题的数值解法，包括Euler 方法、Runge-Kutta 方法、线性多步法以及微分方程组与高阶微分方程的数值解法。
同时，对于求解常微分方程的边值问题中比较常用的打靶法与有限差分法作了一个简单的介绍。
§1 基本概念1.1 常微分方程初值问题的一般提法常微分方程初值问题的一般提法是求解满足如下条件的函数，，b x a x y ≤≤)(=<<=α)(),(a y bx a y x f dxdy， (1.1) 其中),(y x f 是已知函数，α是给定的数值。
通常假定上面所给出的函数),(y x f 在给定的区域},),{(+∞<≤≤=yb x a y x D 上面满足如下条件：(1) 函数),(y x f 在区域D 上面连续；(2) 函数),(y x f 在区域D 上关于变量y 满足Lipschitz(李普希茨)条件：212121,),(),(y y b x a y y L y x f y x f ?≤≤?≤?，， (1.2)其中常数L 称为Lipschitz(李普希茨)常数。
由常微分方程的基本理论可以知道，假如(1.1)中的),(y x f 满足上面两个条件，则常微分方程初值问题(1.1)对于任意给定的初始值α都存在着唯一的解，，b x a x y ≤≤)(并且该唯一解在区间[a,b]上是连续可微的。
微分方程数值解实验报告
微分方程数值解法课程设计报告班级：姓名：学号：成绩：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)。
实验一
常微分方程的数值解法
一、 实验内容 1. 分别用欧拉法和改进的欧拉法（预报—校正格式）和四阶龙格库塔法求 解初值问题。 （h=0.1）
y ′ = −y + x + 1 ,0 ≤ x ≤ 1 (1) y 0 =1
并与真解作比较（列出表格，并画图） 2. 选取一种理论上收敛但是不稳定的算法对问题（1）进行计算，并与真解 作比较（列表，画图） 。 二、 实验结果分析
二 实验结果分析
通过这次实验我掌握了将得到的解进一步精确， 而且要学会比较这几种方法 的精确性，显然，四阶龙格库达比改进欧拉发精确，改进欧拉发比欧拉法精确。
程序结果：
0 1 0.0250000000000000 1.00030990808236 0.0500000000000000 1.00122941879688 0.0750000000000000 1.00274348398340 0.100000000000000 0.125000000000000 0.150000000000000 0.175000000000000 0.200000000000000 0.225000000000000 0.250000000000000 0.275000000000000 0.300000000000000 0.325000000000000 0.350000000000000 0.375000000000000 0.400000000000000 0.425000000000000 0.450000000000000 0.475000000000000 0.500000000000000 0.525000000000000 1.00483741833333 1.00749689930416 1.01070797154688 1.01445701892311 1.01873075361613 1.02351621605355 1.02880077891340 1.03457212180418 1.04081822141212 1.04752735143119 1.05468808618799 1.06228927773127 1.07032004691684 1.07876978334422 1.08762814863660 1.09688505571054 1.10653066070931 1.11655536294553
Euler('fun',0,1,1,10) 其中fun.m为
function z=fun(x,y) z=-y+x+1;
程序结果：
2、改进的欧拉法
编写 TranEuler.m 函数，代码如下
function TranEuler(fun,x0,y0,xn,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; h=(xn-x0)/N; for n=1:N x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n)); y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n))+feval(fun,x(n+1),z0)); end T=[x',y'] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x',y','O',x1',y2');
0.550000000000000 0.575000000000000 0.600000000000000 0.625000000000000 0.650000000000000 0.675000000000000 0.700000000000000 0.725000000000000 0.750000000000000 0.775000000000000 0.800000000000000 0.825000000000000 0.850000000000000 0.875000000000000 0.900000000000000 0.925000000000000 0.950000000000000 0.975000000000000 1
xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓë¸Ä½øµÄÅ·À-·¨Çó½â¶Ô±È'); grid on; legend('¸Ä½øÅ·À-½â','Õæʵ½â');程Leabharlann 结果：3、4 阶 Rk 法
编写 4 阶 Rk 函数，代码如下
clear all; clc; tic fun=inline('-y+x+1','x','y'); [x,y]=ode45(fun,[0,1],1); B=[x,y] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x,y,'O',x1',y2'); xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓë4½×RK·¨Çó½â¶Ô±È'); grid on; legend('4½×RK·¨½â','Õæʵ½â'); toc
1、欧拉法
第一步：编写 Euler.m 函数，代码如下
function Euler(fun,x0,y0,xn,N) x=zeros(1,N+1); y=zeros(1,N+1); x(1)=x0; y(1)=y0; h=(xn-x0)/N; for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end T=[x',y'] y1 = dsolve('Dy=-y+x+1','y(0)=1','x') x1 = 0:0.1:1; y2=x1 + exp(-x1); A=[x1',y2'] plot(x',y','O',x1',y2'); xlabel('x'); ylabel('y(x)'); title('Õæʵ½âÓëÅ·À-·¨Çó½â¶Ô±È'); grid on; legend('Å·À-½â','Õæʵ½â');
1.12694980786900 1.13770486830921 1.14881163717622 1.16026142740887 1.17204577466010 1.18415642032450 1.19658530493382 1.20932456811006 1.22236655099527 1.23570378089426 1.24932896529859 1.26323499184411 1.27741493050958 1.29186201972077 1.30656966094317 1.32153141864356 1.33674102227941 1.35219235372536 1.36787944238047