常微分方程计算

03 拉普拉斯变换的逆变换
定义与性质
定义
逆变换是拉普拉斯变换的逆过程，将 拉普拉斯变换后的函数还原为原函数。
性质
逆变换具有线性性、时移性、微分性、 积分性和相似性等性质，这些性质在 求解常微分方程时具有重要作用。
逆变换的求解方法
表格法
通过查表或计算公式，将拉普拉 斯变换后的函数还原为原函数。 这种方法适用于已知拉普拉斯变 换函数的简单情况。
幂级数法
通过幂级数展开，将拉普拉斯变 换后的函数展开为无穷级数，然 后逐项积分得到原函数。这种方 法适用于较为复杂的拉普拉斯变 换函数。
积分法
通过积分运算，将拉普拉斯变换 后的函数进行积分，得到原函数。 这种方法需要熟练掌握积分运算 和拉普拉斯变换的性质。
04 拉普拉斯变换法的优缺点
优点
高效性
对于一些复杂或难以直接求 解的常微分方程，拉普拉斯 变换法能够提供一种简洁、 高效的求解方法。
普适性
拉普拉斯变换法适用于各种 类型的初值问题，具有广泛 的适用性。
易于计算
拉普拉斯变换的逆变换相对 容易计算，使得求解过程相 对简单。
可处理多变量问题
通过引入偏导数，拉普拉斯 变换法可以处理多变量微分 方程，这是其他方法难以做 到的。
缺点
不易理解物理意义
拉普拉斯变换将原始的微分方程转换为复 平面上的函数，这使得初学者不易理解其
性质
拉普拉斯变换具有线性性、时移性、 微分性、积分性和复共轭性等性质， 这些性质使得求解常微分方程变得更 为简便。
拉普拉斯变换的应用
求解常微分方程
通过拉普拉斯变换，可以将常微分方程转化为代数方程，从而简化求 解过程。
系统分析
在控制工程和信号处理等领域，拉普拉斯变换被广泛应用于系统分析 和系统设计。

计 算 题（每题10分）1、求解微分方程2'22x y xy xe -+=。

2、试用逐次逼近法求方程2y x dxdy+=通过点(0,0)的第三次近似解. 3、求解方程'2x y y y e -''+-=的通解4、求方程组dx dt ydydtx y ==+⎧⎨⎪⎩⎪2的通解5、求解微分方程'24y xy x +=6、试用逐次逼近法求方程2y x dxdy-=通过点(1,0)的第二次近似解。

7、求解方程''+-=-y y y e x '22的通解8、求方程组dxdt x ydydtx y =+=+⎧⎨⎪⎩⎪234的通解9、求解微分方程xy y x '-2=24 10、试用逐次逼近法求方程2y x dxdy-=通过(0,0)的第三次近似解. 11、求解方程''+-=-y y y e x '24的通解12、求方程组dxdtx y dydtx y =+=+⎧⎨⎪⎩⎪2332的通解13、求解微分方程x y y e x (')-=14、试用逐次逼近法求方程22x y dxdy+=通过点(0,0)的第三次逼近解. 15、求解方程''+-=--y y y e x '22的通解16、求解方程x e y y y -=-+''32 的通解17、求方程组⎪⎩⎪⎨⎧-+=-+=yx dt dydtdx x y dt dy dt dx243452的通解 18、解微分方程22(1)(1)0x y dx y x dy -+-= 19、试用逐次逼近法求方程2dyx y dx=-满足初始条件(0)0y =的近似解：0123(),(),(),()x x x x ϕϕϕϕ.20、利用逐次逼近法，求方程22dyy x dx=-适合初值条件(0)1y =的近似解：012(),(),()x x x ϕϕϕ。

# 定义网格密度
arrows=LINE,
# 定义线段类型
axes=NORMAL);
# 定义坐标系类型
在MATLAB的向量场命令为 quiver(x,y,px,py)
回车后Maple就在1 1 的网格点上画出了向量场
44
的图形，并给出了过点(-2, 2) (-2 ,1) (-2,2) 的三
条积分曲线，见下图
M (x,y)co x s2xye , y
N (x,y)co x s2xye x
M(x,y)N(x,y)
y
x
所以方程为全微分方程。
由公式F (x ,y ) 0M (s ,y )d s 0N (0 ,s )d s
x(yc o ss 2 se y)d sy2 d s
0
0
ysinxx2ey2y
或
x
y
F (x ,y )x 0M (s ,y ) d s y 0N (x 0 ,s ) d
s
例：验证方程
( y c o s x 2 x e y ) d x ( s i n x x 2 e y 2 ) d y 0
是全微分方程，并求它的通解。 解：由于 M (x ,y ) y c o sx 2 x e yN (x ,y ) s in x x 2 e y 2
dx
令 zy1n,则 dz(1n)yndy
dx
dx
d z (1 n )P (x )z (1 n )Q (x )
d x
求出此方程通解后, 换回原变量即得伯努利方程的通解。
例 湖泊的污染
设一个化工厂每立方米的废水中含有3.08kg盐酸， 这些废水流入一个湖泊中，废水流入的速率20 立方米每小时. 开始湖中有水400000立方米. 河水 中流入不含盐酸的水是1000立方米每小时, 湖泊 中混合均匀的水的流出的速率是1000立方米每小 时,求该厂排污1年时, 湖泊水中盐酸的含量。 解: 设t时刻湖泊中所含盐酸的数量为 x ( t )

