常微分方程数值解法PPT

偏微分方程数值解法偏微分方程数值解法数理学院数学教研室数理学院数学教研室北京北京中国地质大学中国地质大学China University of Geosciences，Beijing教材教材偏微分方程数值解法偏微分方程数值解法（第二版（第二版）清华大学出版社清华大学出版社 陆金甫陆金甫 关治关治 著著3参考资料参考资料微分方程数值方法微分方程数值方法 (第二版第二版)，胡健伟胡健伟,汤怀民著汤怀民著,科学出版社科学出版社,2007,2参考资料参考资料偏微分方程数值解法偏微分方程数值解法（第二版）（第二版）高等教育出版社高等教育出版社李荣华李荣华 著著偏微分方程数值解法偏微分方程数值解法（第二版）（第二版）科学出版社科学出版社孙志忠孙志忠 著著偏微分方程数值解讲义偏微分方程数值解讲义北京大学出版社北京大学出版社李治平李治平 著著 学习资料信箱学习资料信箱:密码密码:numnum任课教师任课教师:李明霞李明霞考核方式：考核方式：作业作业+期末考试期末考试科学理论科学理论科学实验科学实验科学计算科学计算科学方法科学方法科学计算科学计算自然科学自然科学技术与工程科学技术与工程科学PDE求解求解PDE数值解的应用数值解的应用 挪威气象学家挪威气象学家V.Bjerknes（1904）提出数值预报的思想：通过）提出数值预报的思想：通过求解一组方程的初值问题可以预报将来某个时刻的天气的思想；求解一组方程的初值问题可以预报将来某个时刻的天气的思想；L.F.Richardson(1922)：开创了利用数值积分进行预报天气的：开创了利用数值积分进行预报天气的先例，由于一些原因（如，计算稳定性问题先例，由于一些原因（如，计算稳定性问题“Courant,1928”）并没有取得预期的效果并没有取得预期的效果尝试；尝试；Charney,Fjortoft,and Von Neumann(1950),借助于借助于Princeton大学的的计算机大学的的计算机(ENIAC)，利用一个简单的正压涡度，利用一个简单的正压涡度方程（方程（C.G.Rossby,1940）对天气形式作了）对天气形式作了24小时预报小时预报-成功；成功；1.数值天气预报数值天气预报The Electronic Numerical Integrator and Computer(ENIAC).2.核试验核试验:仪器无法测量变化过程，仪器无法测量变化过程，复杂非线性偏微，无法精确求解；复杂非线性偏微，无法精确求解；数值核试验数值核试验:减少核试验次数，节约经费，减少核试验次数，节约经费，缩短研制周期缩短研制周期.3.风洞实验：设备与实验花费昂贵风洞实验：设备与实验花费昂贵;数值风洞数值风洞:周期短，费用低，容易改变参数周期短，费用低，容易改变参数.4.战争决策：战争决策：海湾战争海湾战争(Navier Stokes方程组方程组)PDE数值解的应用数值解的应用主要内容主要内容常微分方程数值解法：常微分方程数值解法：有限差分方法有限差分方法有限元方法有限元方法有限体积法有限体积法双曲型方程有限差分方法双曲型方程有限差分方法抛物型方程有限差分方法抛物型方程有限差分方法椭圆型方程有限差分方法椭圆型方程有限差分方法2.偏微分方程数值解法：偏微分方程数值解法：单步法，多步法单步法，多步法常微分方程数值解常微分方程数值解数值求解初探数值求解初探常微分方程常微分方程偏微分方程偏微分方程:未知函数是一元函数未知函数是一元函数 ODEODE分类分类:未知函数是多元函数未知函数是多元函数又称数学物理方程又称数学物理方程,PDE,PDE常微分方程的数值解常微分方程的数值解 1963年，美国气象学家年，美国气象学家Lorenz在在研究研究热对流的不稳定问题热对流的不稳定问题时，使时，使用高截断的谱方法，由用高截断的谱方法，由Boussinesq流体的闭合方程组得流体的闭合方程组得到了一个完全确定的三阶常微分到了一个完全确定的三阶常微分方程组，即著名的方程组，即著名的Lorenz系统。系统。例1：自变函数 function xdot=lorenzeq(t,x)xdot=-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(1)*x(2)+28*x(2)-x(3);t_final=100;x0=0;0;1e-10;%t_final为设定的仿真终止时间 t,x=ode45(lorenzeq,0,t_final,x0);plot(t,x),figure;%打开新图形窗口 plot3(x(:,1),x(:,2),x(:,3);axis(10 42-20 20-20 25);%根据实际数值手动设置坐标系 可采用comet3()函数绘制动画式的轨迹。comet3(x(:,1),x(:,2),x(:,3)例2：描述函数：function dx=apolloeq(t,x)mu=1/82.45;mu1=1-mu;r1=sqrt(x(1)+mu)2+x(3)2);r2=sqrt(x(1)-mu1)2+x(3)2);dx=x(2);2*x(4)+x(1)-mu1*(x(1)+mu)/r13-mu*(x(1)-mu1)/r23;x(4);-2*x(2)+x(3)-mu1*x(3)/r13-mu*x(3)/r23;求解：x0=1.2;0;0;-1.04935751;tic,t,y=ode45(apolloeq,0,20,x0);tocelapsed_time=0.8310 length(t),plot(y(:,1),y(:,3)ans=689得出的轨道不正确，默认精度RelTol设置得太大，从而导致的误差传递，可减小该值。改变精度：options=odeset;options.RelTol=1e-6;tic,t1,y1=ode45(apolloeq,0,20,x0,options);tocelapsed_time=0.8110 length(t1),plot(y1(:,1),y1(:,3),ans=1873欧拉法欧拉法折线法折线法1.常微分方程能直接进行积分的是少数，而多数是常微分方程能直接进行积分的是少数，而多数是 借助于计算机来求常微分方程的近似解；借助于计算机来求常微分方程的近似解；2.有限差分法是常微分数值解法中有效的方法；有限差分法是常微分数值解法中有效的方法；3.建立差分算法的两个基本的步骤：建立差分算法的两个基本的步骤：a.对解的存在域剖分；对解的存在域剖分；b.采用不同的算法可得到对微分方程不同的逼近采用不同的算法可得到对微分方程不同的逼近 局部截断误差（相容性）局部截断误差（相容性）；c.数值解对真解的精度数值解对真解的精度整体截断误差（收敛性）整体截断误差（收敛性）；d.数值解收敛于真解的速度数值解收敛于真解的速度（收敛速度）（收敛速度）；e.差分格式的计算差分格式的计算舍入误差（稳定性）舍入误差（稳定性）.将微分方程通过差分方程转化为代数方程解。将微分方程通过差分方程转化为代数方程解。（误差）（误差）在常微分方程差分法中最简单的方法是在常微分方程差分法中最简单的方法是Euler方方法，尽管在计算中不会使用，但从中可领悟到建法，尽管在计算中不会使用，但从中可领悟到建立差分格式的技术路线，下面将对其作详细介绍立差分格式的技术路线，下面将对其作详细介绍:差分方法的基本思想就是差分方法的基本思想就是“以差商代替微商以差商代替微商”23()1111()()()()()()()2!3!nnnu thu tu t hu t hut hut hO hn 23()1111()()()()()()()2!3!nnnu thu tu t hu t hut hut hO hn 考虑如下两个考虑如下两个Taylor公式：公式：（1）（2）从（从（1）得到：）得到：1()()()()iiiu tu tu tO hh 从（从（2）得到：）得到：1()()()()iiiu tu tu tO hh 211()()()()2iiiu tu tu tO hh 从（从（1）减（）减（2）得到：）得到：从（从（1）+（2）得到：）得到：2112()2()()()()iiiiu tu tu tu tO hh 23()1111()()()()()()()2!3!nnnut hutu t hu t hu t hut hOhn 23()1111()()()()()()()2!3!nnnu t hu tu t hu t hu t hut hO hn （1）（2）27()()()lim0nnntu ttu tu tt()()lim0nnhu thu th()()nnu thu th()()()()()232nnnnhu thu thu tutO h()()()nnu thu tO hh由由Taylor展开式展开式总结：总结：()()()()u thu tu tO hh28数值微分公式数值微分公式()()()u tu thO hh()()()22u thu thO hh向前差分向前差分向后差分向后差分中心差分中心差分()()()()()232nnnnhu thu thu tutO h()()()()()232nnnnhu thu thu tutO h()()()()u thu tu tO hh29数值微分公式数值微分公式()()()u tu thO hh()()()22u thu thO hh向前差商向前差商向后差商向后差商中心差商中心差商()()()()()232nnnnhu thu thu tutO h()()()()()232nnnnhu thu thu tutO h对经典的初值问题对经典的初值问题0(,)(0)duf t udtuu (0,)tT 1212(,)(,)f t uf t uL uu 满足满足Lipschitz条件条件保证了方程组的初值问题有唯一解。保证了方程组的初值问题有唯一解。算法构造算法构造：0tuTnt0tit1iitth 1.在求解域上等距离分割：在求解域上等距离分割：2.在在 有：有：1,iit t 1()(,()()iiiiu tfut u thtO h 1(,)iiiiuf thuu 1(,)iiiiuuhf t u 微微分方分方程的精程的精确解确解差差分方分方程的精程的精确解确解(,)duf t udt 3.应用时采用如下递推方式计算：应用时采用如下递推方式计算：0u1u2u3u0t1t2t3t0(0)uu 1(,)iiiiuuhf t u 3311(,)f t u斜斜率率Euler法几何意义及误差法几何意义及误差()(,()u tf t u t 有差分格式(,)1nnnnuuhf t utu0t0u1t2t Nt()u t1()u t2()u t()Nu t11(,)t u(,)00f t u斜斜率率22(,)t u34例例1 利用利用Euler方法计算初值问题方法计算初值问题 22100,(0)0utuu初初值值问问题题的解在的解在t=0.3处的数值解处的数值解.步长步长h=0.1解解:Euler公式为公式为:221(100)nnnnuuh tu 00.1,0hu 由由计计算算得得22100(0)1000uuhhu 22211(1)1000.001uuhhu 22322(2)1000.00501uuhhu 4.例子例子例例2对初值问题对初值问题2(0)1yxyy (0,1)x 5,N 0.2h 用用Euler法求解，用法求解，用即，即，001112223334445(,)1.000,1.200;(,)1.600,1.520;(,)2.320,1.984;(,)3.184,2.621;(,)4.221,3.465f xyyf xyyf xyyf xyyf xyy36例例3 利用利用Euler方法求数值解方法求数值解 1,(0)12uuu初初值值问问题题步长步长h=0.1,解区间解区间0,1 绘制折线，与真解比较绘制折线，与真解比较 37Matlab实现实现h=0.1;u(1)=1;for n=1:10 u(n+1)=u(n)+h*0.5*u(n);endt=0:0.1:1;plot(t,u,ro,Linewidth,2)ut=exp(0.5*t);hold onplot(t,ut,Linewidth,2)380.00.10.20.30.40.50.60.70.80.91精确解精确解ut数值解数值解un 节点节点 ti 1.0000 1.0500 1.1025 1.1576 1.2155 1.2763 1.3401 1.4071 1.4775 1.5513 1.6289 1.0000 1.0513 1.1052 1.1618 1.2214 1.2840 1.3499 1.4191 1.4918 1.5683 1.64872(0)1dyxydxyy (0,1)x 其解析解为：其解析解为：12yx 例例412()nnnnnxyyh yy h=0.2;u(1)=1;x=0:0.2:1;for n=1:5 u(n+1)=u(n)+h*(u(n)-2*x(n)/u(n);endplot(x,u,-ro,Linewidth,2)hold onut=sqrt(1+2*x);plot(x,ut,Linewidth,2)00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.81.9xy42Euler方法的三种解释方法的三种解释1.数值微分：数值微分：用差商来代替导数用差商来代替导数2.数值积分：数值积分：把微分方程变成积分方程把微分方程变成积分方程3.幂级数展开：幂级数展开：将将u(t+h)在在t 做做Taylor展开展开()()(,()u thu thf t u t ()()(,()(,()t htu thu tfudhf t u t 2()()()().2!()(,()hu thu thututu thft u t 一一、局部阶段误差局部阶段误差-相容性相容性0t0tit1iitth在递推的每一步，设定在递推的每一步，设定()iiu tu1()iu t1iu2()O h()iu t(),iiu tu 过点过点(,)iit u作作 的切线，该的切线，该切线的切线的()u t方程为方程为：11()()iiiiiuuu ttt 即：即：1(,)iiiiuuhf t u 方法分析方法分析:44局部截断误差局部截断误差：假设第：假设第i步精确计算的前步精确计算的前提下，数值解提下，数值解 和精确解和精确解 的误差的误差1iR 1iu 1()iu t(,),.,()1000 11iiiiuuht u hiNuu t1111()()()(,(),)iiiiiiiRu tuu tu tht u th ()iiuu t 111222()()()(,()()()(,()()()()()(,()2!()(,()()()(,()2!()2!iiiiiiiiiiiiiiiiiiiiiu tuu tu thf t u tu thu thf t u thu thu tuu thf t u thu thf t u tuu thf t u thu 1(,)iiiiuuhf t u 21()iRO h Euler法法相容性和相容的阶相容性和相容的阶针对的是建立差分格式时由差商代针对的是建立差分格式时由差商代替微商所引起的替微商所引起的局部截断误差局部截断误差.11()(1)qiRO hq Euler法1阶相容21()()(,()()iiiiiRu thu thf t u tO h Euler法：法：1(,)iiiiuuhf t u :若一个离散变量方法的局部截若一个离散变量方法的局部截断误差对任意断误差对任意i满足：满足：整体截断误差是以点整体截断误差是以点0t 的初始值的初始值0u 为出发值，用数值方法推进为出发值，用数值方法推进i+1步到点步到点1it，所得的近似值，所得的近似值 与精确值与精确值 的偏差：的偏差：二二.整体截断误差整体截断误差收敛性收敛性1iu 1()iu t 111()iiiu tu 称为整体截断误差。称为整体截断误差。1(,)iiiiuuhf t u 11()()(,()iiiiiu tu thf t u tR21(,()(,)iiiiiih f t u tf t uCh 21iiihLCh 2210(1)1(1)(1)(1)iiihLChhLhLhL 20(1)(1)1iiChhLhLhL Lipchitz条件0(1)TLTLCheeL 特例，特例，00 若不计初始误差，即若不计初始误差，即则则1(1)TLiCheL 1()iO h 即即欧拉法欧拉法1阶收敛阶收敛2312!3!1(0)xxxxexxex 注：注：20(1)(1)1iiChhLhLhL (1)iihLTLhLee 49收敛性与收敛的阶收敛性与收敛的阶研究的是误差累积产生的研究的是误差累积产生的整体整体 截断误差截断误差.对任意的对任意的t(t0,T，成立，成立若此时，整体截断误差满足若此时，整体截断误差满足 则称方法的收敛为则称方法的收敛为 p阶的阶的.00lim()hntnhtuu t ()pnO h ()iiiu tu 三三.舍入误差舍入误差稳定性稳定性假设一个计算机仅表示假设一个计算机仅表示4个数字（小数点后面），个数字（小数点后面），那么那么10.3333,3 10.11119 计算计算20.333410.11120.11119110.33340.33333xxx 20.33340.33341190.6667133xxxxx 误差大我们的要求是：我们的要求是：最初产生的小误差最初产生的小误差在以后的计算中虽然会传递下去，在以后的计算中虽然会传递下去，但但不会无限制的扩大不会无限制的扩大，这就是稳定性所描述的问题。，这就是稳定性所描述的问题。下面引进稳定性的概念：下面引进稳定性的概念：tu00u0vntnunv设由初值设由初值0u0v得到精确解得到精确解 ，nu由初值由初值得到精确解得到精确解 ，nv若存在常数若存在常数C和充分小的步长和充分小的步长0h使得使得00nnuvC uv 则称数值方法是稳定的。则称数值方法是稳定的。unu1212(,)(,)f t uf t uL uu 四、改进的四、改进的Euler法法将微分方程将微分方程()(,()u tf t u t 在区间在区间1,iit t 上积分，得到上积分，得到11()()(,()iitiitu tu tf t u t dt 用梯形法计算积分的近似值，有用梯形法计算积分的近似值，有111(,()(,()(,()2iitiiiithf t u t dtf tu tf tu t 于是于是1111()()(,()(,()2iiiiiiu tu th f t u tf tu t 1111(,)(,)2iiiiiiuuh f t uf tu 这是一个隐式格式，一般需要用迭代法来求这是一个隐式格式，一般需要用迭代法来求 ，而用显式的而用显式的Euler法提供初值。法提供初值。1iu(0)1(,)iiiiuuhf t u(k 1)(k)111/2(,)(,)iiiiiiuuhf t uf tu 为了简化计算的过程，在此基础上进一步变为如下算法：为了简化计算的过程，在此基础上进一步变为如下算法：111/2(,)(,)iiiiiiuuhf t uf tu1(,)iiiiuuhf t u 此式称为此式称为“改进的改进的Euler法。法。预估预估校正校正其局部截断误差为其局部截断误差为3()O h只迭代一次只迭代一次计算计算隐式法，一般需多隐式法，一般需多次迭代计算次迭代计算1iu 1iu 接下来讨论其几何意义接下来讨论其几何意义tu0it1it0(,)iimf t u()u tiu1iu111(,)iimf tu 122averagemmm 1iu1()iu t00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.81.9xyEuler法、改进的法、改进的Euler法法和解析解的比较和解析解的比较2(0)1dyxydxyy (0,1)x 12yx h=0.2;u(1)=1;x=0:0.2:1;for n=1:5 u(n+1)=u(n)+h*(u(n)-2*x(n)/u(n);endplot(x,u,-ro,Linewidth,2)hold onut=sqrt(1+2*x);plot(x,ut,Linewidth,2)hold onfor n=1:5 z0=u(n)+h*(u(n)-2*x(n)/u(n);u(n+1)=u(n)+h/2*(u(n)-2*x(n)/u(n)+(z0-2*x(n+1)/z0);endplot(x,u,-bs,Linewidth,2)58总结：基本步骤总结：基本步骤 解差分方程，求出格点函数解差分方程，求出格点函数 对区间作分割：对区间作分割：01N:ItttT 求求 u(x)在在tn 上的近似值上的近似值un。由微分方程出发，建立求格点函数的差分方程。由微分方程出发，建立求格点函数的差分方程。这个方程应该满足：这个方程应该满足：A、解存在唯一；、解存在唯一；B、相容、相容；C、稳定，收敛；、稳定，收敛；目的目的关键关键59为了考察数值方法提供的数值解，是否有实用价值，为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论：需要知道如下几个结论：1.1.差分方程对微分方程的逼近程度如何，差分方程对微分方程的逼近程度如何，即即相容性相容性问题。问题。2.2.步长充分小时，所得到的数值解能否逼近步长充分小时，所得到的数值解能否逼近 问题的真解，逼近程度如何，即问题的真解，逼近程度如何，即收敛性收敛性问题。问题。3.3.产生的舍入误差，在以后得各步计算中，是否会产生的舍入误差，在以后得各步计算中，是否会 无限制扩大，即无限制扩大，即稳定性稳定性问题。问题。60相容性（局部截断误差）相容性（局部截断误差）收敛性（整体截断误差）收敛性（整体截断误差）稳定性（舍入误差）稳定性（舍入误差）数值方法的基本问题数值方法的基本问题61 微分方程微分方程差分方程差分方程真解真解u=u(t)真解真解u=un方法分析方法分析:相相容容性性收收敛敛性性 稳定性稳定性 62数值求解微分方程过程示意数值求解微分方程过程示意微微分分方方程程区域剖分区域剖分离散系统的离散系统的性态研究性态研究递推计算或解线递推计算或解线性代数方程组性代数方程组微分方程离散微分方程离散初始和边界条件处理初始和边界条件处理解的存在性、唯一性解的存在性、唯一性解的收敛性和收敛速度解的收敛性和收敛速度解的稳定性解的稳定性得到数值解得到数值解64课堂练习课堂练习1 利用利用Euler方法求数值解方法求数值解 30,(0)1uuu 初初值值问问题题步长步长h=0.1,解区间解区间0,0.5 绘制折线，与真解比较绘制折线，与真解比较 采用采用h=0.05,0.005，再与真解比较，再与真解比较cleara=0;b=0.5;N=5h=(b-a)/N;t=0:h:b;u(1)=1;for n=1:N u(n+1)=u(n)-h*30*u(n);endplot(t,u,-ro,Linewidth,2)ut=exp(-30*t);hold onplot(t,ut,Linewidth,2)1.试对问题试对问题 建立向后建立向后Euler格式和梯形格式格式和梯形格式,并讨论它们的局部截断误差并讨论它们的局部截断误差.作业作业0(,)(0)duf t udtuu (0,)tT 67 刚才的发言，如刚才的发言，如有不当之处请多指有不当之处请多指正。谢谢大家！正。谢谢大家！