数值分析》第五讲：常微分方程数值解.ppt

常微分方程数值解,理学院 陈丽娟,数值分析第五讲,第五章：常微分方程数值解,5.1 引言,1、常微分方程与解,为n阶常微分方程。,如果函数 在区间a,b内n阶可导，称方程,为方程满足定解条件的解。,第五章：常微分方程数值解,解的图示,第五章：常微分方程数值解,本教材重点讨论定解问题(初值问题）,定解条件（初始条件）,是否能够找到定解问题的解取决于,仅有极少数的方程可以通过“常数变易法”、“可分离变量法”等特殊方法求得初等函数形式的解，绝大部分方程至今无法理论求解。,2、数值解的思想,第五章：常微分方程数值解,（1）将连续变量 离散为,（2）用代数的方法求出解函数 在 点的近似值,如果找不到解函数 数学界还关注： 解的存在性 解的唯一性 解的光滑性 解的振动性 解的周期性 解的稳定性 解的混沌性 ,5.2 Euler方法,第五章：常微分方程数值解,第一步：连续变量离散化,第二步：用直线步进,1、Euler格式,18世纪最杰出的数学家之一，13岁时入读巴塞尔大学，15岁大学毕业，16岁获得硕士学位。 1727年-1741年（20岁-34岁）在彼得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、造船业等实际问题。 24岁晋升物理学教授。 1735年（28岁）右眼失明。,第五章：常微分方程数值解,1741年 - 1766（34岁-59岁）任德国科学院物理数学所所长，任职25年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。 1766年应沙皇礼聘重回彼得堡，在1771年（64岁）左眼失明。 Euler是数学史上最多产的数学家，平均以每年800页的速度写出创造性论文。他去世后，人们用35年整理出他的研究成果74卷。,第五章：常微分方程数值解,例 P106,第五章：常微分方程数值解,初值问题,Bernoulli方程,由Bernoulli方程的求解方法可得解析解,Euler格式为,步进计算结果见P106表5.1,第五章：常微分方程数值解,Euler值,Euler格式的误差分析,事实上Euler格式的每一步都存在误差，为了方便讨论算法的好坏，假定第 n 步准确的前提下分析第n+1步的误差，称为局部截断误差。,第五章：常微分方程数值解,由Taylor公式,第五章：常微分方程数值解,Euler格式的误差为,2、后退Euler格式,令,得,令,Euler格式,第五章：常微分方程数值解,后退Euler格式,令,得,令,后退Euler格式的值,Euler格式的值,3、梯形格式,4、改进的Euler格式,第五章：常微分方程数值解,为方便计算，一般用以下改进格式计算,用改进格式计算例5.1的结果见P110表5.2,第五章：常微分方程数值解,5、两步Euler格式,第五章：常微分方程数值解,已知Euler格式,即Euler格式具有一阶精度,如果令,并假定,第五章：常微分方程数值解,则有,记,则,其中,第五章：常微分方程数值解,特别要注意的是：一般,两点预测格式具有二阶精度。,当 时,所以,第五章：常微分方程数值解,考察两点校正格式的精度,则,第五章：常微分方程数值解,比较,得,即梯形格式具有二阶精度，因此两步格式从预测到校正均达到二阶精度。,因此得具有二阶精度的两步Euler格式,预测,校正,第五章：常微分方程数值解,5.3 Lunge-Kutta方法,1、二阶Lunge-Kutta方法（P113-P115),第五章：常微分方程数值解,依据精度要求的待定系数法,令,确定 使,具有二阶精度,第五章：常微分方程数值解,对照,第五章：常微分方程数值解,可解得,得,改进的Euler格式,3、三阶Lunge-Kutta方法,第五章：常微分方程数值解,补充,确定参数 使,具有二阶精度,第五章：常微分方程数值解,使,具有三阶精度,即,分别将,在点 Taylor展开，代入（P116),与 的Taylor展开比较,第五章：常微分方程数值解,得,可解得,第五章：常微分方程数值解,得一个三阶精度的Runge-Kutta格式,4、四阶Lunge-Kutta方法见P117,5.4 几种方法的数值计算,例5.1 P106,第五章：常微分方程数值解,四阶经典Lunge-Kutta方法,例5.1 P106,第五章：常微分方程数值解,几种方法的结果与误差,第五章：常微分方程数值解,参考程序-Euler,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; for k=1:10 y0(k+1)=y0(k)+h*(y0(k)-2*x0(k)/y0(k); end for i=1:10 y1(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; end plot(x,y,b); hold on; plot(x0,y0,or); hold on; plot(x0,y1,*);,第五章：常微分方程数值解,参考程序-Lunge_Kutta,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; for k=1:10 k1=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; end hold on; plot(x,y,b); hold on; plot(x0,y0,or);,第五章：常微分方程数值解,5.5 线性多步方法,1、Adams显式格式,第五章：常微分方程数值解,P27 (2.5.12),第五章：常微分方程数值解,第五章：常微分方程数值解,第五章：常微分方程数值解,r=3 时见P125（5.5.6）,第五章：常微分方程数值解,（5.5.6）式的误差分析,第五章：常微分方程数值解,第五章：常微分方程数值解,例5.1 P106,第五章：常微分方程数值解,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; for k=4:10 y0(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; end hold on; plot(x,y,b); hold on; plot(x0,y0,or);,第五章：常微分方程数值解,几种方法的结果与误差,第五章：常微分方程数值解,2、Adams隐式格式,隐式格式见P125,第五章：常微