常微分方程的simpson公式1. 概述常微分方程是描述自变量的一阶导数和未知函数之间关系的方程。

在实际问题中，常微分方程的求解十分重要，可以描述种种自然现象的演化规律，也是物理、生物、经济等多个学科的基础。

在数值求解常微分方程时，simpson公式是一种被广泛应用的数值积分方法，能够在离散的点上对函数进行数值积分。

2. Simpson公式的原理和应用Simpson公式是一种数值积分方法，其原理基于积分中值定理和牛顿-卡索公式。

对于一个给定的函数 f(x)，我们希望计算其在区间 [a, b] 上的定积分∫(a到b) f(x)dx。

Simpson公式给出的积分近似值为\[ \int_{a}^{b} f(x)dx \approx \frac{h}{3} [f(x_0) + 4f(x_1) + 2f(x_2) + \cdots + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)] \]其中，h 是区间 [a, b] 的步长，即 h=(b-a)/n，n 为等分数。

Simpson公式的应用非常广泛，尤其在解常微分方程时扮演重要角色。

利用数值积分方法，我们可以将常微分方程的求解转化为对一阶导数的数值积分。

具体来说，若已知常微分方程的初值条件，我们可以通过适当的步长和积分方法得到未知函数在一系列离散点的近似值，从而得到整个函数的近似解。

3. Simpson公式的优缺点作为一种数值积分方法，Simpson公式具有一定的优缺点。

它的计算简单直观，易于理解和编程实现，因此在实际应用中非常方便。

与其他数值积分方法相比，Simpson公式具有较高的精度和稳定性，特别是对于光滑的函数和较小的步长。

但是，Simpson公式也存在一些不足之处，例如在处理非光滑函数、函数不连续的情况下可能存在较大误差，需要做适当的修正和优化。

4. 结论Simpson公式是一种简单而实用的数值积分方法，广泛应用于常微分方程的求解和其他科学计算领域。

