081数值计算方法—常微分方程(组)
常微分方程数值解法
用分段的折线逼近函数,此为 “折线法”而非“切线法”， 除第一个点是曲线上的切线，
其它都不是。
2、Euler方法的误差估计
1)局部截断误差。 在一步中产生的误差而非累积误差：
~
T x y y
n1
n1
n1
其中
~
y
是当
y
n
y(
x
)
n
（精确解！）时
n1
由Euler法求出的值，即y 无误差！ n
y y x , y x y 则得： h f
f ,
n1
n2
n
n
n1
n1
同样与Euler法结合，形成迭代算法，对n 0，1，2，
y y x , y 0 hf
n1
n
n
n
y y x , y x y
k 1
推出总体误差与步长的关系。
由微分方程解的存在唯一性,自然假定 （ f x，y）
充分光滑，或满足 Lipschitz条件：
f
x,ny源自xn
f
x
,
n
y
n
L
yxn
y n
第 n 步 的 总 体 截 断 误 差 记 为
en y
xn
y n
则 对 n 1 步：
e x y x y y y T y y ~ ~
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
常微分方程组的解法
常微分方程组的解法常微分方程组是由多个关于未知函数及其导数的方程组成的方程组，它是数学中的重要研究对象。
常微分方程组的解法可以分为解析解法和数值解法两种。
解析解法是指通过数学方法求出常微分方程组的解析表达式。
常微分方程组的解析解法主要包括分离变量法、一阶线性方程法、变量代换法、常数变易法、特殊函数法等。
其中，分离变量法是指将常微分方程组中的各个变量分离出来，然后对每个变量分别积分，最后得到常微分方程组的解析解。
一阶线性方程法是指将常微分方程组转化为一阶线性方程，然后通过求解一阶线性方程来得到常微分方程组的解析解。
变量代换法是指通过合适的变量代换将常微分方程组转化为更简单的形式，然后通过求解简化后的方程组得到常微分方程组的解析解。
常数变易法是指将常微分方程组中的常数作为未知量，然后通过求解常数得到常微分方程组的解析解。
特殊函数法是指通过特殊函数的性质求解常微分方程组，如指数函数、三角函数等。
数值解法是指通过计算机数值计算的方法求出常微分方程组的数值解。
常微分方程组的数值解法主要包括欧拉法、龙格-库塔法、变步长法等。
其中，欧拉法是一种简单的数值解法，它的基本思想是将常微分方程组的解曲线上的点离散化为一系列点，然后通过计算机逐步求解得到常微分方程组的数值解。
龙格-库塔法是一种高阶数值解法，它通过计算机采用多个不同的计算公式来逼近常微分方程组的解曲线，从而得到更为准确的数值解。
变步长法是一种自适应数值解法，它通过计算机根据误差大小自动调整步长大小，从而得到更为准确的数值解。
常微分方程组的解法包括解析解法和数值解法两种，每种方法都有其适用的范围和优缺点。
在实际应用中，需要根据具体问题的性质和求解要求选择合适的解法来求解常微分方程组。
常微分方程的数值解法
常微分方程的数值解法…………江南大学信计1203柯恒一、前言对于很多微分方程，我们很难求出解析解，这时我们需要采取数值手段求解。
在数值分析这门课中，老师讲到dy dt =f (t ,y )的数值解法，我们采用了欧拉格式，后退欧拉格式，梯形公式，改进欧拉格式及4阶龙格-库塔方法求解。
而老师没讲关于微分方程组和高阶微分方程的解法。
现在我来粗虐的说一下微分方程组及高阶微分方程的数值解，这里主要是借助matlab 中的相关函数进行计算并仿真出图像。
二、相关问题初值问题问题1，微分方程组问题描述：给定一个微分方程，并且告诉我们初始状态。
我们便可求出整个过程的解。
给定下面方程（L 表示省略号）：11122112112(,,,,)(,,,,)(,,,,)n n n n dy f x y y y dx dy f x y y y dx dy f x y y y dx ⎧=⎪⎪⎪=⎪⎨⎪⎪⎪=⎪⎩而且初始状态y ′i (x 0)已知记成y(0).问题解决：四阶龙格库塔方法：K i1=f i (x n ,y 1n ,y 2n ,……,y nn );K i2=f i (x n +h 2,y 1n +h 2∗k 11,y 2n +h 2∗K 21,……,y nn +h 2∗K n1) K i3=f i (x n +h ,y 1n +h ∗K 12,y 2n +h ∗K 22,……,y nn +h ∗K n2) K i4=f i (x n +h 2,y 1n +h 2∗K 13y 2n +h 2∗K 32,……,y nn +h 2∗K n3),y i,n+1=y i,n+h6∗(K i1+2∗K i2+2∗k i3+2∗K i4)通过龙格库塔方法，我们可以计算后面点的值，在matlab中调用ode45来实现四阶龙格库塔方法的调用。
应用举例：有一个同步地球卫星，现在加速到4km/s进行变轨，试分析变轨后卫星运动的轨迹。
已知地球的质量M=5.97e24,引力常量的值为6.672e-11问题建模：我们已经知道行星的运动是在一个平面上的。
常微分方程中的数值方法
常微分方程中的数值方法常微分方程是数学中的一个重要分支。
它主要研究的对象是随时间变化的函数。
在实际应用中，我们需要求解这些函数的解析解，但通常情况下，解析解并不容易得到，甚至是不可能得到。
因此，我们需要使用数值方法来求解这些函数的数值近似解。
在本文中，我们将介绍常微分方程中的数值方法。
一、欧拉法欧拉法是常微分方程数值解法中最基本的一种方法。
它是根据欧拉公式推导而来的。
具体地，我们可以将一阶常微分方程dy/dt=f(t,y)写成如下形式：y(t+h)=y(t)+hf(t,y(t))其中，h是步长，f(t,y)是t时刻y的导数。
欧拉法就是通过上面的公式进行逐步逼近，然后得到最终的数值解。
欧拉法的计算过程非常简单，但所得到的解可能会出现误差。
这是因为欧拉法忽略了f(t+h,y(t+h))和f(t,y(t))之间的变化。
因此，我们需要使用更为精确的数值方法来解决这个问题。
二、改进欧拉法为了解决欧拉法中的误差问题，我们可以使用改进欧拉法。
改进欧拉法又称作四阶龙格-库塔法。
它的基本思想是对欧拉法公式进行改进，以提高计算精度。
具体地，根据龙格-库塔公式，可将改进欧拉法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)改进欧拉法的计算过程比欧拉法要复杂些，但所得到的数值解比欧拉法更精确。
这种方法适用于一些特殊的问题，但在求解一些更为复杂的问题时，还需要使用其他的数值方法。
三、龙格-库塔法龙格-库塔法是求解常微分方程中数值解的常用方法之一。
它最常用的是四阶龙格-库塔法。
这种方法的基本思想是使用四个不同的斜率来计算数值解。
具体地，我们可以将四阶龙格-库塔法表示为：y(t+h)=y(t)+1/6(k1+2k2+2k3+k4)其中，k1=h*f(t,y)k2=h*f(t+h/2,y+k1/2)k3=h*f(t+h/2,y+k2/2)k4=h*f(t+h,y+k3)与改进欧拉法相比，龙格-库塔法的计算复杂度更高，但所得到的数值解更为精确。
微分方程组的数值求解方法
微分方程组的数值求解方法微分方程组数值求解方法微分方程组是数学中非常重要的一个分支，它描述了许多自然界和社会生活中的现象，例如电路的运行、天体的运行、生命体的生长等等。
我们需要对微分方程组进行求解，才能够得到它们的解析解，从而更好地理解和应用它们。
然而，大多数微分方程组不可能用解析法求解，因此，我们需要采用数值方法来求解微分方程组。
常见的微分方程组数值求解方法包括欧拉法、龙格库塔法和变步长法等。
下面，我们将逐一介绍它们的基本原理和优缺点。
一、欧拉法欧拉法是微分方程组数值求解方法中最简单的一种。
它的基本思想是将微分方程组中的各个变量离散化，然后根据微分方程组的导数计算每一步的值。
具体来讲，欧拉法的数值求解公式为：\begin{aligned} &x_{n+1}=x_n+hf_n(x_n,y_n,z_n),\\&y_{n+1}=y_n+hf_n(x_n,y_n,z_n),\\&z_{n+1}=z_n+hf_n(x_n,y_n,z_n), \end{aligned}其中，$x(t)$，$y(t)$，$z(t)$是微分方程组的解，$f_n(x_n,y_n,z_n)$是微分方程组导数在点$(x_n,y_n,z_n)$处的值，$h$为时间步长。
欧拉法的优点是简单易懂，方便实现，缺点是误差较大，计算不够精确。
因此，在实际应用中，往往需要采用更加精确的数值方法。
二、龙格库塔法龙格库塔法是微分方程组数值求解方法中比较常用的一种。
它的基本思想是通过多次计算微分方程组中的导数，以获得更加精确的数值解。
具体来讲，龙格库塔法的求解公式为：\begin{aligned}&k_{1x}=hf_n(x_n,y_n,z_n),k_{1y}=hf_n(x_n,y_n,z_n),k_{1z}=hf_n (x_n,y_n,z_n),\\&k_{2x}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{2},z_n+\frac{k_ {1z}}{2}),k_{2y}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{2},z_n+ \frac{k_{1z}}{2}),k_{2z}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{1y}}{ 2},z_n+\frac{k_{1z}}{2}),\\&k_{3x}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{2},z_n+\frac{k_ {2z}}{2}),k_{3y}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{2},z_n+ \frac{k_{2z}}{2}),k_{3z}=hf_n(x_n+\frac{h}{2},y_n+\frac{k_{2y}}{ 2},z_n+\frac{k_{2z}}{2}),\\&k_{4x}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3z}),k_{4y}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3z}),k_{4z}=hf_n(x_n+h,y_n+k_{3y},z_n+k_{3 z}),\\&x_{n+1}=x_n+\frac{k_{1x}}{6}+\frac{k_{2x}}{3}+\frac{k_{3x}}{ 3}+\frac{k_{4x}}{6},\\&y_{n+1}=y_n+\frac{k_{1y}}{6}+\frac{k_{2y}}{3}+\frac{k_{3y}}{ 3}+\frac{k_{4y}}{6},\\&z_{n+1}=z_n+\frac{k_{1z}}{6}+\frac{k_{2z}}{3}+\frac{k_{3z}}{ 3}+\frac{k_{4z}}{6}, \end{aligned}其中，$k_{1x}$，$k_{1y}$，$k_{1z}$，$k_{2x}$，$k_{2y}$，$k_{2z}$，$k_{3x}$，$k_{3y}$，$k_{3z}$，$k_{4x}$，$k_{4y}$，$k_{4z}$是微分方程组中导数的值。
常微分方程数值解法
第七章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。
另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。
因此，有必要探讨常微分方程初值问题的数值解法。
本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。
第一节 欧拉法求解常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(y x y y x f dxdy（1）的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<<<n x x x x 210 上的近似值 ,,,,,210n y y y y{}n y 称为问题的数值解，数值解所满足的离散方程统称为差分格式，1--=i i ix x h 称为步长，实用中常取定步长。
显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。
定理 设函数()y x f ,在区域+∞≤≤-∞≤≤y b x a D ,:上连续，且在区域D 内满足李普希兹(Lipschitz)条件，即存在正数L ，使得对于R 内任意两点()1,y x 与()2,y x ,恒有()()2121,,y y L y x f y x f -≤-则初值问题(1)的解()x y 存在并且唯一。
一、欧拉法(欧拉折线法)若将函数)xy 在点nx处的导数()n x y '用两点式代替， 即()hx y x y x y n n n )()(1-≈'+，再用n y 近似地代替()n x y ,则初值问题(1)变为⎩⎨⎧==++=+ ,2,1,0),()(001n x y y y x hf y y n n n n(2)(2)式就是著名的欧拉(Euler)公式。
以上方法称 为欧 拉法或欧拉折线法。
常微分方程组的数值解
dx xxi
x xi
由数值微分公式，我们有 向前差商公式
yi1 h
yi
f (xi , yi )
yi1 yi h f (xi , yi )
可以看到，给出初值，就可以用上式求出所有的 {yi}
求的是在一系列离散点列上，求未知函数y在这些点上的值的近 似。称为数值离散方法。
基本步骤如下：
① 对区间作分割：I : a x0 x1 xn b
(1)
y(a) y0
由常微分方程理论可知：只要上式中的函数f(x,y)在区域
G={a≤x≤b,－∞＜y＜∞}内连续，且关于 y 满足Lipschitz条件，即存 在与 x, y 无关的常数L，使
f (x, y1) f (x, y2 ) L y1 y2
对G内任意两个y1, y2 都成立，则方程(1) 的解y y(x)必定存在且唯一。
这类形式的方法也称为差分方法。
当假定 yi为准确值，即在 yi y(xi )的前提下来估计误差
y y ( xi1 )
i 1
，这种截断误差称为局部截断误差。
由（2）、（3）知Euler公式在 xi处的局部截断误差为：
Ri
y(xi1)
yi1
h2 2
y"( i ), (xi
i
xi1), (i
0,1,n 1)
① 收敛性问题 步长充分小时，所得到的数值解能否逼近问题的真解；
② 误差估计
③ 稳定性问题 舍入误差，在以后各步的计算中，是否会无限制扩大；
对于初值问题（1），先将其离散化，即把[a,b]区间n等分， 得各离散节点
xi a ih
(i 0,1,2,n)
其中h b a . n
设y yx为方程的解，则y(xi1 )在(xi , yi )处
常微分方程数值解
p
定义：一个方法的总体截断误差若为O h ， 则称之为 p阶方法。
一般地，方法的总体截断误差阶越高，精度 也越高。
误差分析表
Euler方法 局部截 断误差 向前Euler O(h2) 方法 向后Euler O(h2) 方法 O(h3) 梯形公式 总体截 断误差 O(h)
O(h) O(h2)
微分方程离散化常用方法
用差商代替微商 数值积分 Taylor展开
A 用差商代替微商 y xn1 y xn dy f ( xn , y( xn )) dx x n, y n xn1 xn
yn 1 yn 则 f xn , yn h yn 1 yn hf xn , yn n 0, 1, 2,
迭代收敛 条件
0<hL<1 0<hL<2
(L为Lip常数）
向后Euler 方法收敛条件与截 断误差
局部截断误差
T
O n 1
h
2
，
整体截断误差 O h (当 0 hL 1 时)
y y
n 1
k 1
k
n 1
h f
k 1
n 1
x
,y n 1
k
取左端点的矩形公式
则有Leabharlann xn 1xnf ( x, y )dx h f ( xn , yn )
yn 1 yn hf ( xn , yn ) (n 0,1,)
用yn 1 , yn 代替y ( xn 1 ), y ( xn ), 对右端积分采用 取右端点的矩形公式
则有
xn 1
其中hk=xk+1-xk , 如是等距节点h=(b-a)/n ,
数值分析中的常微分方程数值求解
数值分析中的常微分方程数值求解常微分方程是自然科学中一类最为普遍的数学模型，涉及到热力学、物理、化工等多个领域。
然而，解常微分方程并非易事。
尤其是当我们面对一些复杂、非线性、多维的方程组时，常微分方程数值求解成为了一个十分关键的问题。
因此，数值求解方法成为了常微分方程研究中的重要组成部分。
本文将介绍一些数值解常微分方程的常见方法和应用。
1. 一般线性方法一般线性方法（general linear methods）是经典的常微分方程数值解法之一。
它以一种特殊的形式给出步进公式：$$ y_{n+1}=\sum_{i=0}^{s-1}\alpha_i y_{n-i}+h\sum_{i=0}^{s-1}\beta_i f(t_{n-i},y_{n-i}) $$ 其中，$y_{n}$为第$n$步的项值，$f(t_n,y_n)$为时间$t_n$处函数$y(t)$的导数。
$\alpha_i$和$\beta_i$是常数，可以通过确定如下特征方程来选择：$$ \sum_{i=0}^{s-1}\alpha_i\ lambda^{i}=0,~(\lambda\in C) $$ 与此同时，也可以通过选择$\beta_i$来使方法达到一定的准确性和稳定性。
2. Runge-Kutta方法比一般线性方法更为流行的方法是Runge-Kutta方法。
通常附加一个或多个修正以获得更好的数值稳定性和误差控制。
第1阶Runge-Kutta方法仅使用导数$f(t_n,y_n)$估算下一个项的值：$$y_{n+1}=y_n+hf(t_n,y_n)$$ 许多高阶方法可以使用中间的“插值”来更准确地估计下一个步骤：$$y_{n+1}=y_n+h\sum_{i=1}^kb_ik_i$$$$k_i=f(t_n+c_ih,y_n+h\sum _{j=1}^{i-1}a_{ij}k_j)$$ $k_i$是第$i$台车的估计值，$a_{ij}$和$b_i$在经典Runge-Kutta方法和其他变体中具有不同的取值。
常微分方程解
8.2 改进Euler方法
对于(2)计算yn ,由于迭代工作量较大 , 一般只 迭代一次, 构成一类预估 校正算法, 即
( p) yn y n hf ( xn , y n ) h (c) ( p) y y [ f ( x , y ) f ( x , y n ` n n n n n )] (c) 并取yn yn 。
常微分方程数值解基本思想
数值方法的基本思想是：在解的存在区间上取 n + 1 个节点
a x0 x1 x2 xn b
这里差 hi xi 1 xi ，i = 0,1, …, n 称为由 xi 到 x i+1 的步长。这些 hi 可以 相等，但一般取成相等的，这时 h
常微分方程描写的物理现象
镭的衰变规律 单摆的运动 伯努利方程 RLC振荡电路 物理场计算
常微分方程
我们考虑一阶常微分方程初值问题
dy f ( x, y ) dx y( x ) y
在区间[a, b]上的解，其中 f (x, y) 为 x, y 的已知 函数，y0 为给定的初始值，将上述问题的精确 解记为 y(x)。
8.2 改进Euler方法
上式还常写成 yn yn (k k ) k hf ( xn , yn ) k hf ( x h, y k ) (n ,,,...) n n 该式称为改进 Euler方法, 亦可写成 h yn yn [ f ( xn yn ) f ( xn , yn hf ( xn , yn ))]
这就是隐式的 `Euler公式或向后Euler方法，它与显式 的不同在于，它每算一 步要解函数方程 (2)才能得到 yn 1。
常微分方程的数值解
欧拉方法的实现
确定步长和初始值
根据问题的性质和精度要求，选择合适的步长 和初始值。
迭代计算
根据欧拉方法的公式，迭代计算下一个点的值。
终止条件
当达到预设的迭代次数或误差范围时，停止迭代。
常微分方程的应用
总结词
常微分方程在自然科学、工程技术和社会科学等领域有广泛应用。
详细描述
在物理学中，常微分方程可以描述物体的运动规律、电磁波的传播等；在化学中，可以描述化学反应 的动力学过程；在社会科学中，可以用于研究人口增长、经济趋势等。此外，常微分方程还在控制工 程、航空航天等领域有广泛应用。
确定步长和初始值
在应用龙格-库塔方法之前，需要 选择合适的步长和初始值。步长 决定了迭代的精度，而初始值则 用于启动迭代过程。
编写迭代公式
根据选择的步长和初始值，编写 龙格-库塔方法的迭代公式。该公 式将使用已知的函数值和导数值 来计算下一步的函数值。
迭代求解
按照迭代公式进行迭代计算，直 到达到所需的精度或达到预设的 最大迭代次数。
欧拉方法的误差分析
截断误差
由于欧拉方法只使用了微分方程的一次项， 因此存在截断误差。
全局误差
全局误差是实际解与近似解之间的最大偏差。
局部误差
由于每一步都使用了上一步的结果，因此存 在局部误差。
稳定性
欧拉方法是稳定的，但步长和初始值的选择 会影响其稳定性和精度。
04 龙格-库塔方法
龙格-库塔方法的原理
常用的数值解法包括欧拉方法、龙格-库塔方法、改进的欧拉方法、预估 校正方法和步进法等。
常微分方程的数值解法
数值计算方法
都是一次的,则y称它, y是线, 性的, ,y否(n则) 称为非线性的。
在高等数学中，对于常微分方程的求解，给出 了一些典型方程求解析解的基本方法，如可分离变 量法、常系数齐次线性方程的解法、常系数非齐次 线性方程的解法等。但能求解的常微分方程仍然是 有限的，大多数的常微分方程是不可能给出解析解。 譬如
y x2 y2
这个一阶微分方程就不能用初等函数及其积 分来表达它的解。
再如，方程
y y
y
(0)
1
的解 y e x ,虽然有表可查,但对于表 上没有给出 e x 的值,仍需插值方法来
计算
从实际问题当中归纳出来的微分方程，通常主要依
靠数值解法来解决ቤተ መጻሕፍቲ ባይዱ本章主要讨论一阶常微分方程
初值问题
y f (x, y)
y
(
x0
)
y0
( 7.1 )
在区间a ≤ x ≤ b上的数值解法。
可以证明,如果函数在带形区域 R=a≤x≤b,
-∞＜y＜∞}内连续，且关于y满足李普希兹
(Lipschitz)条件，即存在常数L(它与x,y无关)使
f (x, y1) f (x, y2 ) L y1 y2
对R内任意两个 y1, y2 都成立,则方程( 7.1 )的解 y y(x) 在a, b上存在且唯一。
数值计算方法
常微分方程的数值解法
包含自变量、未知函数及未知函数的导数或微 分的方程称为微分方程。在微分方程中, 自变量的 个数只有一个, 称为常微分方程.。自变量的个数 为两个或两个以上的微分方程叫偏微分方程。微分 方程中出现的未知函数最高阶导数的阶数称为微分 方程的阶数。如果未知函数y及其各阶导数
第8章 常微分方程数值解法 本章主要内容： 1.欧拉法
第8章 常微分方程数值解法本章主要内容：1．欧拉法、改进欧拉法. 2．龙格-库塔法。
3．单步法的收敛性与稳定性。
重点、难点一、微分方程的数值解法在工程技术或自然科学中，我们会遇到的许多微分方程的问题，而我们只能对其中具有较简单形式的微分方程才能够求出它们的精确解。
对于大量的微分方程问题我们需要考虑求它们的满足一定精度要求的近似解的方法，称为微分方程的数值解法。
本章我们主要讨论常微分方程初值问题⎪⎩⎪⎨⎧==00)(),(yx y y x f dx dy的数值解法。
数值解法的基本思想是：在常微分方程初值问题解的存在区间[a,b]内，取n+1个节点a=x 0＜x 1＜…＜x N =b （其中差h n = x n –x n-1称为步长，一般取h 为常数，即等步长），在这些节点上把常微分方程的初值问题离散化为差分方程的相应问题，再求出这些点的上的差分方程值作为相应的微分方程的近似值（满足精度要求）。
二、欧拉法与改进欧拉法欧拉法与改进欧拉法是用数值积分方法对微分方程进行离散化的一种方法。
将常微分方程),(y x f y ='变为()*+=⎰++11))(,()()(n xn x n n dtt y t f x y x y1．欧拉法（欧拉折线法）欧拉法是求解常微分方程初值问题的一种最简单的数值解法。
欧拉法的基本思想：用左矩阵公式计算（＊）式右端积分，则得欧拉法的计算公式为：Nab h N n y x hf y y n n n n -=-=+=+)1,...,1,0(),(1 欧拉法局部截断误差11121)(2++++≤≤''=n n n n n x x y h R ξξ或简记为O （h 2）。
我们在计算时应注意欧拉法是一阶方法，计算误差较大。
欧拉法的几何意义：过点A 0（x 0，y 0），A 1（x 1，y 1），…，A n （x n ，y n ），斜率分别为f （x 0，y 0），f （x 1，y 1），…，f （x n ，y n ）所连接的一条折线，所以欧拉法亦称为欧拉折线法。
数值方法常微分方程
数值方法常微分方程数值方法是一种近似求解常微分方程（ODEs）的方法，它是通过将连续问题离散化为离散问题来实现的。
常微分方程是数学中常见的用于描述动态系统的工具，它描述了未知函数与其导数之间的关系。
求解常微分方程对于预测系统的行为和发展非常重要。
在许多现实问题中，解析求解常微分方程是非常困难甚至不可能的。
而数值方法则提供了一种近似求解常微分方程的有效和可行的途径。
数值方法基于将微分方程中的函数在离散的点上进行近似，通过计算函数的离散解来预测函数在给定时间和空间范围内的行为。
常用的数值方法包括欧拉方法、隐式欧拉方法、梯形规则、龙格-库塔方法（RK4）、多步法等。
在数值方法中，最简单的方法是欧拉方法（Euler method）。
该方法将微分方程中的导数用差分代替，通过迭代逼近函数的解。
该方法的基本思想是通过将微分方程近似为差分方程，在离散的时间点上计算函数的值。
欧拉方法的计算公式是：y[i+1]=y[i]+h*f(t[i],y[i])其中，y[i]是在时间点t[i]处的函数的近似值，h是时间步长，f(t[i],y[i])是在时间t[i]和函数y[i]处的导数值。
尽管欧拉方法是最简单的数值方法之一，但它有一些局限性。
首先，它对步长的选择非常敏感，步长选择过大或过小都可能导致数值解的不稳定性。
其次，欧拉方法的精度较低，由于使用了一阶近似，所以在一些情况下可能会产生较大的误差。
因此，为了提高数值解的精度和稳定性，我们需要使用更高阶的数值方法。
龙格-库塔方法（RK方法）是一种常用的、高阶的数值方法。
它是一系列积分方法的集合，其中RK4是最常用的方法之一、RK4可以通过使用连续的斜率来计算函数值的改变量，在四个时间点上进行计算。
该方法的计算公式为：k1=h*f(t[i],y[i])k2=h*f(t[i]+h/2,y[i]+k1/2)k3=h*f(t[i]+h/2,y[i]+k2/2)k4=h*f(t[i]+h,y[i]+k3)y[i+1]=y[i]+1/6*(k1+2*k2+2*k3+k4)与欧拉方法相比，RK4方法具有更高的精度。
第10章_常微分方程数值解
迭代法的收敛条件：
计
算 方
|
y y (k 1) i 1
(k ) i 1
|
h| 2
f
( xi1,
y(k) i 1
)
f
( xi1,
y (k 1) i 1
)
|
法 课
h 2
L|
y(k) i 1
y (k 1) i 1
|
件
其中L为Lipschitz常数。
当
0 h L 1
2
时迭代过程收敛。即迭代过程中误差没有扩大。
称为欧拉中点公式.容易看出,中点公式计算yi+1时,不仅需要yi的 值，还需要yi-1的值.
8
结束
凡是从已知(或已算出)的y0,y1,…,yi能直接从公式算出yi+1的公 式称显式，否则称隐式公式(8.4),(8.6)都是显式，(8.5)是隐式
在计算yi+1时，只需要yi的值，则称公式为单步法；若除yi之外 还需要以前的yi-1等多个值，则称多步法公式.(8.4)、(8.5)是单 计 步法，(8.6)是多步法，确切说是二步法.
0.2 xi yi
计 算
yi
1
yi
0.05
yi
2xi yi
y(0) i 1
2 xi 1 y (0)
i 1
方
y0 1
法 课
件 结果（略）
15
结束
§8.3 龙格-库塔(Runge-Kutta)法
欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出
高精度的显式一步法：龙格-库塔法,简称R-K法. 计 算 8.3.1 二阶R-K法
8.2.1 欧拉方法的导出
把区间［a,b］分为n个小区间，取步长h=(b-a)/ n ，节点 xi=x0+ih,i=0,1,2,…，n，其中x0=a，又设y (x)为上述问题的解.
常微分方程(组)的数值解法
刚性常微分方程组求解
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];
求解初值问题：
2x y' y y y ( 0) 1
(0 x 1)
y ' f ( x, y ) y (a) y 0
ode输入函数 function f=fun(x,y) f=y-2*x/y; 输出变量为因变 量导数的表达式
自变量在前，因变 量在后
2.3 常微分方程（组）的数值解法
知识要点
常微分方程初值问题---ode45,0de23
微分方程在化工模型中的应用
•间歇反应器的计算 •活塞流反应器的计算
•全混流反应器的动态模拟
•定态一维热传导问题
•逆流壁冷式固定床反应器一维模型
•固定床反应器的分散模型
Matlab常微分方程求解问题分类
边值问题：
高阶微分方程odefile的编写
求解： y"a(t )( y' ) 2 b(t ) y et cos 2t
a(t ) e t cos 2te 2t , b(t ) cos(2t )
数值分析-计算方法-常微分方程-精品文档
Ri 的主项Байду номын сангаас
/* leading term */
7.1 Euler’s Method
欧拉公式的改进：
隐式欧拉法 /* implicit Euler method */
向后差商近似导数
y(xi+1)≈yi+hf(xi+1, y(xi+1))
y ( x ) y ( x i 1 i) y ( x ) i 1 xi h
考虑一阶常微分方程的初值问题 /* Initial-Value Problem */:
dy f (x , y) x [a ,b ] dx 0 y(a) y
只要 f (x, y) 在[a, b] R1 上连续，且关于 y 满足 Lipschitz 条 f ( x , y ) f ( x , y ) | L | y y | 1 2 1 2 件，即存在与 x, y 无关的常数 L 使 | 对任意定义在 [a, b] 上的 y1(x) 和 y2(x) 都成立，则上述IVP存 在唯一解。 要计算出解函数 y(x) 在一系列节点 a = x0< x1<…< xn= b y ( x ) ( i 1 , ... , n ) 处的近似值 y i i 节点间距hi=xi+1-xi称为步长，当 hi = h 为常数时称为等步长。
3 h R y ( x ) y y ( x ) O ( h ) i i 1 i 1 i 2
2
即隐式欧拉公式具有 1 阶精度。
7.1 Euler’s Method
梯形公式 /* trapezoid formula */
— 显、隐式两种算法的平均
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
科学计算—理论、方法及其基于MATLAB 的程序实现与分析微分方程（组）数值解法§1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。
如揭示质点运动规律的Newton 第二定律：()()()⎪⎩⎪⎨⎧'='==000022x t x x t x t F dt xd m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：()()00y y t f ay dtdy=+= （２） 的解为：()()()τττd f e y e t y tt a at⎰-+=00 （３）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。
一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：()()000,y t y t t t y t F dtdyf=≤≤= （７）其中()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=t y t y t y t y n 21 （８）()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=y t f y t f y t f y t F n ,,,,21（９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。
§１.1 常微分方程（组）的Cauch 问题数值解法概论假设要求在点（时刻）k t ，n k ,,2,1 =处初值问题（７）的解()00,,y t t y y =的（近似）值,如果已求得k t 时刻的值()k t y 或它的近似值k y (如0t 时刻的值()0t y )，那么将式（７）的两端在区间[]1,+k k t t 上积分()()()()dt t y t f t y t y dt dtdyk kk kt t k k t t ⎰⎰++=-=+11,1 (10)可得()()()()dt t y t f t y t y k kt t k k ⎰++=+1,1 (11)或()()()dt t y t f y y t y k kt t k k k ⎰++=≈++1,11 (12)显然,为了利用式（11）或（12）求得()1+k t y 的精确值（近似值1+k y ），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解, 可以说,不同的式值积分方法将给出不同的Cauch 问题的数值解法。
§1.2 最简单的数值解法——Euler 方法假设要求在点（时刻）kh t t k +=0，nt t h f 0-=，n k ,,2,1 =处初值问题（７）的解()t y y =的近似值。
首先对式（７）的两端积分，得()()()()dt t y t f t y t y dt dtdyk k k kt t k k t t ⎰⎰++=-=+11,1 （13） 对于式（13）的右边，如果用积分下限k t t =处的函数值()()k k t y t f ,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有()()()()()[]k k t tk k t y t hf dt t y t f t y t y k k,,11≈=-⎰++ （14）进而得到下式给出的递推算法—Euler 方法()()001,,2,1,y t y nk y t hf y y k k k k ==+=+ （15）例１ 用Euler 方法解如下初值问题，取3.0=h ，()1030sin 23=≤≤-=y t t y y dtdy解：由（15）得()1010,,2,1sin 6.03.131==-=+y k t y y y kk k k结果如下：如果取1.0=h ，其结果如下图所示：Euler_Method§1.3 改进的Euler 方法对于（15）的右边，如果被积函数用积分限k t t =和1+=k t t 处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有()()()()()()()()[]111,,2,1++++≈=-⎰+k k k k t t k k t y t f t y t f hdt t y t f t y t y k k（16） 进而得到下式给出的递推算法：()()[]()00111,,2,1,,2y t y nk y t f y t f hy y k k k k k k ==++=+++ （17）通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求1+k y 时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了１） 尽可能地保持算法（17）精度高的优点； ２） 尽可能地利用Euler 方法计算简单的长处； 人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 ()k k k k y t hf y y ,01+=+ （18） 修正 ()()[]n k y t f y t f h y y k k k k k k ,,2,1,,2111 =++=+++ （19）()00y t y =例２ Euler 方法与改进的Euler 方法的比较 下图是当3.0=h 时比较的结果：open Improved_Euler_Method.m§1.4 Euler 方法和改进的Euler 方法的误差分析由Taylor 公式()()()()()()()()()()22111,!21h O h t y t f t y t t t y t t t y t y t y k k k k k k k k k k k ++=+-''+-'+=+++（19） 说明Euler 方法的截断误差是()2h O ，类似地，由()()()()()()321!21,h O h t y h t y t f t y t y k k k k k +''++=+ （20） ()()()()()()321111!21,h O h t y h t y t f t y t y k k k k k +''+-=++++ （21） 以及()()()()()()h O t y t y h t y t y t y k k k k k +''=''⇒+'''+''=''++11 (22)让式（20）的两端减式（21）的两端，可得()()()()()()[]()()()()()()()[]()31132111,,241,,2h O t y t f t y t f ht y h O h O h t y t f t y t f ht y t y k k k k k k k k k k k +++=++++=+++++ （23）从上述推导Euler 方法、改进的Euler 方法的过程以及例１、例２容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于： 1 在推导Euler 方法时，我们是用待求解函数()t y y =在一点处的变化 率()()k k t y t f ,代替()t y y =在区间],[1+k k t t 上的平均变化率：()()()()t y t hf dt t y t f k kt t ,,1=⎰+ （24）2 而在推导改进的Euler 方法时，我们是用待求解函数()t y y =在两点处变化率的平均值()()()()[]11,,21+++k k k k t y t f t y t f 代替()t y y =在区间],[1+k k t t 上的平均变化率；显然，通常()()()()[]11,,21+++k k k k t y t f t y t f 比()()k k t y t f ,更接近于()t y y =在区间],[1+k k t t 上的平均变化率。
由此启发人们：适当地选取区间],[1+k k t t 上函数()t y y =若干点处的变化率，用它们加权平均值代替()t y y =在区间],[1+k k t t 上的平均变化率，近似解的精度应更高。
下面将要介绍的Runge —Kutta 法就是基于上述想法得到的。
§2 Runge —Kutta 法Runge —Kutta 法是按选取区间],[1+k k t t 上函数()t y y =变化率()y t f ,的个数的多少和截断误差的阶数()m h O 来区分的一系列方法，如 １ 二阶的Runge —Kutta 法（改进的Euler 方法）()()()⎪⎪⎩⎪⎪⎨⎧++=+==++21111212,,K K hy y hK y t f K y t f K k k k k k k （25） ２ 三阶的Runge —Kutta 法()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=++=⎪⎭⎫ ⎝⎛++==+32112312146,2,2,K K K h y y hK y h t f K K h y h t f K y t f K k k k k k k k k （26） ３ 四阶的Runge —Kutta 法 1) 古典形式()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++++=++=⎪⎭⎫ ⎝⎛++=⎪⎭⎫ ⎝⎛++==+432113423121226,2,22,2,K K K K hy y hK y h t f K K h y h t f K K h y h t f K y t f K kk k k k k k k k k （27） 2) Gill 公式(具有减小舍入误差的优点)()()()()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+++-++=⎪⎪⎭⎫ ⎝⎛++-+=⎪⎪⎭⎫ ⎝⎛-+-++=⎪⎭⎫ ⎝⎛++==+432113242131212222622222,222212,22,2,K K K K hy y hK hK y h t f K hK hK y h t f K K h y h t f K y t f K k k k k k k k k k k （28） 4 Runge —Kutta 法的一般形式()∑=++=+=ni i i k k k k k K a h y h y t h y y 11,,φ (29)其中()h y t k k ,,φ称为增量函数(Increment Function)以及()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++++=+++=++==-----11,111,1122212123111121,,,,n n n n k n k n k k k k k k hK q hK q y h p t f K hK q hK q y h p t f K hK q y h p t f K y t f K （30）需要特别指出的是：在确定（29）、（30）中的参数i a ，i p 和ij q 时，应该使（29）的右端和适当阶的（如n 阶）Taylor 展式一致，这样，至少对于底阶的R-K 法来说，参与加权的斜率个数n 与方法的阶数是一致的。