第十章 常微分方程(组)的求解-2014120120190541
第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)为上述问题的解.
常微分方程的解法总结总结
常微分方程的解法总结前言常微分方程（Ordinary Differential Equation，ODE）是研究一阶或高阶导数与未知函数之间关系的数学方程。
在物理学、工程学和计算机科学等领域，常微分方程扮演着重要的角色。
解决常微分方程是这些领域中许多问题的关键。
本文将总结常用的常微分方程解法方法，帮助读者加深对常微分方程的理解并提供解决问题的思路。
一、可分离变量法可分离变量法是一种常见且简单的求解常微分方程的方法。
它适用于形如dy/dx = f(x)g(y)的一阶常微分方程。
解题思路：1.将方程写成dy/g(y) = f(x)dx的形式，将变量进行分离。
2.两边同时积分得到∫(1/g(y))dy = ∫f(x)dx。
3.求出积分后的表达式，并整理得到解 y 的表达式。
使用这种方法解决常微分方程的步骤相对简单，但要注意确认分母不为零以及选取合适的积分常数。
二、特殊方程类型的求解除了可分离变量法，常微分方程还存在一些特殊的方程类型，它们可以通过特定的方法进行解决。
1. 齐次方程齐次方程是指形如dy/dx = F(y/x)的方程。
其中，F(t) 是一个只有一个变量的函数。
解题思路：1.令 v = y/x，即 y = vx。
将方程转化为dy/dx = F(v)。
2.对于dv/dx = F(v)/x这个方程，可以使用分离变量法进行求解。
3.求出 v(x) 后，将其代入 y = vx 得到完整的解。
2. 齐次线性方程齐次线性方程是指形如dy/dx + P(x)y = 0的方程。
解题思路：1.使用积分因子法求解，将方程乘以一个积分因子，使得左边变成一个可积的形式。
2.求积分因子的方法是根据公式μ = e^(∫P(x)dx)，其中 P(x) 是已知的函数。
3.通过乘积的方式求解完整的方程。
3. 一阶线性常微分方程一阶线性常微分方程是指形如dy/dx + P(x)y = Q(x)的方程。
解题思路：1.使用积分因子法，将方程乘以一个积分因子，使得左边变成一个可积的形式。
如何求解常微分方程
如何求解常微分方程？常数变易法、积分因子法，函数变换法。
大致与微积分同时产生。
事实上，求y′＝f(x)的原函数问题便是最简单的微分方程。
I.牛顿本人已经解决了二体问题：在太阳引力作用下，一个单一的行星的运动。
他把两个物体都理想化为质点，得到3个未知函数的3个二阶方程组，经简单计算证明，可化为平面问题，即两个未知函数的两个二阶微分方程组。
用现在叫做“首次积分”的办法，完全解决了它的求解问题。
17世纪就提出了弹性问题，这类问题导致悬链线方程、振动弦的方程等等。
总之，力学、天文学、几何学等领域的许多问题都导致微分方程。
在当代，甚至许多社会科学的问题亦导致微分方程，如人口发展模型、交通流模型……。
因而微分方程的研究是与人类社会密切相关的。
当初，数学家们把精力集中放在求微分方程的通解上，后来证明这一般不可能，于是逐步放弃了这一奢望，而转向定解问题：初值问题、边值问题、混合问题等。
但是，即便是一阶常微分方程，初等解(化为积分形式)也被证明不可能，于是转向定量方法（数值计算）、定性方法，而这首先要解决解的存在性、唯一性等理论上的问题。
方程对于学过中学数学的人来说是比较熟悉的；在初等数学中就有各种各样的方程，比如线性方程、二次方程、高次方程、指数方程、对数方程、三角方程和方程组等等。
这些方程都是要把研究的问题中的已知数和未知数之间的关系找出来，列出包含一个未知数或几个未知数的一个或者多个方程式，然后取求方程的解。
但是在实际工作中，常常出现一些特点和以上方程完全不同的问题。
比如：物质在一定条件下的运动变化，要寻求它的运动、变化的规律；某个物体在重力作用下自由下落，要寻求下落距离随时间变化的规律；火箭在发动机推动下在空间飞行，要寻求它飞行的轨道，等等。
物质运动和它的变化规律在数学上是用函数关系来描述的，因此，这类问题就是要去寻求满足某些条件的一个或者几个未知函数。
也就是说，凡是这类问题都不是简单地去求一个或者几个固定不变的数值，而是要求一个或者几个未知的函数。
常微分方程的解法及应用(常见解法及举实例)
常微分方程的解法及应用 （常见解法及举实例）
课 程 名 称： 高等数学（2） 专 业 班 级： 成 员 组 成：
联 系 方 式：
2012 05月25 年日
摘要
常微分方程是微积分学的重要组成部分，广泛用于具体问题的 研究中。求解常微分的问题，常常通过变量分离、两边积分， 如果是高阶的则通过适当的变量代换，达到降阶的目的来解决 问题。本文就是对不同类型的常微分方程的解法的系统总结： 先对常微分方程
解法：
若得其解为则 原方程通解为
2.4二阶线性微分方程解的结构
形如： 若时，（方程一）称为：二阶线性齐次微分方程。
若时，（方程二）称为：二阶非齐次微分方程
2.4.1 二阶线性齐次微分方程解的结构
定理1 :如果函数与是方程(5.2)的两个解, 则
也是（方程一)的解，其中是任意常数.
定理2 : 如果与是方程(5.2)的两个线性无关的特解，则
2.1.4 伯努利方程
形如：
当时， 一阶线性微分方程（公式法） 当时， 可分离变量微分方程 求通解过程： 作变量代换
（积分因子公式法）
2.2 一阶微分方程的应用举例
例1细菌的增长率与总数成正比。如果培养的细菌总数在24h内 由100增长为400、那么前12h后总数是多少? 分析：
例2。。某人的食量是2500 cal／天，其中1200 cal用于基本的 新陈代谢(即自动消耗)。在健身训练中，他所消耗的大约是16 cal/kg/天，乘以他的体重(kg)。假设以脂肪形式贮藏的热量 100%的有效，而1kg脂肪含热量10,000 cal。求出这人的体重是 怎样随时间变化的。 输入率=2500 cal／天
定义及一般解法做简单阐述，然后应用变量替换法解齐次性微 分方程，降阶法求高阶微分方程，讨论特殊的二阶微分方程， 并且用具体的实例分析常微分方程的应用。
常微分方程的求解
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 "。
第十章.常微分方程数值解法
n+1 n k +1 n+1 n
(
n
n
k
n+1
n+1
)
k = 0, ,2,L 1
（二）向后Euler 法的稳定性 对y = λy 用 后 向 Euler法： = y + hλ y y 误 公 ： 差 式 ρ = ρ + hλ ρ ρ = 1 = 1 则 1− hλ ρ [(1−λh)(1−λh)]
( 0)
xn
(
(
))
(
)
(
)
梯形公式的稳定性
梯形公式比Euler 方法的局部与总体误差均高一阶， = 0(h )， e
2 N
但每次迭代均多算一次函数值—提高精度的计算代价
ρ
n+ 1
= ρn +
λh
2
(ρ + ρ )
n n+ 1 1 2
计算方法常微分方程求解
x1
y ( x1 ) y0 h f ( x1 , y ( x1 ))
y i 1 y i h f ( x i 1 , y i 1 ) ( i 0, ... , n 1)
由于未知数 yi+1 同时出现在等式的两边，不能直接得到，故 称为隐式 欧拉公式，而前者称为显式 欧拉公式。
1 1 2 1 , 2 p 2
存在无穷多个解。所有满足上式的格式统称为2阶龙格 - 库 塔格式。注意到，p 1, 1 2 1 就是改进的欧拉法。
2
为获得更高的精度，进一步推广
yi1 yi h[ 1 K1 2 K2 ... m Km] K1 f ( xi , yi ) K2 f ( xi 2 h, yi 21 hK1 ) K3 f ( xi 3 h, yi 31 hK1 32 hK2 )
yi (1 h)i y0
对任意固定的 x = xi = i h ，有 yi y0 (1 h) xi / h y0[(1 h)1/ h ]xi
y0e xi y( xi )
稳定性
定义 若某算法在计算过程中任一步产生的误差在以后的计
算中都不会增长，则称该算法是稳定的 。
... ... Km f ( xi m h, y m1 hK1 m2 hK2 ... m m1 hKm1 )
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定 系数，确定这些系数的 步骤与前面相似。
欧拉公式：
y( x1 ) y( x0 ) y( x)dx y( x0 ) hy( x0 ) y0 h f ( x0 , y0 )
常微分方程的常见解法
实例解析
实例1
求解一阶线性常微分方程 $y' + p(x)y = q(x)$，通过引入参数 $lambda$，可以将方程转化为 $lambda y = q(x)$，从而简化求解过程。
实例2
求解二阶常微分方程 $y'' + y' + y = 0$，通过引入参数 $lambda$，可以将方程转化为 $lambda^2 + lambda + 1 = 0$，从而求解出 $lambda$ 的值，进一步得到原方程的解。
当 (M(x)) 和 (N(x)) 均为非零函数时，该方法适用。
实例解析
1. 确定积分因子
选择积分因子为 (e^x)
5. 解出原方程
将 (e^x y = frac{1}{3} e^{3x} + C) 代入 原方程，解得 (y = frac{1}{3} x^2 + Ce^{-x})
4. 解方程
对两边积分，得到 (e^x y = frac{1}{3} e^{3x} + C)
04 积分因子法
定义与特点
定义
积分因子法是一种通过引入一个因子来简化微分方程的方法。
特点
通过乘以一个适当的因子，可以将微分方程转化为可分离变量的形式，从而简化求解过程。
适用范围
适用于形如 (M(x)y' + N(x)y = f(x)) 的线性微分方程，其中 (M(x)) 和 (N(x)) 是 已知函数，(f(x)) 是给定的函数。
实例2
考虑一阶常微分方程 (dy/dx = xy)，其中 (x > 0) 且 (y > 0)。通过分离变量法， 我们可以得到 (dy/y = xdx)，进一步求解得到 (ln|y| = frac{1}{2}x^2 + C)，其 中 (C) 是积分常数。
常微分方程解法
常微分方程解法常微分方程是数学中的一门重要分支，研究描述自然界和社会现象中变化规律的方程。
解常微分方程的方法多种多样，下面将介绍常见的几种解法。
一、分离变量法分离变量法适用于形如dy/dx=f(x)g(y)的一阶常微分方程。
解题步骤如下：1. 将方程写成dy/g(y)=f(x)dx的形式，将变量分离。
2. 对两边同时积分，得到∫dy/g(y)=∫f(x)dx。
3. 左边的积分可以通过换元或者使用常见函数的积分公式进行计算。
4. 右边的积分可以通过与左边的积分结果进行比较来判断是否需要使用特殊的积分技巧。
5. 对左右两边同时积分后，解出方程中的积分常数。
6. 将积分常数代回原方程中，得到完整的解。
二、常数变易法常数变易法适用于形如dy/dx+p(x)y=q(x)的一阶常微分方程。
解题步骤如下：1. 先求出对应的齐次方程dy/dx+p(x)y=0的通解。
2. 假设原方程的特解为y=u(x)v(x)，其中u(x)是一个待定的函数，v(x)是齐次方程的通解。
3. 将y=u(x)v(x)代入原方程中，整理后得到关于u(x)和v(x)的方程。
4. 解出关于u(x)的方程，得到u(x)的值。
5. 将u(x)的值代入v(x)中，得到特解。
6. 特解与齐次方程的通解相加，即得到原方程的完整解。
三、二阶齐次线性方程解法二阶齐次线性方程的一般形式为d^2y/dx^2+p(x)dy/dx+q(x)y=0。
解题步骤如下：1. 求解对应的齐次方程d^2y/dx^2+p(x)dy/dx+q(x)y=0的特征方程r^2+p(x)r+q(x)=0，其中r为未知数。
2. 求解特征方程得到两个不同的根r1和r2。
3. 根据r1和r2的值，得到齐次方程的通解y=c1e^r1x+c2e^r2x，其中c1、c2为任意常数。
四、变量替换法变量替换法适用于形如dy/dx=f(y/x)的一阶常微分方程。
解题步骤如下：1. 进行变量替换，令u=y/x，即y=ux。
常微分方程的经典求解方法
R1 + R 2 α ≠ R1 R 2 c
v (0 ) = v (0 )
−
+
v0 (t) = Ae
R +R2 − 1 t R R2C 1
R2 E −αt + e R1 + R2 − R1R2cα
R2E v0 (0) = 0 = A+ R + R2 − R R2Cα 1 1 R2E ∴A = − R + R2 − R R2Cα 1 1
•经典法不足之处
•若微分方程右边激励项较复杂，则难以处理。 •若激励信号发生变化，则须全部重新求解。 •若初始条件发生变化，则须全部重新求解。 •这种方法是一种纯数学方法，无法突出系统响 应的物理概念。
•经典时域分析方法 微分方程的全解即系统的完全响应, 由齐次解 和特解组成 y (t ) = y (t ) + y (t )
h p
齐次解 yh (t) 的形式由齐次方程的特征根确定 特解 y p (t) 的形式由方程右边激励信号的形式 确定
齐次解yh(t)的形式
(1) 特征根是不等实根s1, s2, …, sn
−αt
dv0 (t ) R1 + R2 1 + v0 (t ) = e(t ) dt R1 R2 c R1c
−
R2
V0(t)
e(t)
Ae
R1 + R2 t R1R2C
因激励信号为 则：
u (t )
−αt
P46.表2—2若
R2E B= R +R2 −R R2cα 1 1
B ( t ) = Be
R1 + R2 −αt E −αt −αBe + Be = e R1R2c R1c
常微分方程(组)的数值解法
刚性常微分方程组求解
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 )
常微分方程常见形式及解法
()
A．江南制造总局的汽车
B．洋人发明的火车
C．轮船招商局的轮船
D．福州船政局的军舰
[解析] 由材料信息“19世纪七十年代，由江苏沿江居民 到上海”可判断最有可能是轮船招商局的轮船。
[答案] C
[题组冲关]
1．中国近代史上首次打破列强垄断局面的交通行业是 ( )
A．公路运输
B．铁路运输
C．轮船运输
常微分方程 毕文彬
4
简易微分方程的求解方法
01
一阶线性常微分方程
02
二阶常系数齐次常微分方程
常微分方程 毕文彬
5
01一阶线性常微分方程
对于一阶线性常微分方程，常用的方法是常 数变易法： 对于方程：
可知其通解：
然后将这个通解代回到原式中，即可求出 C(x)的值
常微分方程 毕文彬
6
02二阶常系数齐次常微分方程
齐次二阶线性微分方程：
描述谐振子的齐次二阶常系数线性微分方程：
非齐次一阶非线性微分方程：
描述长度为L的单摆的二阶非线性微分方程：
常微分方程 毕文彬
3
微分方程的解
微分方程的解通常是一个函数表达式（含一 个或多个待定常数，由初始条件确定）。例如： dy/dx=sinx， 的解是 y=-cosx+C， 其中C是待定常数； 例如，如果知道 y=f(π)=2， 则可推出 C=1， 而可知 y=-cosx+1，
出口价格同步变动的现象。与这一现象直接相关的近代事
业是
()
A．电报业
B．大众报业
C．铁路交通业 D．轮船航运业
[解析] 材料主要反映了信息交流的快捷，故选A。
[答案]A
[题组冲关]
常微分方程的常见解法
Euler折线法
近似导数
y(x0)
y(x1) h
y( x0 )
记为
y( x1 ) y( x0 ) hy( x0 ) y0 h f ( x0 , y0 )
解:设t时刻雪球的体积为
，表面积为 ，
由题得
球体与表面积的关系为
引入新常数
再利用题中的条件得
分离变量积分得方程得通解为
再利用条件 确定出常数C和r代入关系式得 t的取值在 之间。
方程为全微分方程的充要条件
定理2.1 设函数
和
在一个矩形区域
中连续且有连续的一阶偏导数，则
是全微分方程的充要条件为:
（2.3.3)
nan (x x0 )n1
f
x,
an
(
x
x0
)n
n1
n0
展开后比较两端同次幂的系数确定
an ,
y
y0
N n1
cn1 (x n
x0 )n
例：用待定系数法求
dy x2 y2 ,
的近似解。
dx
y(0) 1
解： 令 y a n (x x0 )n, 由 y (0) 1 得 a0 1 n0
([diff(y(x),x)=-y(x)],y(x),
# 定义微分方程
x=-2..2,
# 指定x范围
[[y(-2)=2],[y(-2)=1],[y(-2)=-2]], # 给出3个初始值
dirgrid=[17,17],
常微分方程的基本概念与解法
常微分方程的基本概念与解法常微分方程是数学中的一门重要分支，用于描述自然界中的各种变化规律。
本文将介绍常微分方程的基本概念和常见的解法。
一、常微分方程的概念常微分方程是关于未知函数的导数和自变量之间的关系式，其中自变量通常表示时间。
一般形式为dy/dx = f(x, y)，其中y是未知函数，f(x, y)是已知函数。
常微分方程可分为一阶常微分方程和高阶常微分方程两种。
1. 一阶常微分方程一阶常微分方程是指未知函数的导数只涉及到一阶导数的方程。
一阶常微分方程的一般形式为dy/dx = f(x, y)，也可以写成f(x, y)dx - dy = 0。
其中f(x, y)是已知函数，x是自变量，y是未知函数。
2. 高阶常微分方程高阶常微分方程是指未知函数的导数涉及到高阶导数的方程。
高阶常微分方程的一般形式为d^n y/dx^n = f(x, y, dy/dx, d^2 y/dx^2, ..., d^(n-1) y/dx^(n-1))，其中n为正整数，f是已知函数，x是自变量，y是未知函数。
二、常微分方程的解法解常微分方程的方法多种多样，根据方程的类型和特点选择不同的解法。
1. 可分离变量法当方程可以写成dy/dx = g(x)h(y)的形式时，可以使用可分离变量法解方程。
这种方法的关键是将变量分离，即将含有y的项移到方程的一边，含有x的项移到方程的另一边，然后分别积分得到x和y的表达式。
2. 线性常微分方程的求解线性常微分方程是指方程可以写成dy/dx + P(x)y = Q(x)的形式。
对于线性常微分方程，可以使用积分因子法求解。
首先找到一个函数u(x)，使得dy/dx + P(x)y = Q(x)乘以u(x)后变为全导数，则原方程可以写成d(uy)/dx = Q(x)u(x)的形式。
然后对等式两边进行积分并解得y的表达式。
3. 齐次线性常微分方程的求解齐次线性常微分方程是指方程可以写成dy/dx = f(y/x)的形式。
常系数微分方程组的解法
将高阶线性微分方程转化为幂级数形式，然后通过幂 级数的性质求解方程。
高阶非线性微分方程的解法
分离变量法
将非线性微分方程转化为多个一阶微分方程 ，然后分别求解。
迭代法
通过迭代公式逐步逼近非线性微分方程的解。
数值解法
利用数值计算方法求解非线性微分方程的近 似解，如欧拉法、龙格-库塔法等。
05
解决微分方程组对于理解复杂系统的 行为和预测未来发展趋势具有重要意 义。
常系数微分方程组的定义
常系数微分方程组是指方程中的系数 为常数的一类微分方程组。
常系数微分方程组的一般形式为 dy/dx = f(x, y)，其中 f(x, y) 是已知 的函数。
02
线性常系数微分方程组的解法
特征根法
总结词
神经传导
在神经传导过程中，微分方程组可以用来描述神 经信号的传递速度和传导通路的建立。
生态系统的稳定性
微分方程组可以用来分析生态系统的稳定性，如 物种之间的相互作用和生态平衡的维持。
THANKS
感谢观看
特征根法是一种通过解方程的特征方程来求解线性常系数微 分方程组的方法。
详细描述
特征根法的基本思想是，对于形如$y'' + py' + qy = 0$的一阶 线性常系数微分方程，通过求解其特征方程$lambda^2 + plambda + q = 0$，得到其特征根$lambda_1$和 $lambda_2$，然后利用这些特征根来求解原微分方程。
线性微分方程的方法。
02
通过将多个变量分离，可以将一个复杂的微分方程组
分解为多个简单的微分方程，从而简化求解过程。
03
第十章常微分方程（组）求解
第⼗章常微分⽅程（组）求解⼀、求微分⽅程的解（⼀）相关函数（命令）及简介1， dsolve('equ1','equ2',…):Matlab 求微分⽅程的解析解。
equ1,equ2,…为⽅程（或条件）。
写⽅程（或条件）时⽤Dy 表⽰y 关于⾃变量的⼀阶导数，⽤D2y 表⽰y 关于⾃变量的⼆阶导数，依次类推。
2， simplify(s):对表达式s 使⽤maple 的化简规则进⾏化简。
例如： syms xsimplify(sin(x)^2+cos(x)^2) ans=13,[r,how]=simple(s):由于Matlab 提供了多种化简规则，simple 命令就是对表达式s ⽤各种规则进⾏化简，然后⽤r 返回最简形式，how 返回形成这种形式所⽤的规则。
例如： syms x[r,how]=simple(cos(x)^2-sin(x)^2) r=cos(2*x) how=combine4,[T,Y]=solver(odefun,tspan,y0),求微分⽅程的数值解。
（1）其中的solver 为命令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb 之⼀。
（2）odefun 是显式常微分⽅程：00(,),().dyf t y dty t y ?==? （3）在积分区间0[,]f tspan t t =上，从0t 到f t ，⽤初始条件0y 求解。
（4）要获得问题在其他指定时间点012,,,...t t t 上的解，则令012[,,,...,]f tspan t t t t =（要求是单调的）（5）因为没有⼀种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采⽤不同的solver 。
求解器solver ODE 类型特点说明ode45 ⾮刚性单步算法：4，5阶的Runge-Kutta ⼤部分场合的⾸选算法⽅程；累计截断误差达3()x ?ode23 ⾮刚性单步算法：2，3阶的Runge-Kutta 使⽤于精度较低的情形⽅程；累计截断误差达3()x ?ode113 ⾮刚性多步法：Adams 算法；⾼低精度均计算时间⽐ode45短可到3610~10--。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
再用命令 >>clear; close; t=0:0.1:1; >>y=t+exp(-t); plot(t,y); %化解析解的图形 >>hold on; %保留已经画好的图形，如果下面再 画图，两个图形和并在一起 >>[t,y]=ode45('fun8',[0,1],1); >>plot(t,y,'ro'); %画数值解图形，用红色小圈画 >>xlabel('t'),ylabel('y')
线性常微分方程的解满足叠加原理,从而他们 的求解可归结为求一个特解和相应齐次微分方 程的通解.一阶变系数线性微分方程总可用这 一思路求得显式解。 高阶线性常系数微分方程可用特征根法求得 相应齐次微分方程的基本解，再用常数变异法 求特解。
一阶常微分方程与高阶微分方程可以互化，已给 一个 n 阶方程
的数值解.不妨取
则上面方程可化为 " 9.8 sin , (0) 15, ' (0) 0 解：先看看有没有解析解.运行MATLAB代码 >>clear; >>s=dsolve('D2y=9.8*sin(y)','y(0)=15' ,'Dy(0)=0','t') >>simplify(s) 知原方程没有解析解
运行MATLAB代码 >>clear; close; >>[t,y]=ode45('fun9',[0,10],[15,0]); >>plot(t,y(:,1)); %画随时间变化图,y(:2) 则表示的值 >>xlabel('t'),ylabel('y1')
16.5
16
y1
15.5 15 0
y ( n) f (t, y' , y",, y ( n1) )
设
y1 y, y2 y' ,, yn y (n1)，可将上式化为
y1 ' y 2 y ' y 2 3 y' y n n 1 y n ' f (t , y1 , y 2 , , y n )
一阶方程组
反过来，在许多情况下，一阶微分方程组也 可化为高阶方程。所以一阶微分方程组与高阶 常微分方程的理论与方法在许多方面是相通的， 一阶常系数线性微分方程组也可用特征根法求 解。
三、常微分方程的数值解法
除常系数线性微分方程可用特征根法求解， 少数特殊方程可用初等积分法求解外，大部分 微分方程无限世界，应用中主要依靠数值解法。 考虑一阶常微分方程初值问题
四、解微分方程的MATLAB命令
MATLAB中主要用dsolve求符号解析解， ode45,ode23,ode15s求数值解。
s=dsolve(‘方程1’, ‘方程2’,…,’初始条 件1’，’初始条件2’ …,’自变量’)
用字符串方程表示，自变量缺省值为t。导数用D表示，2阶导数 用D2表示，以此类推。S返回解析解。在方程组情形，s为一个符 号结构。
y' (t ) f (t , y(t )),t 0 t t f y(t 0 ) y0
其中 y ( y1 , y2 ,, ym )', f ( f1 , f 2 ,, f m )', y0 ( y10 , y20 ,, ym0 )'.
所谓数值解法，就是寻求 y(t ) 在一系列离散节点
求数值解
令
y1 , y2 '
可将原方程化为如下方程组
y1 ' y 2 y 2 ' 9.8 sin( y1 ) y (0) 15, y (0) 0 2 1
建立M文件fun9.m如下 %M文件fun9.m function f=fun9(t,y) f=[y(2), 9.8*sin(y(1))]'; %f向量必须为一 列向量
若上式中的系数 无关，称之为常系数。
ai (t ),i 1,2,, n
均与
t
二、常微分方程的解析解
有些微分方程可直接通过积分求解.例如，一
dy dy dt ， y 1 可化为 阶常系数常微分方程 y 1 dt
两边积分可得通解为
y cet 1
。其中 为任意常数。
c
有些常微分方程可用一些技巧,如分离变量法, 积分因子法,常数变异法,降阶法等可化为可积分 的方程而求得解析解.
数学实验
第十章 常微分方程（组） 的求解
一、微分方程的概念
未知的函数以及它的某些阶的导数连同自变 量都由一已知方程联系在一起的方程称为微 分方程。 如果未知函数是一元函数，称为常微分方程。 常微分方程的一般形式为
F (t, y, y' , y",, y ) 0
( n)
如果未知函数是多元函数，成为偏微分方 程。
t 0 t1 t n t f 上的近似值 y k , k 0,1,, n ，称
hk t k 1 t k 为步长，通常取为常量
h
， 最简单的
数值解法是Euler法。
Euler法的思路极其简单：在节点出用差商 近似代替导数
y (t k 1 ) y (t k ) y ' (t k ) h
这样导出计算公式（称为Euler格式）
yk 1 yk hf (t k , yk ), k 0,1,2,
他能求解各种形式的微分方程。Euler法也称折线法。
Euler方法只有一阶精度，改进方法有二阶 Runge-Kutta法、四阶Runge-Kutta法、五 阶Runge-Kutta-Felhberg法和先行多步法 等，这些方法可用于解高阶常微分方程（组） 初值问题。边值问题采用不同方法，如差分法、 有限元法等。数值算法的主要缺点是它缺乏物 理理解。
方程（2）求解的MATLAB代码为： >>clear; >>s=dsolve('D2y=sin(2*x)y','y(0)=0','Dy(0)=1','x') >>simplify(s) %以最简形式显示s。
结果为 s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x) ans =-2/3*sin(x)*cos(x)+5/3*sin(x)。
[tout,yout]=ode45(‘yprime’,[t0,tf],y0)
采用变步长四阶Runge-Kutta法和五阶Runge-Kutta-Felhberg法 求数值解，yprime是用以表示f(t,y)的M文件名，t0表示自变量的 初始值，tf表示自变量的终值，y0表示初始向量值。输出向量 tout表示节点(t0,t1, …,tn)T,输出矩阵yout表示数值解，每一 列对应y的一个分量。若无输出参数，则自动作出图形。
例２ 求解微分方程
先求解析解，再求数值解，并进行比较。 y' y t 1, y(0) 1, 解：由 >>clear; >>s=dsolve('Dy=-y+t+1','y(0)=1','t') >>simplify(s) 可得解析解为 y t e t
下面再求其数值解，先编写M文件fun8.m %M函数fun8.m function f=fun8(t,y) f=-y+t+1;
ode45是最常用的求解微分方程数值解的命 令，对于刚性方程组不宜采用。ode23与 ode45类似，只是精度低一些。ode12s用来 求解刚性方程组，是用格式同ode45。可以用 help dsolve, help ode45查阅有关这些命令 的详细信息.
例1 求下列微分方程的解析解
（1） （3）
方程（3）求解的MATLAB代码为： >>clear; >>s=dsolve('Df=f+g','Dg=gf','f(0)=1','g(0)=1') >>simplify(s.f) %s是一个结构 >>simplify(s.g)
结果为 ans =exp(t)*cos(t)+exp(t)*sin(t) ans =-exp(t)*sin(t)+exp(t)*cos(t)
联系一些未知函数的一组微分方程组称为微分方程组。 微分方程中出现的未知函数的导数的最高阶解数称为 微分方程的阶。 若方程中未知函数及其各阶导数都是一次的，称为线 性常微分方程，一般表示为
y ( n) a1 (t ) y ( n1) an1 (t ) y'an (t ) y b(t )
1
2
3
4
5 t
6
7
8
9
10
实验内容
书上课后： 1、（2）（4）（6）（8）
2、
1.4 1.35 1.3 1.25 1.2 1.15 1.1 1.05 1
y
0
0.2
0.4 t
0.6
0.8
1
由图10.1可见，解析解和数值解吻4; mgsin , (0) 0 , ' (0) 0