常微分方程数值解法课件.pptx

1 引言1.0 基本概念1.常微分方程的初值问题：称为具有初值(1.2)的常微分方程.若f(x,y)在a x b,|y|+上连续，且关于y满足Lip条件：常数L使|f(x,y1)f(x,y2)|L|y1 y2|则初值问题(1.1)(1.2)存在唯一连续可微解y(x).注：以下总假设f 满足Lip条件.第1页/共73页1 引言1.0 基本概念1.常微分方程的初值问题：称为具有初值(1.2)的常微分方程.(1.1)(1.2)等价于微分方程：(1.3)注：一般无初等解(解析解)，即使有形式也复杂.第2页/共73页1 引言1.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)的计算值.第3页/共73页1 引言1.0 基本概念3.数值方法 两种转化：由微分出发的数值方法.由积分 出发的数值方法.计算方法 步进法：从初始条件出发，逐步求y1,y2,yn.又有两种：单步法，多步法.注：采用等距节点：第4页/共73页1 引言1.1 基于数值微分的求解公式.(1.6)第5页/共73页1 引言1.1 基于数值微分的求解公式.1.前进欧拉公式 (1.6)的前半部分为：令 yi+1=yi+hf(xi,yi)(1.7)其中yi=y(xi),则yi+1 y(xi+1)第6页/共73页1 引言1.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 第7页/共73页1 引言1.1 基于数值微分的求解公式.2.后退欧拉公式 (1.6)的后半部分令 yi+1=yi+hf(xi+1,yi+1)(1.9)其中yi=y(xi),则yi+1 y(xi+1)第8页/共73页1 引言1.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)第9页/共73页1 引言1.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)为关于函数方程称隐式公式.第10页/共73页1 引言1.1 基于数值微分的求解公式.【例1】取h=0.1求解初值问题：(1.11).解：，xi=ih=0.1 i,(i=0,1,2,10)欧拉法：第11页/共73页1 引言1.1 基于数值微分的求解公式.【例1】取h=0.1求解初值问题：(1.11).解：，xi=ih=0.1 i,(i=0,1,2,10)后退欧拉法：第12页/共73页1 引言1.1 基于数值微分的求解公式.注：为避免求解函数方程，采用显式与隐式结合的方法：此方法称为 预测校正系统.求解过程为：第13页/共73页1 引言1.1 基于数值微分的求解公式.预测校正系统：【例2】利用预测校正系统求解例1.第14页/共73页1 引言1.1 基于数值微分的求解公式.预测校正系统：注：显式比隐式方便，但有时隐式效果比显式好.(4介绍).第15页/共73页1 引言1.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)做估计.第16页/共73页1 引言1.2 截断误差定义1.2 设ei(h)(i=1,2,k)为求解公式第i步的局部截断误差.称为该求解公式在点上的整体截断误差.注：局部截断误差ek(h)与yk有关.整体截断误差Ek(h)与y1,y2,yk有关.所有ek(h)都与h有关.第17页/共73页1 引言1.2 截断误差定义1.3 若局部截断误差e(h)=O(hp+1)，则称该求解公式具有p阶精度.注：欧拉法具有一阶精度.(精度越高越好)第18页/共73页1 引言作业 P208 1，2，3.第19页/共73页1 引言1.3 基于数值积分的求解公式 (1.13)若已知y(xk)=yk,则计算积分可求出y(xk+1).如用矩形公式求积分则有y(xk+1)=y(xk)+hf(xk,yk)令yk+1=y(xk)+hf(xk,yk)即为欧拉公式.故欧拉公式又称矩形法.第20页/共73页1 引言1.3 基于数值积分的求解公式 (1.13)考虑1.梯形公式记 (1.14)第21页/共73页1 引言1.3 基于数值积分的求解公式1.梯形公式记 (1.14)称(1.14)为梯形(求解)公式.简称梯形法.第22页/共73页1 引言1.3 基于数值积分的求解公式1.梯形公式梯形(求解)公式,简称梯形法:(1.14)注：梯形公式的余项：故是二阶精度.第23页/共73页1.3 基于数值积分的求解公式1.梯形公式 (1.14)梯形公式为隐式公式.预测校正系统 (1.15)称(1.15)为改进的欧拉公式，也可记为1 引言第24页/共73页1 引言1.3 基于数值积分的求解公式1.梯形公式 (1.14)可以证明，改进欧拉公式也具有二阶精度.第25页/共73页1 引言1.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 第26页/共73页1 引言1.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，得方程第27页/共73页1 引言1.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方程第28页/共73页1 引言1.3 基于数值积分的求解公式2.辛卜生公式 记 (1.17)第29页/共73页1 引言1.3 基于数值积分的求解公式2.辛卜生公式记 (1.17)其余项第30页/共73页1 引言1.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)第31页/共73页1 引言1.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)的误差：第32页/共73页1 引言1.3 基于数值积分的求解公式2.辛卜生公式记 (1.17)(1.19)称(1.19)为辛卜生求解公式，其中fi2=f(xi2,y(xi2)，fi1=f(xi1,y(xi1)，fi=f(xi,y(xi)注：隐式(需显化)多步将在3中讨论.第33页/共73页2 Runge-Kutta法2.0 原理 其中K =f(,y()=y()称为y在xi1,xi上的平均斜率.欧拉法：改进欧拉法：(2.1)第34页/共73页2 Runge-Kutta法2.0 原理 其中K =f(,y()=y()称为y在xi1,xi上的平均斜率.对(1.17)显化：辛卜生：(2.4)第35页/共73页2 Runge-Kutta法2.0 原理其中K =f(,y()=y()称为y在xi1,xi上的平均斜率.设想：在中多计算(预测)几个点上的值然后可加权取平均值作为的近似值可能构成更高阶的公式.一阶二阶三阶第36页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式 (*)其中0 j 1，yi1+jh是y(xi1+jh)的预测值.称(*)为R-K公式注：(2.1)(2.4)分别称为二阶，三阶R-K公式.j，j，j为待定系数.使(*)的阶数尽量高.第37页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式参数的确定，以m=2为例.欲求 1，2，2.原则:使ei(h)=y(xi)yi的阶数尽可能高第38页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式展开展开 原则:使ei(h)=y(xi)yi的阶数尽可能高第39页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式 原则:使ei(h)=y(xi)yi的阶数尽可能高第40页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式 欲求截断误差ei(h)=y(xi)yi关于h的阶数尽可能高，应使无穷多解，从而有许多2阶R-K公式第41页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式应使注：取 1=2=1/2，2=1，即为改进欧拉公式.第42页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式应使注：取 1=0，2=1，2=1/2，即为中点公式第43页/共73页2 Runge-Kutta法2.1 Runge-Kutta公式应使注：二阶R-K公式的截断误差为故为二阶方法.相仿可得更高阶的R-K公式.第44页/共73页2 Runge-Kutta法2.2 经典R-K公式 在4解R-K公式中最重要的是经典R-K公式.(2.6)注：(2.6)为4阶方法.第45页/共73页2 Runge-Kutta法2.2 经典R-K公式 在4解R-K公式中最重要的是经典R-K公式.(2.6)注：R-K法对4阶以上不一定能提高整数阶.第46页/共73页2 Runge-Kutta法2.2 经典R-K公式【例4】使用三阶，四阶R-K法求解初值问题：的部分计算值y1，y2，y3，其中h=0.1.解 使用三阶R-K法第47页/共73页2 Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：的部分计算值y1，y2，y3，其中h=0.1.解 使用三阶R-K法第48页/共73页2 Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：的部分计算值y1，y2，y3，其中h=0.1.解 使用四阶R-K法第49页/共73页2 Runge-Kutta法【例4】使用三阶，四阶R-K法求解初值问题：的部分计算值y1，y2，y3，其中h=0.1.解 使用四阶R-K法第50页/共73页2 Runge-Kutta法注 使用R-K法要求具备较好的光滑性，否则效果不如低阶的.作业P209 8 9，10.第51页/共73页3 线性多步法单步法的优点：简单，计算yk+1只用yk.缺点:没有充分利用前面的信息且计算y(xk+h)较困难回顾Simpson：(1.19)考虑：(3.1)两种插值求积：将xk1,xk增加内部节点，改为xk2,xk导出的公式称为闭型求解公式.线性多步第52页/共73页3 线性多步法考虑：(3.1)两种插值求积：将xk1,xk增加内部节点，改为xk2,xk导出的公式称为闭型求解公式.在xk1,xk外增加插值节点，导出的公式称为开型求解公式.开型有显和隐，闭型也有显和隐.第53页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式 取节点xk3,xk2,xk1，在xk3,xk上作F(x)=f(x,y(x)的插值多项式.第54页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式 取节点xk3,xk2,xk1，在xk3,xk上记xki=xk ih,x=xk+th，则 第55页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式 取节点xk3,xk2,xk1，记xki=xk ih,x=xk+th，则 代入(3.1)得第56页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式 取节点xk3,xk2,xk1，记xki=xk ih,x=xk+th，则 令 (3.4)称(3.4)为亚当斯显式求解公式(线性多步).第57页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式 取节点xk3,xk2,xk1，记xki=xk ih,x=xk+th，则余项：从而(3.4)具有3阶精度.称为3阶亚当斯求解公式.第58页/共73页3 线性多步法3.1 开型求解公式1.亚当斯显式求解公式类似地取xk4,xk3,xk2,xk1 在xk4,xk上作F(x)=f(x,y(x)的插值多项式，可导出4阶亚当斯显式求解公式：(3.6)(3.7)4阶精度第59页/共73页3 线性多步法3.1 开型求解公式2.亚当斯隐式求解公式 取xk3,xk2,xk1,xk，在xk3,xk上作F(x)=f(x,y(x)的插值多项式用上述方法可导出：(3.8)(3.9)称为亚当斯隐式求解公式.第60页/共73页3 线性多步法3.1 开型求解公式2.亚当斯隐式求解公式 (3.8)(3.9)称为亚当斯隐式求解公式.注：利用4阶公式(3.6)显化之：(3.10)称(3.10)为亚当斯预测校正系统.第61页/共73页3 线性多步法3.2 闭型求解系统 将xk1,xk扩充为xk4,xk，取xk4，xk3，xk2，xk1为节点，作F(x)=f(x,y(x)的牛顿前插多项式.则 第62页/共73页3 线性多步法3.2 闭型求解系统 将xk1,xk扩充为xk4,xk，取xk4，xk3，xk2，xk1为节点，作F(x)=f(x,y(x)的牛顿前插多项式.则 令x=xk+(t 4)h 则 第63页/共73页3 线性多步法3.2 闭型求解系统令x=xk+(t 4)h 则由第64页/共73页3 线性多步法3.2 闭型求解系统令 (3.11)称(3.11)为米尔恩求解公式(Miline).余项：第65页/共73页3 线性多步法3.2 闭型求解系统 (3.11)称(3.11)为米尔恩求解公式(Miline).余项：分别在x=xk处展开y(xk 4h)和y(xk ih)为台劳级数并整理得 (3.12)第66页/共73页3 线性多步法3.2 闭型求解系统 (3.11)余项：(3.12)注：辛卜生公式求解 (3.13)第67页/共73页3 线性多步法3.2 闭型求解系统 (3.11)余项：(3.12)注：多步法的开始须用单步法，使用高精度多步法时，应用高精度单步法计算开始值.第68页/共73页3 线性多步法3.2 闭型求解系统【例5】初值问题设h=0.2.已知y1=0.181，分别用3阶亚当斯和米尔恩公式计算x3=0.6和x4=0.8的计算值.解 用3阶亚当斯计算y3 y(0.6)3阶亚当斯公式：第69页/共73页3 线性多步法3.2 闭型求解系统【例5】初值问题设h=0.2.已知y1=0.181，分别用3阶亚当斯和米尔恩公式计算x3=0.6和x4=0.8的计算值.解 用3阶亚当斯计算y3 y(0.6)先用3阶R-K法求y2 y(0.4)，然后求y3第70页/共73页3 线性多步法3.2 闭型求解系统【例5】初值问题设h=0.2.已知y1=0.181，分别用3阶亚当斯和米尔恩公式计算x3=0.6和x4=0.8的计算值.解 用米尔恩求y4 y(0.8)4阶米尔恩公式第71页/共73页3 线性多步法3.2 闭型求解系统【例5】初值问题设h=0.2.已知y1=0.181，分别用3阶亚当斯和米尔恩公式计算x3=0.6和x4=0.8的计算值.解 用米尔恩求y4 y(0.8)4阶米尔恩公式第72页/共73页感谢您的观看！第73页/共73页