数 值 分 析第5章 常微分方程数值解法1111111111112222222222223333333333331 引言v1.0 基本概念1. 常微分方程的初值问题：称为具有初值(1.2)的常微分方程. 若f(x,y)在axb, |y|+上连续，且关于y满足Lip条件：常数L使| f(x, y1) f(x, y2)| L|y1 y2|则初值问题(1.1)(1.2)存在唯一连续可微解y(x).注：以下总假设f 满足Lip条件.1111111111112222222222223333333333331 引言v1.0 基本概念1. 常微分方程的初值问题：称为具有初值(1.2)的常微分方程. (1.1)(1.2)等价于微分方程： (1.3)注：一般无初等解(解析解)，即使有形式也复杂.1111111111112222222222223333333333331 引言v1.0 基本概念2. 初值问题的数值解 设(1.1)(1.2)的解y(x)在节点xi处的近似解值为 yi y(xi), a x1 x2 xn = b则称yi (i = 1, 2, , n)为(1.1)(1.2)的数值解，又称y(xi)的计算值.1111111111112222222222223333333333331 引言v1.0 基本概念3. 数值方法 两种转化： 由微分出发的数值方法. 由积分 出发的数值方法. 计算方法 步进法：从初始条件出发，逐步求y1, y2, , yn. 又有两种：单步法，多步法.注：采用等距节点：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式. (1.6)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.1. 前进欧拉公式 (1.6)的前半部分为：令 yi+1 = yi + hf(xi, yi) (1.7)其中yi = y(xi) , 则yi+1 y(xi+1)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.1. 前进欧拉公式 令 yi+1 = yi + hf(xi, yi) (1.7)其中yi = y(xi) , 则yi+1 y(xi+1)记 (1.8)则称(1.7)为前进欧拉求解公式. 简称为欧拉公式或欧拉法. (1.8)称为欧拉公式的余项：ei+1(h) = y(xi+1) yi+1 1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.2. 后退欧拉公式 (1.6)的后半部分令 yi+1 = yi + hf(xi+1, yi+1) (1.9)其中yi = y(xi), 则yi+1 y(xi+1) 1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.2. 后退欧拉公式令 yi+1 = yi + hf(xi+1, yi+1) (1.9)其中yi = y(xi), 则yi+1 y(xi+1) 注：(1.9)中f(xi+1, yi+1) f(xi+1, y(xi+1) 余项 (1.10)1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.2. 后退欧拉公式令 yi+1 = yi + hf(xi+1, yi+1) (1.9)其中yi = y(xi), 则yi+1 y(xi+1) 注： 称(1.9)为后退欧拉公式(后退欧拉法). 称(1.10)为后退欧拉法的误差近似值. 欧拉法与后退欧拉公式的区别：(1.7)为直接计算公式称显式公式.(1.9)为关于函数方程称隐式公式.1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.【例1】取h=0.1求解初值问题： (1.11).解： ，xi = ih = 0.1i, (i = 0, 1, 2, , 10) 欧拉法：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.【例1】取h=0.1求解初值问题： (1.11).解： ，xi = ih = 0.1i, (i = 0, 1, 2, , 10) 后退欧拉法： 1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.注：为避免求解函数方程，采用显式与隐式结合的方法： 此方法称为 预测校正系统. 求解过程为：1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.预测校正系统：【例2】利用预测校正系统求解例1.1111111111112222222222223333333333331 引言v1.1 基于数值微分的求解公式.预测校正系统：注：显式比隐式方便，但有时隐式效果比显式好.(4介绍).1111111111112222222222223333333333331 引言v1.2 截断误差定义1.1 称ek(h) = y(xk) yk为计算yk的公式第k步的局部截断误差. 注：“局部”是指在计算第k步时，假定前面yi = y(xi) (i k).而yk y(xk) 欧拉法. 后退欧拉法.一般根据y(xk)对y(k), y(k)做估计.1111111111112222222222223333333333331 引言v1.2 截断误差定义1.2 设ei(h) (i = 1, 2, , k)为求解公式第i步的局部截断误差.称为该求解公式在点上的整体截断误差.注：局部截断误差ek(h)与yk有关. 整体截断误差Ek(h)与y1, y2, , yk有关.所有ek(h)都与h有关.1111111111112222222222223333333333331 引言v1.2 截断误差定义1.3 若局部截断误差e(h)=O(hp+1)，则称该求解公式具有p阶精度.注：欧拉法具有一阶精度.(精度越高越好)1111111111112222222222223333333333331 引言作业 P208 1，2，3.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式 (1.13)若已知y(xk) = yk, 则计算积分可求出y(xk+1) . 如用矩形公式求积分则有y(xk+1) = y(xk) + hf(xk, yk)令yk+1 = y(xk) + hf(xk, yk)即为欧拉公式. 故欧拉公式又称矩形法.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式 (1.13)考虑1. 梯形公式记 (1.14)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式1. 梯形公式记 (1.14)称(1.14)为梯形(求解)公式. 简称梯形法.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式1. 梯形公式梯形(求解)公式, 简称梯形法: (1.14)注：梯形公式的余项： 故是二阶精度.111111111111222222222222333333333333v1.3 基于数值积分的求解公式1. 梯形公式 (1.14) 梯形公式为隐式公式.预测校正系统 (1.15)称(1.15)为改进的欧拉公式，也可记为1 引言1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式1. 梯形公式 (1.14) 可以证明，改进欧拉公式也具有二阶精度.1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解取h=0.1.计算到x=0.5.解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(Euler法) 求解公式：yk =yk1+h(xk1yk1+1)= hxk1+(1 h)yk1 + h = 0.1xk1+0.9yk1+0.1 1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(梯形法)求解公式：yk=yk1+h(xk1yk1+1)+(xkyk+1)/2解出yk，得方程1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式【例3】用欧拉法，梯形法以及改进欧拉法求解解：f(x, y) = xy + 1, a = x0 = 0, b = 0.5, y0 = 1, n = 5(改进Euler法)求解公式：yk=yk1+h(xk1yk1+1) + xk (yk +h(xkyk+1)+1/2得=0.905yk1+0.045xk1+0.05xk+0.095方程1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式2. 辛卜生公式 记 (1.17)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式2. 辛卜生公式记 (1.17)其余项1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式2. 辛卜生公式记 (1.17)将xk1, xk 对分： 调整下标为xi2, xi ：xi2 = xk1, xi1 = xk1+h1, xi = xk1+2h1= xk则(1.17)化为 (1.19)称(1.19)为辛卜生求解公式，其中fk2= f(xk2, y(xk2)，fk1 = f(xk1, y(xk1)，fk = f(xk, y(xk)1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式2. 辛卜生公式记 (1.17) (1.19)称(1.19)为辛卜生求解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注： (1.19)的误差：1111111111112222222222223333333333331 引言v1.3 基于数值积分的求解公式2. 辛卜生公式记 (1.17) (1.19)称(1.19)为辛卜生求解公式，其中fi2= f(xi2, y(xi2)，fi1 = f(xi1, y(xi1)，fi = f(xi, y(xi)注： 隐式(需显化)多步将在3中讨论.1111111111112222222222223333333333332 Runge - Kutta法v2.0 原理 其中K = f(, y() = y()称为y在xi1, xi上的平均斜率.欧拉法：改进欧拉法：(2.1)1111111111112222222222223333333333332 Runge - Kutta法v2.0 原理 其中K = f(, y() = y()称为y在xi1, xi上的平均斜率.对(1.17)显化：辛卜生： (2.4)1111111111112222222222223333333333332 Runge - Kutta法v2.0 原理其中K = f(, y() = y()称为y在xi1, xi上的平均斜率.设想：在中多计算(预测)几个点上的值然后可加权取平均值作为的近似值可能构成更高阶的公式.一阶二阶三阶1