下载文档原格式
115§6常微分方程的求解一、知识背景常微分方程在物理学，工程技术中运用非常广泛，相当重要。
许多物质运动的过程用常微分方程来描述，如：质点的加速运动、谐振动、电容充放电过程及电感通电断电等过程，因此，求解常微分方程成为很多物理问题求解的一种常用方法。
而有时很难求出解析解，但可求出常微分方程的数值解以逼近解析解，以完成对摸型的研究。
求解常微分方程数值解通常有欧拉法和龙格-贝塔法。
欧拉法解法的基本思想是在小区间上用差商代替导数，并通过把小区间不断地划分求极限，从而最终得到数值解。
龙格-贝塔方法基本思想也是用差商代替导数，但是该方法是在小区间内运用了微分中值定理，在小区间内多取点，再取加权平均值来构造精度更高的计算式。
在MATLAB 系统中，主要采用龙格-贝塔法来计算常微分方程的数值解。
图 6-10d d =+c q tq R或 cqt q R-=d d 我们也可写成右面方程组形式：),(d d t q f tq =，00)(q t q =这即是一阶常微分方程初值问题的一般形式。
二、 计算指令 —— ode23，ode45语句格式（以ode23为例）：[ t, y ] = ode23 (‘f ’,tspan,y0,tol )语句中各符号意义如下f ：求解的常微分方程的文件名，把方程写成函数形式并存储于m 文件中。
方程形式为y'=f(t,y)。
举例：如图6-1中所示电路，先将开关k 掷向“1”端，待电容器c 充完电后，将开关k 掷向“2”端，电容器开始放电。
放电过程满足下面方程：116tspan ：输入[t0,tf],分别为自变量的初始值和最终值，为单调递增（减）的积分区间。
y0：函数的初始值矢量。
Tol ：误差范围，(缺省值为0.000001)[t,f]：t 是输出的时间列矢量，矩阵y 的每个列矢量是解的一个分量。
例1：用求数值解方法,求解常微分方程：3'x x -=，初始值x(0)=1。
第7章常微分方程初值问题数值解法本章探讨常微分方程特解的常用数值方法的构造和原理，主要介绍求常微分方程初值问题的常用方法和有关知识。
重点论述Euler方法、Runge-Kutta方法和线性多步法的原理、构造、局部截断误差和稳定性等内容。
1201217.1 实际案例工程技术里某些振动问题可以表示为单摆的运动，其运动规律的微分方程为：sin 0(0)30(0)0gl θθθθ⎧+=⎪⎪⎨=⎪⎪=⎩怎样求出其特解()t θ？该微分方程不能用通常的解析方法来求解！ 怎样解不能用解析方法求解的微分方程特解问题是本章要解决的问题。
1227.2 问题的描述和基本概念1、常微分方程初值问题 ● 一般形式0(,)()y f x y y a y '=⎧⎨=⎩()a x b << （7.1） 式中(,)f x y 已知，0()y a y =称为初值条件。
● 初值问题的数值方法和数值解求函数()y y x =在若干离散点k x 上的近似值(0,1,,)k y k n = 的方法称为初值问题的数值方法，而称(0,1,,)k y k n = 为初值问题的数值解。
1232. 建立数值解法的思想与方法微分方程初值问题的数值解法是用离散化方法将初值问题化为差分方程后再求解的方式。
设区间[a,b ]上的一组节点为01n a x x x b =<<<=距离1k k k h x x +=-称为步长。
求数值解一般是从0y 开使逐次顺序求出12,,y y 。
初值问题的解法有单步法和多步法两种。
单步法：计算1k y +时只用到k y 一个值；多步法：计算1k y +时要用1,,,k k k l y y y +- 多个值。
数值解法的显格式和隐格式。
124● 微分方程基本的离散化方法数值微分法，数值积分法和Taylor 展开法。
● 常微分方程初值问题化为差分方程1、 用离散方法去掉方程中的导数'y 得到近似离散化方程；2、 在近似离散化方程中用k y 代替()k y x ；3、 在近似离散化方程中将近似号“ ”用等号“=”代替。
1251) 数值微分法由初值问题（7.1）有'()(,())k k k y x f x y x =，用数值微分的2点前差公式代替'()y x ，得近似离散化方程记1k k h x x +=-，做()k k y y x →，“”，得差分方程1(,)k kk k y y f x y h+-= 写容易计算的形式1(,)k k k k y y hf x y +=+ (0,1,2,,1)k n =- （7.2） （Euler 公式）由初值条件0()y y a =及式（7.2）可求出（7.1）的数值解12,,,n y y y 。
公式（7.2）是显式单步法。
1262）数值积分法在1[,]k k x x +上对'(,)y f x y =两边取定积分，得111()()'(,())k k kkx x k k x x y x y x y dx f x y x dx +++-==⎰⎰对右端积分采用梯形公式（数值积分公式）得近似离散化方程：于是得到求初值问题（7.1）的梯形方法(0,1,,1)k n =-该公式是隐式单步法。
1273）Taylor 展开法因为初值问题中函数(,)f x y 是已知函数，由(,)y f x y '=，可以计算''y ，'''y ，…， 于是有函数()y y x =在x 处的Taylor 展式取上式右端前若干项，得近似离散化方程。
例如取前两项有1()()(,())k k k k y x y x hf x y x +≈+于是又得到Euler 公式：1(,)k k k k y y hf x y +=+1287.3 数值解法的误差、阶与绝对稳定性单步法数学描述为 显式：1(,,)k k k k y y h x y h ϕ+=+ （7.4）隐式：11(,,,)k k k k k y y h x y y h ϕ++=+ （7.5）其中ϕ为与(,)f x y 有关的函数，称为增量函数。
129显式单步法的一些概念定义7.1 设()y x 是问题（7.1）的解，k y是经过式（7.4）求出的的计算解，则称111()k k k e y x y+++=- 为显式单步法（7.4）在节点1k x +的整体截断误差，而称11()()(,(),)k k k k k T y x y x h x y x h ϕ++=-- （7.6） 为在1k x +点的局部截断误差。
()k y x 表示解()y x 在k x 的值，是准确值，没有误差；k y 表示由求数值解公式得出()k y x 的近似值，是数值解，有截断误差；k y 表示用计算机计算k y 给出的计算解，有舍入误差。
130局部截断误差1k T +的理解假设在计算()k y x 时没有误差（()k k y y x =）下，由式（7.4）计算出的1k y +（1()(,(),)k k k k y y x h x y x h ϕ+=-）与1()k y x +的误差()111k k k T y x y +++=-。
整体截断误差还要加上()k y x 与k y 的误差。
考察初值问题解法的优劣，引入阶的概念。
定义7.2 如果初值问题（7.1）对某种数值解法的局部截断误差为11()P k T O h ++=则称该方法具有p 阶精度或该方法是p 阶方法。
方法的阶越高，方法越好。
131主局部截断误差或局部截断误差的主项如果某方法是p 阶方法，若其局部截断误差11()P k T O h ++=按h 展开为1121()(,())()P P P k k k T O h g x y x h O h ++++==+则称1(,())P k k g x y x h +为主局部截断误差。
在同阶方法中，主局部截断误差越小，方法越好。
例7.1常微分方程初值问题00(,)()y f x y y x y '=⎧⎨=⎩的单步法为1111[(,)2(,)],3n n n n n n k k hy y f x y f x y h x x ++++=++=-试求其局部截断误差主项并回答它是几阶精度的？132解 该单步公式的局部截断误差是1111()()[(,())2(,())]3n n n n n n n hT y x y x f x y x f x y x ++++=--+()11()()()2()3n n n n hy x y x y x y x ++''=--+234(4)()()()()()2!3!4!n n n n n h h h y x hy x y x y x y x ''''''=+++++22()()()()()332!n n n n n h h y x y x h y x hy x y x ⎛⎫'''''''---+++ ⎪⎝⎭231212(1)()()()()3323n n hy x h y x O h '''=--+-+2321()()()6n h y x O h O h ''=-+=故局部截断误差主项是)(61''2n x y h -，方法是一阶的。
133求阶p 的另一方法 因为1()()(,(),)k k k k k T y x h y x h x y x h ϕ+=+--去掉下标，有()()()(,(),)T x y x h y x h x y x h ϕ=+--若将()y x h +在x 点展开有1()()P T x O h +=则知该方法的阶是p 。
134例如，对Euler 方法1(,)k k k k y y hf x y +=+，有(,(),)(,())x y x h f x y x ϕ=那么()()()(,())T x y x h y x hf x y x =+--将()y x h +在x 点展开，有2()()'()''()2!h y x h y x hy x y x +=+++2()(,())''()2!h y x hf x y x y x =+++故有22()''()()2h T x y x O h =+=因此，Euler 方法是一阶方法。
135定义7.3 设用某种数值方法求初值问题（7.1）在任意节点k x 的数值解k y 时，则称该数值方法是绝对稳定的。
这里k ε是计算机计算k y 时得出的计算解k y的舍入误差，k k k yy ε=+ 。
通常用试验方程'y y λ= （λ为复数）来讨论求解初值问题的数值方法绝对稳定性；对具体初值问题，可取(,)k k x y fyλ∂=∂。
稳定性常与步长h 有关。
136定义7.4某方法在试验方程中绝对稳定的复平面h μλ=范围称为该方法的绝对稳定域，它与复平面实轴的交称为该方法的绝对稳定区间。
绝对稳定域包含复平面左半平面(){}|Re 0h h λλ<的方法称为是A-稳定的。
绝对稳定域越大，对应的方法绝对稳定性越好。
1377.4 Euler 方法的有关问题用Euler 公式1(,)k k k k y y hf x y +=+，0()y a y =求解初值问题（7.1）数值解的方法称为Euler 方法。
1 )Euler 方法的几何意义Euler 方法常称为折线法。
2) Euler 方法的误差设1k y+ 为1k y +的计算解，(,)f x y 满足 (,())(,)()k k k k k k f x y x f x y L y x y -≤-其中0(0,1,2,,)L k n >= 。
则有Euler 方法的局部截断误差2()k T O h =；138Euler 方法的总体截断误差1111111()()k k k k k k k e y x yy x y y y +++++++=-=-+- 由()()()1,k k k k y y x hf x y x +=+，1(,)k k k k yy hf x y +=+ ,有 11(1)k k k e T hL e ++≤++2111(1(1)kk k T T hL T +-≤++++因为对任意m2110(1)()(1)k km m k k m m m e hL T O h hL ++-==≤+=+∑∑12(1)1()()11k hL O h O h hL ++-==+- 由k 的任意性，可知Euler 方法的总体截断误差为()k e O h = 当h 足够小时，由Euler 方法计算所得数值解可以很好地逼近准确解，从而Euler 方法是收敛的。
