第5讲：常微分方程数值解

常微分方程数值解《数值分析》第五讲第五章：常微分方程数值解§5.1引言1、常微分方程与解为n阶常微分方程。如果函数在区间[a,b]内n阶可导，称方程满足方程的函数称为微分方程的解。则如为任意常数）为方程的解,一般称为方程的通解。如果则有为方程满足定解条件的解。第五章：常微分方程数值解方程的通解满足定解条件的解微分关系（方程）解的图示第五章：常微分方程数值解本教材重点讨论定解问题(初值问题）定解条件（初始条件）是否能够找到定解问题的解取决于仅有极少数的方程可以通过“常数变易法”、“可分离变量法”等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解。如等等2、数值解的思想第五章：常微分方程数值解（1）将连续变量离散为（2）用代数的方法求出解函数在点的近似值*数学界关注工程师关注如果找不到解函数数学界还关注：解的存在性解的唯一性解的光滑性解的振动性解的周期性解的稳定性解的混沌性解的分岔性……§5.2Euler方法第五章：常微分方程数值解第一步：连续变量离散化第二步：用直线步进·····1、Euler格式Euler格式例P106第五章：常微分方程数值解初值问题Bernoulli型方程常数变易法第五章：常微分方程数值解第五章：常微分方程数值解令将代入Euler格式计算得x精确值Euler方法Euler方法误差01.00000001.00000000.00000000.11.09544511.10000000.00455490.21.18321601.19181820.00860220.31.26491111.27743780.01252680.41.34164081.35821260.01657180.51.41421361.43513290.02091940.61.48323971.50896630.02572660.71.54919331.58033820.03114490.81.61245151.64978340.03733190.91.67332011.71777930.044459311.73205081.78477080.0527200第五章：常微分方程数值解Euler值2、Euler格式的误差分析第五章：常微分方程数值解事实上Euler格式的每一步都存在误差，为了方便讨论算法的好坏，假定第n步准确的前提下分析第n+1步的误差，称为局部截断误差。定义1定义2第五章：常微分方程数值解Euler格式的误差即Euler格式具有一阶精度将在点Taylor展开的计算格式第五章：常微分方程数值解3、Euler格式与后退Euler格式Euler格式的值后退Euler格式的值隐式格式第五章：常微分方程数值解4、后退格式的精度局部误差分析的要求第五章：常微分方程数值解对照即后退Euler格式具有1阶精度第五章：常微分方程数值解3、梯形格式得梯形格式Euler格式的值后退Euler格式的值第五章：常微分方程数值解梯形格式几何解释第五章：常微分方程数值解4、梯形格式的精度局部误差分析的要求第五章：常微分方程数值解梯形格式5、改进的Euler格式预测校正第五章：常微分方程数值解隐式格式为方便计算，一般用以下改进格式计算用改进格式计算例5.1的结果为x精确值Euler方法改进Euler方法Euler方法误差改进Euler误差01.00000001.00000001.00000000.00000000.00000000.11.09544511.10000001.09590910.00455490.00046400.21.18321601.19181821.18409660.00860220.00088060.31.26491111.27743781.26620140.01252680.00129030.41.34164081.35821261.34336020.01657180.00171940.51.41421361.43513291.41640190.02091940.00218840.61.48323971.50896631.48595560.02572660.00271590.71.54919331.58033821.55251410.03114490.00332080.81.61245151.64978341.61647480.03733190.00402320.91.67332011.71777931.67816640.04445930.004846311.73205081.78477081.73786740.05272000.0058166第五章：常微分方程数值解第五章：常微分方程数值解6、预测精度的改进—两步Euler格式第五章：常微分方程数值解如果令第五章：常微分方程数值解整理得两步Euler格式具有2阶精度。因此可得具有2阶精度的Euler改进格式预测校正第五章：常微分方程数值解问题：精度还能提高吗？§5.3Lunge-Kutta方法1、二阶Lunge-Kutta方法（P113-P115)第五章：常微分方程数值解依据精度要求的待定系数法令加权平均斜率

在点的斜率