称为局部截断误差。 显然，这个误差在逐 步计算过程中会传播， 积累。因此还要估计 这种积累
yn +1 = yn + hf ( xn , yn )
步计算是精确的前提下， ， 定义 在假设 yi = y(xi)，即第 i 步计算是精确的前提下，考虑 称为局部截断误差 的截断误差 Ri = y(xi+1) − yi+1 称为 局部截断误差 /* local
当 β−1≠0 时，为隐式公 则为显式公式 显式公式。 式; β−1=0 则为显式公式。
f n = f ( xn , y n )
基于数值积分的构造法 上积分， 将 y′ = f ( x , y ) 在[ xn− p , xn+1 ] 上积分，得到
xi +1 xn− p
y ( xn +1 ) − y ( xn − p ) = ∫
f ( x, y ( x))dx
x n +1
n− p
只要近似地算出右边的积分 I k ≈ ∫x f ( x, y ( x )) dx ，则可通 只要近似地算出右边的积分 过 yn+1 = yn− p + I k近似y(xn+1) 。而选用不同近似式 Ik，可得到不 同的计算公式。 同的计算公式。
为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论： ① 步长充分小时，所得到的数值解能否逼近问题得真解；即收敛性问题 ② 误差估计 ③ 产生得舍入误差，在以后得各步计算中，是否会无限制扩大；稳定性问题
8.1 Euler公式 公式
做等距分割，利用数值微分代替导数项，建立差分方程。
为了使解存在唯一，一般，要加限制条件在f上，要求f对y满足Lipschitz条件：

18—1 常微分方程数值解法2§1 引言§2 Euler 方法§3 Runge －Kutta 方法§4 单步法的收敛性与稳定性§5 线性多步法§6 方程组与高阶方程的情况§7 边值问题的数值解法3§1 引言微分方程：关于一个未知函数的方程，方程中含有未知函数的（偏）导数，以及自变量等，其中关于未知函数导数的最高次数称为微分方程的阶数.例如：0)()(')()(''=++−x c y x b y x a x y4实际中，很多问题的数学模型都是微分方程. 常微分方程作为微分方程的基本类型之一，在理论研究与工程实际上应用很广泛. 很多问题的数学模型都可以归结为常微分方程. 很多偏微分方程问题，也可以化为常微分方程问题来近似求解.微分方程的应用情况5对于一个常微分方程：'(,) ,[,]dy y f x y x a b dx==∈为了使解存在，一般要对函数f 施加限制条件，例如要求f 对y 满足Lipschitz 条件：1212(,)(,)f x y f x y L y y −≤−6同时，一个有解的微分方程通常会有无穷多个解例如cos() sin(),dyx y x a a R dx=⇒=+∀∈为了使解唯一，需要加入一个限定条件. 通常会在端点出给出，如下面的初值问题：(,),[,]()dyf x y x a b dx y a y ⎧=∈⎪⎨⎪=⎩7常微分方程的解是一个函数，但是，只有极少数特殊的方程才能求解出来，绝大多数是不可解的.并且计算机没有办法对函数进行运算. 一般考虑其近似解法，一种是近似解析法，如逼近法、级数解法等，另一种是本章介绍的数值解法.8§2 Euler 方法92－1 Euler 公式对常微分方程初值问题：⎩⎨⎧==00')(),(y x y y x f y 数值求解的关键在于消除其中的导数项——称为离散化. 利用差商近似逼近微分是离散化的一个基本途径.10现在假设求解节点为),,1,0(m i ih a x i "=+=，其中ma b h −=为步长，这些节点相应的函数值为)(,),(1m x y x y ". 在点n x 处，已知))(,()('n n n x y x f x y =用n x 的向前差商nn n n x x x y x y −−++11)()(近似代替)('n x y ，如§1，则得到所谓的Euler 公式1(,)n n n n y y hf x y +=+——单步、显式格式11Euler 公式的局部截断误差：假设)(n n x y y =情况下，11)(++−n n y x y 称为局部截断误差.'''2311''23()()()()()2()(,()(()))2n n n n n n n n n y x y x y y x hy x h O h y x h y x f x y x h O h ++−=+++−−=+故有)(2)(''211n n n x y h y x y ≈−++. 122－2 后退的Euler 公式同样对常微分方程初值问题，在1+n x 点，已知))(,()(111'+++=n n n x y x f x y ，如果用向后差商hx y x y n n )()(1−+代替)(1'+n x y ，则得到后退的Euler 公式：111(,)n n n n y y hf x y +++=+——单步、隐式格式13相对于以上可以直接计算1+n y 的Euler 公式（显式），上式是隐式公式. 一般来讲，显式容易计算，而隐式具有更好的稳定性.求解上述公式，通常使用迭代法：对于给定的初值)0(1+n y，计算(1)()111(,)(0,1,)k k n n n n y y f x y k ++++=+="， 如果)(1lim k n k y +∞→收敛，则其极限必满足上述后退Euler 公式.14局部截断误差：假设)(n n x y y =，则),()(111++++=n n n n y x hf x y y .由于)]()[,())(,(),(1111111+++++++−+=n n n y n n n n x y y x f x y x f y x f η且''''2111(,())()()()()n n n n n f x y x y x y x hy x O h +++==++15则有'2''31111(,)[()]()()()()n y n n n n n n y hf x y y x y x hy x h y x O h η++++=−++++将此式减去式2'''31()()()()()2n n n n h y x y x hy x y x O h +=+++ 可得，2''311111()(,)[()]()()2n n y n n n n h y x y hf x y x y y x O h η+++++−=−−+16考虑到21111(,)()1(,)y n y n hf x O h hf x ηη++=++−，则有22''3''11()()()()22n n n n h h y x y y x O h y x ++−=−+≈−172－3 梯形公式由于上述两个公式的局部截断误差绝对值相等，符号相反，故求其算术平均得到梯形公式：111[(,)(,)]2n n n n n n hy y f x y f x y +++=++——单步、隐式格式18梯形法同样是隐式公式，可用下列迭代公式求解：(0)1(1)()111(,)[(,)(,)]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 +++++⎧=+⎪⎨=++⎪⎩局部截断误差：类似于后退Euler ，可计算出)(12)('''311n n n x y h y x y −≈−++192－4 改进的Euler 公式上述用迭代法求解梯形公式虽然提高了精度，但计算量也很大. 实际上常采用的方法是，用Euler 公式求得初始值（预测），然后迭代法仅施行一次（校正）——改进的Euler 公式：1111(,)[(,)(,)]2n n n n n n n n n n y y f x y hy y f x y f x y ++++⎧=+⎪⎨=++⎪⎩20估计上式中第二式当1+n y 为准确值时的局部截断误差：''11113(3)()()(()[()()])2()12n n n n n n n hy x y y x y x y x y x hy x ++++−=−++≈−212－5 Euler 两步公式如果用中心差商hx y x y n n 2)()(11−+−代替)('n x y ，则得Euler 两步公式112(,)n n n n y y hf x y +−=+——两步、显式格式22假设1−n y 及n y 均为准确值，利用Taylor 展式容易计算Euler 两步公式的局部截断误差为：11113(3)()()(()2(,()))()3n n n n n n n y x y y x y x hf x y x h y x +++−−=−+≈23此式与梯形公式相结合，得到如下的预测－校正公式：111112(,)[(,)(,)]2n n n n n n n n n n y y hf x y hy y f x y f x y −++++⎧=+⎪⎨=++⎪⎩假设第一式中的1−n y 及n y ，以及第二式中的n y 及1+n y 均是准确值，则有，2441)()(1111−≈−−++++n n n n y x y y x y 从而可得以下的事后估计式，111111114()()51()()5n n n n n n n n y x y y y y x y y y ++++++++⎧−≈−−⎪⎪⎨⎪−≈−⎪⎩25可以期望，以上式估计的误差作为计算结果的补偿，可以提高计算精度.以n p 及n c 分别表示第n 步的预测值和校正值，则有以下的“预测－改进－校正－改进”方案（其中在1+n p 与1+n c 尚未计算出来的前提下，以n n c p −代替11++−n n c p ：26预测：'112n n n hy y p +=−+预测的改进：)(5411n n n n c p p m −−=++计算：),(11'1+++=n n n m x f m校正：)(2'1'1++++=n n n n m y hy c校正的改进：)(511111++++−+=n n n n c p c y计算：),(11'1+++=n n n y x f y27例 用Euler 方法求解初值问题2'[0,0.6](0)1y y xy x y ⎧=−−∈⎨=⎩取0.2h =，要求保留六位小数. 解：Euler 迭代格式为2210.2()0.80.2k k k k k k k k y y y x y y x y +=+−−=−因此2821000(0.2)0.80.20.8y y y x y ≈=−= 22111(0.4)0.80.20.6144y y y x y ≈=−=23222(0.6)0.80.20.461321y y y x y ≈=−=29例 用改进的Euler 方法求解初值问题2'sin 0[0,0.6](0)1y y y x x y ⎧++=∈⎨=⎩取0.2h =，求(0.2),(0.4)y y 的近似值，要求保留六位小数.解：改进的Euler 格式为212211110.2(sin )0.2(sin sin )2k k k k k k k k k k k k k y y y y x y y y y x y y x +++++⎧=+−−⎪⎨=+−−−−⎪⎩30即，222110.820.08sin 0.1(0.80.2sin )sin k k k k k k k k y y y x y y x x ++=−−−则有1(0.2)0.807285y y ≈=，2(0.4)0.636650y y ≈=31§3 Runge －Kutta 方法Def.1如果一种方法的局部截断误差为)(1+p h O ，则称该方法具有p 阶精度. 323－2 Runge —Kutta 方法的基本思想上述的Taylor 级数法虽然可得到较高精度的近似公式，但计算导数比较麻烦. 这里介绍不用计算导数的方法.))(,()()()('1h x y h x f h x y hx y x y n n n n n θθθ++=+=−+——平均斜率.33如果粗略地以),(n n y x f 作为平均斜率，则得Euler 公式；如果以221K K +作为平均斜率，其中),(1n n y x f K =，),(112hK y x f K n n +=+，则得改进的Euler 公式.343－3 二阶的Runge －Kutta 方法对点n x 和)10(≤<+=+p ph x x n p n ，用这两点斜率的线性组合近似代替平均斜率，则得计算公式：11122121()(,)(,)n n n n n p n y y h K K K f x y K f x y phK λλ++⎧=++⎪=⎨⎪=+⎩35现确定系数p ,,21λλ，使得公式具有二阶精度. 因为，取n y 为()n y x ，则'1(,)(,())'()n n n n n nK f x y f x y x y x y === 再把2K 在),(n n y x 处展开，有36'21(,)(,)n p n n n n K f x y phK f x ph y phy +=+=++代入可得，'2''31122()()n n n n y y hy ph y O h λλλ+=++++'2(,)(,)(,)()n n x n n y n n n f x y f x y ph f x y phy O h =+⋅+⋅+'2(')(,)()n x y n n y ph f f y x y O h =+⋅+⋅+'''2()n n y ph y O h =+⋅+37相比较二阶Taylor 展开''2'12n n n n y h hy y y ++=+，有，⎪⎩⎪⎨⎧==+211221p λλλ满足此条件的公式称为二阶Runge －Kutta 公式.38可以验证改进的Euler 公式属于二阶Runge －Kutta 公式. 下列变形的Euler 公式也是二阶Runge －Kutta 公式：12121(,)(,)22n n n n n n y y hK K f x y h h K f x y K +⎧⎪=+⎪=⎨⎪⎪=++⎩393－4 三阶Runge －Kutta 公式同二阶Runge －Kutta 公式，考虑三点,,(01)n n p n q x x x p q ++≤≤≤试图用它们的斜率321,,K K K 的线性组合近似代替平均斜率，即有如下形式的公式：1112233121312()(,)(,)(,())n n n n n n n n y y h K K K K f x y K f x ph y phK K f x qh y qh rK sK λλλ+=+++⎧⎪=⎪⎨=++⎪⎪=+++⎩40把32,K K 在),(n n y x 处展开，通过与)(1+n x y 在n x 的直接Taylor 展式比较，可确定系数s r q p ,,,,,,321λλλ，满足下式，从而使得上述公式具有三阶精度，41特别地，2,1,1,21,32,61231=−======s r q p λλλ是其一特例.123232223311213161p q p q pqs r s λλλλλλλλ++=⎧⎪⎪+=⎪⎪⎪+=⎨⎪⎪=⎪⎪+=⎪⎩423－5 四阶Runge －Kutta 公式相同的方法，可以导出下列经典的四阶Runge －Kutta 公式：112341213243(22)6(,)(,)22(,)22(,)n n n n n n n n n n h y y K K K K K f x y h h K f x y K h h K f x y K K f x h y hK +⎧=++++⎪⎪=⎪⎪⎪=++⎨⎪⎪=++⎪⎪=++⎪⎩43例 用经典四阶Runge —Kutta 方法求解初值问题'83[0,0.4](0)1y y x y =−⎧∈⎨=⎩，取0.2h =，求(0.4)y 的近似值，要求保留六位小数.解：四阶Runge —Kutta 格式为44112341211123122241330.2(22)6(,)830.2(,)83(0.1) 5.6 2.120.2(,)83(0.1) 6.32 2.372(,0.2)83(0.2) 4.208 1.578k k k k k k k k k k k kk k k k ky y K K K K K f x y y K f x y K y K yK f x y K y K y K f x y K y K y ++++⎧=++++⎪⎪==−⎪⎪⎪=+=−+=−⎨⎪⎪=+=−+=−⎪⎪⎪=+=−+=−⎩则10.5494 1.2016k k y y +=+，45故12(0.2) 2.3004,(0.4) 2.4654y y y y ≈=≈=.注：由准确解382()33xy x e −=−可得(0.2) 2.300792,(0.4) 2.465871y y ==46§5 线性多步法基本思想：在计算1+i y 之前，已计算出一系列的近似值i y y ,,1"，如果充分利用这些已知信息，可以期望会获得更高精度的)(1+i x y 的近似值1+i y .基本方法：基于数值积分与基于Taylor 展开的构造方法.475－1 基于数值积分的构造方法对方程),('y x f y =两边从i x 到1+i x 积分，则得∫++=+1),()()(1i ix x i i dxy x f x y x y 设)(x P r 是f (x , y )的插值多项式，由此可得以下的一般形式的计算公式：∫++=+1)(1i ix x r i i dxx P y y 48例 取线性插值))(,())(,()(11111+++++−−+−−=i i i i ii i i i i r x y x f x x x x x y x f x x x x x P ，则得到梯形法：)],(),([2111+++++=i i i i i i y x f y x f hy y495－2 Adams 显式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(11r i r i i i i i f x f x f x −−−−"构造插值多项式)(x P r ，由牛顿后插公式（注意到：j i j i j f f −Δ=∇）j i jrj j i r f j t th x P −=Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑0)1()(其中!)1()1(j j s s s j s +−−=⎟⎟⎠⎞⎜⎜⎝⎛". 50可得10rj i i rj i jj y y h f αΔ+−==+∑——Adams 显式公式其中1(1)j j t dt j α−⎛⎞=−⎜⎟⎝⎠∫，它可写成：∑=−++=rj ji rj i i f h y y 01β515－3 Adams 隐式公式在区间],[1+i i x x 上利用r ＋1个数据点),(,),,(),,(1111+−+−++r i r i i i i i f x f x f x "构造插值多项式)(x P r ，由牛顿后插公式101)1()(+−=+Δ⎟⎟⎠⎞⎜⎜⎝⎛−−=+∑j i jrj ji r f j t th x P 可得*11rj i i rj i j j y y h f α+−+==+Δ∑——Adams 隐式公式52其中01(1)jj t dt j −−⎛⎞α=−⎜⎟⎝⎠∫，它又可写成： *11ri i rj i j j y y h f β+−+==+∑535－4 Adams 预测－校正公式以r ＝3时的Adams 显式与隐式公式为例. 此时，显式公式为)9375955(243211−−−+−+−+=i i i i i i f f f f hy y 利用Taylor 展式，容易计算局部截断误差为)(720251)5(5i x y h . 54)5199(242111−−+++−++=i i i i i i f f f f hy y 同样利用Taylor 展开可得，其局部截断误差为5(5)19()720i h y x −. 隐式公式为55⎪⎩⎪⎨⎧+−++=−+−+=−−+++−−−+)519),(9(24)9375955(24211113211i i i i i i i i i i i i i f f f y x f hy y f f f f h y y 注 利用2－5节的相同作法同样可以构造更精确的计算过程.可构造利用显式预测，隐式校正的计算公式：56§6 方程组与高阶方程的情形6－1 一阶方程组常微分方程初值问题为⎩⎨⎧==00)(),('y x y y x f y 此时T m y y y ),,(1"=，Tm f f f ),,(1"=. 此时上述的一切方法均可使用，只是注意y 与f 此时为向量.576－2 化高阶方程为一阶方程组解下列的m 阶方程()(1)'(1)(1)000000(,,',,)(),'(),,()m m m m y f x y y y y x y y x y yx y −−−⎧=⎨===⎩""令)1(21,,',−===m m y y y y y y "，则有58'12'23'1'12(,,,,)m m m m y y y y y yy f x y y y −⎧=⎪=⎪⎪⎨⎪=⎪⎪=⎩#"初始条件为：)1(00'002001)(,,)(,)(−===m m y x y y x y y x y "。

求常微分方程的数值解一、背景介绍常微分方程（Ordinary Differential Equation，ODE）是描述自然界中变化的数学模型。

常微分方程的解析解往往难以求得，因此需要寻找数值解来近似地描述其行为。

求解常微分方程的数值方法主要有欧拉法、改进欧拉法、龙格-库塔法等。

二、数值方法1. 欧拉法欧拉法是最简单的求解常微分方程的数值方法之一。

它基于导数的定义，将微分方程转化为差分方程，通过迭代计算得到近似解。

欧拉法的公式如下：$$y_{n+1}=y_n+f(t_n,y_n)\Delta t$$其中，$y_n$表示第$n$个时间步长处的函数值，$f(t_n,y_n)$表示在$(t_n,y_n)$处的导数，$\Delta t$表示时间步长。

欧拉法具有易于实现和理解的优点，但精度较低。

2. 改进欧拉法（Heun方法）改进欧拉法又称Heun方法或两步龙格-库塔方法，是对欧拉法进行了精度上提升后得到的一种方法。

它利用两个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

改进欧拉法的公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\Delta t,y_n+k_1\Delta t)$$$$y_{n+1}=y_n+\frac{1}{2}(k_1+k_2)\Delta t$$改进欧拉法比欧拉法精度更高，但计算量也更大。

3. 龙格-库塔法（RK4方法）龙格-库塔法是求解常微分方程中最常用的数值方法之一。

它通过计算多个斜率来近似函数值，并通过加权平均来计算下一个时间步长处的函数值。

RK4方法是龙格-库塔法中最常用的一种方法，其公式如下：$$k_1=f(t_n,y_n)$$$$k_2=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_1\Delta t}{2})$$ $$k_3=f(t_n+\frac{\Delta t}{2},y_n+\frac{k_2\Delta t}{2})$$ $$k_4=f(t_n+\Delta t,y_n+k_3\Delta t)$$$$y_{n+1}=y_n+\frac{1}{6}(k_1+2k_2+2k_3+k_4)\Delta t$$三、数值求解步骤对于给定的常微分方程，可以通过以下步骤求解其数值解：1. 确定初值条件：确定$t=0$时刻的函数值$y(0)$。

常微分方程的数值解算法常微分方程的数值解算法是一种对常微分方程进行数值计算的方法，这可以帮助我们更好地理解和研究自然现象和工程问题。

在本文中，我们将介绍一些常用的数值解算法，探讨它们的优缺点和适用范围。

常微分方程（ODE）是描述自然现象和工程问题的重要数学工具。

然而，对于许多ODE解析解是无法求出的，因此我们需要通过数值方法对其进行求解。

常微分方程可以写作：y' = f(t, y)其中，y是函数，f是给定的函数，表示y随t的变化率。

这个方程可以写成初始值问题（IVP）的形式：y'(t) = f(t,y(t))，y(t0) = y0其中，y（t0）=y0是方程的初始条件。

解决IVP问题的典型方法是数值方法。

欧拉方法欧拉方法是最简单的一阶数值方法。

在欧拉方法中，我们从初始条件开始，并在t = t0到t = tn的时间内，用以下公式逐步递推求解：y n+1 = y n + hf (t n, y n)其中，f（t n，y n）是点（t n，y n）处的导数， h = tn - tn-1是时间间隔。

欧拉方法的优点是简单易懂，容易实现。

然而，它的缺点是在整个时间段上的精度不一致。

程度取决于使用的时间间隔。

改进的欧拉方法如果我们使用欧拉方法中每个时间段的中间点而不是起始点来估计下一个时间点，精度就会有所提高。

这个方法叫做改进的欧拉方法（或Heun方法）。

公式为：y n+1 = y n + h½[f(t n, y n)+f(tn+1, yn + h f (tn, yn))]这是一个二阶方法，精度比欧拉方法高，但计算量也大一些。

对于易受噪声干扰的问题，改进的欧拉方法是个很好的选择。

Runge-Kutta方法Runge-Kutta方法是ODE计算的最常用的二阶和高阶数值方法之一。

这个方法对定义域内的每个点都计算一个导数。

显式四阶Runge-Kutta方法（RK4）是最常用的Runge-Kutta方法之一，并已得到大量实践的验证。

故方程的通解为 y Y y1 c1ex c2e-2x e-x
8、解：由方程解出 y，得 y x 2 1 p x ， 代入 dx 1 dy 得 dx dp 即 p cx
2x 2p
p
xp
故通解为 y c (x2 1) 1
2
2c
9、解：方程化为 y 2 y 2x3 x
对应的齐次方程 y 2 y 0 的通解为 y=cx2 (4¹) x
6、试用逐次逼近法求方程 dy x y 2 通过点(1,0)的第二次近似解。 dx
7、求解方程 y y 2y 2ex 的通解
8、求方程组
dx dt dy dt
2x y 3x 4y
的通解
9、求解微分方程 xy 2y 2x4
10、试用逐次逼近法求方程 dy x y 2 通过(0,0)的第三次近似解. dx
75、利用代换 y u 将方程 ycos x 2ysin x 3y cos x ex 化简，并求出原方程的 cos x
通解。
76、求下列线性微分方程组
dx
dt dy
2x 2x
4y 2y
4e2t
(1) (2)
dt
dy1
dx
2 y1
2 y2
77、解下列微分方程组
dy2 dx
y2
58、 x x sin at, a 0
59、 2y 5y cos2 x
60、 y 4y xsin 2x
61、 y 2y 3 4sin 2x
62、 y 2y 2y 4ex cos x
63、 y 9y 18cos3x 30sin 3x
64、 x x sin t cos 2t
11、求解方程 y y 2y 4ex 的通解

01
扰动值满足原来的差分方程，如果原差分方程的解是不增长的，即有
03
从而需要
02
这时就能保证Euler方法的稳定性。
04
Euler格式条件稳定
隐式Euler格式是恒稳定(无条件稳定)的
隐式Euler方法
由于λ<0，从而有 与 恒成立。
1
则：
2
而
3
显然：
4
校正后的误差
从而有：
事后估计式
令pn和cn分别代表第n步的预报值和校正值， 和
可作为pn+1和cn+1的改进值。在cn+1未确定前，可用pn-cn来代替pn+1-cn+1进行计算。
改进后的公式
Exercises 习题3的第13题。
设xn-x0=nh≤T(T为常数)，则
从而
显然，如果初值准确，则有h→0,en → 0.
1
Euler格式收敛。
2
04
03
01
02
稳定性
每一步的计算并不严格准确，存在计算误差的传播问题——扰动。
若
则称为稳定的。
Euler格式和隐式Euler格式
稳定性问题的讨论
Euler格式 设在节点值yn上有一扰动值εn，它的传播使节点值yn+1上产生大小为εn+1的扰动值。假设Euler方法的计算过程不再引入新的误差，则扰动值满足：
改进的思路：
01
先用欧拉方法求得一个初步的近似值，记为 (预报值)，代替右侧的yn+1直接计算，得到校正值yn+1。
02
改进的Euler公式
03
或如下平均化形式
例题
精度分析

-179-第十五章 常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。

如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如22x y dxdy+=，于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。

§1 常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是⎪⎩⎪⎨⎧=≤≤=0)(),(y a y bx a y x f dx dy(1)在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条件，即存在常数L ，使得|||),(),(|y y L y x f y x f −≤−这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。

所谓数值解法，就是求问题（1）的解)(x y 在若干点b x x x x a N =<<<<=L 210处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解，n n n x x h −=+1称为由n x 到1+n x 的步长。

今后如无特别说明，我们总取步长为常量h 。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：（i ）用差商近似导数若用向前差商hx y x y n n )()(1−+代替)('n x y 代入（1）中的微分方程，则得),1,0())(,()()(1L =≈−+n x y x f hx y x y n n n n 化简得))(,()()(1n n n n x y x hf x y x y +≈+如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ，则有),1,0(),(1L =+=+n y x hf y y n n n n （2）这样，问题（1）的近似解可通过求解下述问题⎩⎨⎧==+=+)(),1,0(),(01a y y n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出L ,,21y y 。

1
拉普拉斯变换法用于求解常微分方程的基本思路：
对常微分方程进行拉氏变换法，得代数方程，求解 再反变换获取原方程的解
问题： 1. 什么是拉氏变换 2. 拉氏变换的基本性质 3. 什么是拉氏逆变换 4. 如何用拉氏变换求解微分方程
2
1拉普拉斯变换定义(简称拉氏变换)
对于在 [0, ) 上有定义的函数 f (t)
sx0(n2)
x (n1) 0
16
x(n) a1x(n1) an1x an x f (t)
给(4.32)两端施行Laplace Transform
sn
X
(s)
s n1 x0
sn2 x0
sx0(n2)
x (n1) 0
a1[s n1 X
(s)
sn2 x0
s n3 x0
x (n2) 0
]
an1[sX (s) x0 ] an X (s) F (s)
F (s) test f (t)dt
0
F (n) (s) (1)n t nest f (t)dt
0
F (n) (s) (1)n L[tn f (t)]
10
§3 拉普拉斯逆变换 已知象函数，求原函数
L1[F (s)] f (t)
也具有线性性质
L1[c1F1(s) c2F2 (s)] c1L1[F1(s)] c2L1[F2 (s)]
(sn a1sn1 an1s an ) X (s) F (s) B(s)
X (s) F(s) B(s) A(s)
x(t) L1[ X (s)] L1[ F (s) B(s)] A(s)
17
用拉氏变换求微分方程实例
例5 求 dx x e2t 满足初始条件 x(0) 0的特解

12
7.1.1 欧拉法及其截断误差
4、欧拉公式的截断误差是O(h2)，公式是1 阶的。
因为
yi+ 1 ? yi
1
h f ( x i , y i ) = y ( x i ) + h y ¢( x i )
1 2
1
2n ) ( n y ( x ) y (( x)i 1 ) ( y i()x i ) y y ( ) ( x h i ) ( )yh ( ( x x i ) xi y x ) ( xi ) x y y 2 n! 2
y0 y( x0 )
i1 i i i
(欧拉公式)
9
7.1.1 欧拉法及其截断误差
例 取步长 h=0.1，用欧拉法求解初值问题
ì y ¢= x + y ï ï í ï y (0) = 1 ï î
y i 1 y i h f ( x i , y i ) , i 0 ,1 , 2 , y0 y( x0 )
y f ( x , y ), y( x0 ) y0
x [a , b ]
23
7.1 欧拉法和改进的欧拉法
欧拉公式
y i 1 y i h f ( x i , y i ) , i 0 ,1 , 2 , y0 y( x0 )
h y i 1 y i [ f ( x i , y i ) f ( x i 1 , y i 1 )] , i 0 ,1 , 2 , 2 y y( x ) 0 0
( p)
1.2 ? 1.24
1.528
y 2 = y 1 + 0 .1[( x1 + y 1 ) + ( x 2 + y 2 )] = 1 .2 4 + 0 .1(0 .2 + 1 .2 4 + 0 .4 + 1 .5 2 8) = 1 .5 7 6 8

-1-第十五章 常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。

如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如22dyy x dx=+，于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。

§1 常微分方程的离散化下面主要讨论一阶常微分方程的初值问题，其一般形式是(,)()dyf x y a x bdxy a y ⎧=≤≤⎪⎨⎪=⎩ (1) 在下面的讨论中，我们总假定函数(,)f x y 连续，且关于y 满足李普希兹(Lipschitz)条件，即存在常数L ，使得|(,)(,)|||f x y f x y L y y -≤-这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。

所谓数值解法，就是求问题（1）的解 y (x )在若干点012N a x x x x b =<<<<=处的近似值(1,2,,)n y n N =的方法，(1,2,,)n y n N = 称为问题（1）的数值解，1n n n h x x +-=称为由n x 到1n x +的步长。

今后如无特别说明，我们总取步长为常量h 。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商()()1n n y x y x h+-代替()n y x '代入（1）中的微分方程，则得()()()()1,(0,1,)n n n n y x y x f x y x n h+-≈=化简得()()()()1,n n n n y x y x hf x y x +≈+如果用()n y x 的近似值n y 代入上式右端，所得结果作为()1n y x +的近似值，记为1n y +， 则有()1,(0,1,)n n n n y y hf x y n +=+=(2)这样，问题（1）的近似解可通过求解下述问题()10,(0,1,)()n n n n y y hf x y n y y a +⎧=+=⎨=⎩ (3) 得到，按式（3）由初值0y 可逐次算出1y ,2y ,…。

常微分方程的解法常微分方程（Ordinary Differential Equations, ODE）是数学中的一个重要分支，它研究的是包含未知函数及其导数的方程。

在科学和工程领域中，常微分方程被广泛应用于描述自然现象和系统行为的数学模型。

解常微分方程是研究ODE的核心问题，本文将介绍几种常见的常微分方程解法。

一、分离变量法对于某些可分离变量的常微分方程，我们可以通过将未知函数和变量分离来求解方程。

具体步骤如下：1. 将方程变形，将所有含有未知函数及其导数的项移到等式的一侧；2. 将含有未知函数的项移到一侧，含有变量的项移到另一侧；3. 对两边同时积分，得到解的形式。

例如，考虑求解以下常微分方程：$$\frac{{dy}}{{dx}} = x^2$$将方程分离变量并进行积分，得到：$$\int{1}\ dy = \int{x^2}\ dx$$积分后得到：$$y = \frac{{x^3}}{{3}} + C$$其中C为积分常数，代表无穷多个可能的解。

二、线性线性常微分方程是指方程中的未知函数及其导数项构成一个线性组合的方程。

对于形如$${{d^n y}\over{dx^n}} + a_{n-1}{{d^{n-1} y}\over{dx^{n-1}}} + \ldots + a_1{{dy}\over{dx}} + a_0y = f(x)$$的线性常微分方程，其中$f(x)$为已知函数，我们可以使用特征方程来求解。

1. 求解特征方程$${{d^n r}\over{dr^n}} + a_{n-1}{{d^{n-1}r}\over{dr^{n-1}}} + \ldots + a_1{{dr}\over{dr}} + a_0r = 0.$$特征方程的解为$r_1, r_2, \ldots, r_n$；2. 如果特征方程的解都是实数，则对应的齐次解为$$y_c(x) =C_1e^{r_1x} + C_2e^{r_2x} + \ldots + C_ne^{r_nx}$$其中$C_1, C_2,\ldots, C_n$为常数；3. 如果特征方程的解包含复数，则对应的齐次解为$$y_c(x) =e^{\alpha x}(C_1 \cos(\beta x) + C_2 \sin(\beta x))$$其中$\alpha$和$\beta$是复数，$C_1$和$C_2$是常数；4. 采用常数变易法，设待求的解可以表示为$$y_p(x) =u_1(x)e^{r_1x} + u_2(x)e^{r_2x} + \ldots + u_n(x)e^{r_nx}$$将$u_1(x),u_2(x), \ldots, u_n(x)$代入原方程得到未知常数的方程组，并解此方程组得到$u_1(x), u_2(x), \ldots, u_n(x)$；5. 根据待定系数法，将所有齐次解$y_c(x)$和特解$y_p(x)$相加，得到原方程的通解$y(x) = y_c(x) + y_p(x)$。

常微分方程的解法及应用常微分方程是数学中的一个重要分支，广泛应用于各个领域，例如物理学、生物学、经济学等。

本文将介绍常微分方程的解法和应用。

一、常微分方程的解法常微分方程是描述物理现象和自然现象的重要数学工具，例如天文学、电子学、量子力学、流体力学、热力学、生物学、化学等。

常微分方程主要分为初值问题和边值问题两种。

1.初值问题初值问题是指在某个初始时刻$t_0$，系统的状态已知，求在此后的任意时间$t$内该系统的状态。

其一般形式如下：$$\frac{dy}{dt}=f(y,t), \ \ \ \ y(t_0)=y_0$$其中，$y$是未知的函数，$f$是已知的函数，$y_0$是已知的常数。

2.边值问题边值问题是指在某个区间$[a,b]$内，系统的状态已知，求满足某个条件的函数$y(t)$。

其一般形式如下：$$\frac{d^2y}{dt^2}=f(y,t), \ \ \ \ y(a)=y_A, \ \ \ \ y(b)=y_B$$其中，$y_A$和$y_B$是已知的常数。

3.解法常微分方程的解法有多种方法，下面介绍比较常用的两种方法：欧拉法和四阶龙格-库塔法。

（1）欧拉法欧拉法是常微分方程求解的一种最简单的数值方法，它的基本思想是将微分方程转化为差分方程，利用差分方程求解。

假设在时间t时，y的值为$y(t)$，而在时间$t+h$时的y的值可以用下式计算：$$y(t+h)=y(t)+h\times f(y(t),t)$$其中，$f(y,t)$是微分方程的右端函数，$h$是每次迭代的步长。

（2）四阶龙格-库塔法四阶龙格-库塔法是常微分方程求解的一种较为精确的数值方法，其基本思想是采用区间加权平均法对微分方程进行求解。

四阶龙格-库塔法是由四个步骤组成，分别为：1）计算斜率$k_1=f(y_i,t_i)$2）计算斜率$k_2=f(y_i+\frac{h}{2}k_1,t_i+\frac{h}{2})$3）计算斜率$k_3=f(y_i+\frac{h}{2}k_2,t_i+\frac{h}{2})$4）计算斜率$k_4=f(y_i+hk_3,t_i+h)$将这四个斜率加权平均后即得到四阶龙格-库塔法的解式：$$y_{i+1}=y_i+\frac{1}{6}(k_1+2k_2+2k_3+k_4)$$二、常微分方程的应用常微分方程广泛应用于各个领域，本节将介绍三个常微分方程的应用：自然增长模型、振动模型和物理模型。

已知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