在点的斜率用Euler格式预测第五章：常微分方程数值解对照第五章：常微分方程数值解可解得得改进的Euler格式3、三阶Lunge-Kutta方法第五章：常微分方程数值解补充第五章：常微分方程数值解分别将第五章：常微分方程数值解由得一个三阶精度的Runge-Kutta格式4、四阶Lunge-Kutta方法第五章：常微分方程数值解补充点第五章：常微分方程数值解同理可得一个经典四阶精度的Runge-Kutta格式§5.4几种方法的数值计算例5.1P106改进的Euler格式Euler格式第五章：常微分方程数值解四阶经典Lunge-Kutta方法例5.1P106第五章：常微分方程数值解x精确值Euler方法改进Euler方法四阶Lunge-KuttaEuler方法误差改进Euler误差四阶Lunge-Kutta误差01.00000001.00000001.00000001.00000000.00000000.00000000.00000000.11.09544511.10000001.09590911.09544550.00455490.00046400.00000040.21.18321601.19181821.18409661.18321670.00860220.00088060.00000080.31.26491111.27743781.26620141.26491220.01252680.00129030.00000120.41.34164081.35821261.34336021.34164240.01657180.00171940.00000160.51.41421361.43513291.41640191.41421560.02091940.00218840.00000200.61.48323971.50896631.48595561.48324220.02572660.00271590.00000250.71.54919331.58033821.55251411.54919650.03114490.00332080.00000310.81.61245151.64978341.61647481.61245530.03733190.00402320.00000380.91.67332011.71777931.67816641.67332470.04445930.00484630.000004611.73205081.78477081.73786741.73205640.05272000.00581660.0000056几种方法的结果与误差第五章：常微分方程数值解参考程序-Eulerx=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10y0(k+1)=y0(k)+h*(y0(k)-2*x0(k)/y0(k));endfori=1:10y1(1)=1.0;y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i));y1(i+1)=y1(i)+h*(y1(i)-2*x0(i)/y1(i)+y1(i+1)-2*x0(i+1)/y1(i+1))/2;endplot(x,y,'b');holdon;plot(x0,y0,'or');holdon;plot(x0,y1,'*');第五章：常微分方程数值解参考程序-Lunge_Kuttax=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;fork=1:10k1=y0(k)-2*x0(k)/y0(k);k2=y0(k)+h*k1/2-(2*x0(k)+h)/(y0(k)+h*k1/2);k3=y0(k)+h*k2/2-(2*x0(k)+h)/(y0(k)+h*k2/2);k4=y0(k)+h*k3-2*(x0(k)+h)/(y0(k)+h*k3);y0(k+1)=y0(k)+h*(k1+2*k2+2*k3+k4)/6;endholdon;plot(x,y,'b');holdon;plot(x0,y0,'or');第五章：常微分方程数值解§5.5线性多步方法1、Adams显式格式第五章：常微分方程数值解Euler公式梯形公式····参看P27(2.5.12)第五章：常微分方程数值解第五章：常微分方程数值解第五章：常微分方程数值解r=3时第五章：常微分方程数值解当r=1时，（5.5.6）式的误差分析第五章：常微分方程数值解第五章：常微分方程数值解例5.1P106第五章：常微分方程数值解Adams程序x=0:0.01:1;y=sqrt(1+2.*x);a=0.0;b=1.0;n=10;h=(b-a)/n;x0=a:h:b;y0(1)=1.0;y0(2)=1.0954;y0(3)=1.1832;y0(4)=1.2649;fork=4:10y0(k+1)=y0(k)+h*(55*(y0(k)-2*x0(k)/y0(k))-59*(y0(k-1)-2*x0(k-1)/y0(k-1))+37*(y0(k-2)-2*x0(k-2)/y0(k-2))-9*(y0(k-3)-2*x0(k-3)/y0(k-3)))/24;endholdon;plot(x,y,'b');holdon;plot(x0,y0,'or');第五章：常微分方程数值解x精确值四阶Lunge-KuttaAdams四阶Lunge-Kutta误差Adams误差01.00000001.00000001.00000000.00000000.00000000.11.09544511.09544551.09540000.00000040.00004510.21.18321601.18321671.18320000.00000080.00001600.31.26491111.26491221.26490000.00000120.00001110.41.34164081.34164241.34153280.00000160.00010800.51.41421361.41421561.41402400.00000200.0001896