常微分方程数值解法.ppt

偏微分方程数值解法,数理学院数学教研室,北京中国地质大学ChinaUniversityofGeosciences，Beijing,教材,偏微分方程数值解法（第二版）清华大学出版社陆金甫关治著,3,参考资料,微分方程数值方法(第二版)，胡健伟,汤怀民著,科学出版社,2007,2,参考资料,偏微分方程数值解法（第二版）高等教育出版社李荣华著,偏微分方程数值解法（第二版）科学出版社孙志忠著,偏微分方程数值解讲义北京大学出版社李治平著,学习资料信箱:pdenum密码:numnum任课教师:李明霞mathlmx考核方式：作业+期末考试,科学理论,科学实验,科学计算,科学方法,科学计算,自然科学,技术与工程科学,PDE求解,PDE数值解的应用,挪威气象学家V.Bjerknes（1904）提出数值预报的思想：通过求解一组方程的初值问题可以预报将来某个时刻的天气的思想；L.F.Richardson(1922)：开创了利用数值积分进行预报天气的先例，由于一些原因（如，计算稳定性问题“Courant,1928”）并没有取得预期的效果尝试；Charney,Fjortoft,andVonNeumann(1950),借助于Princeton大学的的计算机(ENIAC)，利用一个简单的正压涡度方程（C.G.Rossby,1940）对天气形式作了24小时预报-成功；,1.数值天气预报,TheElectronicNumericalIntegratorandComputer(ENIAC).,2.核试验:仪器无法测量变化过程，复杂非线性偏微，无法精确求解；数值核试验:减少核试验次数，节约经费，缩短研制周期.,3.风洞实验：设备与实验花费昂贵;数值风洞:周期短，费用低，容易改变参数.,4.战争决策：海湾战争(NavierStokes方程组),PDE数值解的应用,主要内容,常微分方程数值解法：,有限差分方法,有限元方法,有限体积法,双曲型方程有限差分方法,抛物型方程有限差分方法,椭圆型方程有限差分方法,2.偏微分方程数值解法：,单步法，多步法,常微分方程数值解,数值求解初探,常微分方程,偏微分方程,:未知函数是一元函数ODE,分类,:未知函数是多元函数,又称数学物理方程,PDE,常微分方程的数值解,1963年，美国气象学家Lorenz在研究热对流的不稳定问题时，使用高截断的谱方法，由Boussinesq流体的闭合方程组得到了一个完全确定的三阶常微分方程组，即著名的Lorenz系统。,例1：自变函数functionxdot=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(1042-2020-2025);%根据实际数值手动设置坐标系,可采用comet3()函数绘制动画式的轨迹。comet3(x(:,1),x(:,2),x(:,3),例2：,描述函数：functiondx=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.8310length(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.8110length(t1),plot(y1(:,1),y1(:,3),ans=1873,欧拉法折线法,1.常微分方程能直接进行积分的是少数，而多数是借助于计算机来求常微分方程的近似解；2.有限差分法是常微分数值解法中有效的方法；3.建立差分算法的两个基本的步骤：1)建立差分格式，包括：a.对解的存在域剖分；b.采用不同的算法可得到对微分方程不同的逼近局部截断误差（相容性）；c.数值解对真解的精度整体截断误差（收敛性）；d.数值解收敛于真解的速度（收敛速度）；e.差分格式的计算舍入误差（稳定性）.2)差分格式求解将微分方程通过差分方程转化为代数方程解。（误差）,在常微分方程差分法中最简单的方法是Euler方法，尽管在计算中不会使用，但从中可领悟到建立差分格式的技术路线，下面将对其作详细介绍:,差分方法的基本思想就是“以差商代替微商”,考虑如下两个Taylor公式：,（1）,（2）,从（1）得到：,从（2）得到：,从（1）减（2）得到：,从（1）+（2）得到：,（1）,（2）,27,由Taylor展开式,总结：,28,数值微分公式,向前差分,向后差分,中心差分,29,数值微分公式,向前差商,向后差商,中心差商,对经典的初值问题,满足Lipschitz条件,保证了方程组的初值问题有唯一解。,算法构造：,0,t,u,T,1.在求解域上等距离分割：,2.在有：,微分方程的精确解,差分方程的精确解,3.应用时采用如下递推方式计算：,33,Euler法几何意义及误差,34,例1,利用Euler方法计算初值问题,的解在t=0.3处的数值解.步长h=0.1,解:Euler公式为:,4.例子,例2,对初值问题,用Euler法求解，用,即，,36,例3,利用Euler方法求数值解,步长h=0.1,解区间0,1,绘制折线，与真解比较,37,Matlab实现h=0.1;u(1)=1;forn=1:10u(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);holdonplot(t,ut,Linewidth,2),38,0.00.10.20.30.40.50.60.70.80.91,精确解ut,数值解un,节点ti,1.00001.05001.10251.15761.21551.27631.34011.40711.47751.55131.6289,1.00001.05131.10521.16181.22141.28401.34991.41911.49181.56831.6487,其解析解为：,例4,h=0.2;u(1)=1;x=0:0.2:1;forn=1:5u(n+1)=u(n)+h*(u(n)-2*x(n)/u(n);endplot(x,u,-ro,Linewidth,2)holdonut=sqrt(1+2*x);plot(x,ut,Linewidth,2),42,Euler方法的三种解释,1.数值微分：用差商来代替导数2.数值积分：把微分方程变成积分方程3.幂级数展开：将u(t+h)在t做Taylor展开,一、局部阶段误差-相容性,0,t,在递推的每一步，设定,过点,作的切线，该切线的,方程为：,即：,方法分析:,44,局部截断误差,局部截断误差：假设第i步精确计算的前提下，数值解和精确解的误差,Euler法,相容性和相容的阶,相容性针对的是建立差分格式时由差商代替微商所引起的局部截断误差.,Euler法1阶相容,Euler法：,q阶相容:若一个离散变量方法的局部截断误差对任意i满足：,整体截断误差是以点,的初始值,为出发值，用数值方法推进i+1步到点,，所得的近似值与精确值的偏差：,二.整体截断误差收敛性,称为整体截断误差。,Lipchitz条件,特例，,若不计初始误差，即,则,即,欧拉法1阶收敛,注：,49,收敛性与收敛的阶,收敛性:研究的是误差累积产生的整体截断误差.收敛:对任意的t(t0,T，成立收敛阶:若此时，整体截断误差满足则称方法的收敛为p阶的.,三.舍入误差稳定性,假设一个计算机仅表示4个数字（小数点后面），,那么,计算,误差大,我们的要求是：最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大，这就是稳定性所描述的问题。下面引进稳定性的概念：,t,u,0,设由初值,得到精确解，,由初值,得到精确解，,若存在常数,和充分小的步长,使得,则称数值方法是稳定的。,u,四、改进的Euler法,将微分方程,在区间,上积分，得到,用梯形法计算积分的近似值，有,于是,这是一个隐式格式，一般需要用迭代法来求，而用显式的Euler法提供初值。,为了简化计算的过程，在此基础上进一步变为如下算法：,此式称为“改进的Euler法。,预估,校正,其局部截断误差为,只迭代一次计算,隐式法，一般需多次迭代计算,接下来讨论其几何意义,t,u,0,Euler法、改进的Euler法和解析解的比较,h=0.2;u(1)=1;x=0:0.2:1;forn=1:5u(n+1)=u(n)+h*(u(n)-2*x(n)/u(n);endplot(x,u,-ro,Linewidth,2)holdonut=sqrt(1+2*x);plot(x,ut,Linewidth,2)holdonforn=1:5z0=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,总结：基本步骤,解差分方程，求出格点函数,对区间作分割：,求u(x)在tn上的近似值un。,由微分方程出发，建立求格点函数的差分方程。这个方程应该满足：,A、解存在唯一；B、相容；C、稳定，收敛；,目的,关键,59,为了考察数值方法提供的数值解，是否有实用价值，需要知道如下几个结论：,1.差分方程对微分方程的逼近程度如何，即相容性问题。2.步长充分小时，所得到的数值解能否逼近问题的真解，逼近程度如何，即收敛性问题。3.产生的舍入误差，在以后得各步计算中，是否会无限制扩大，即稳定性问题。,60,相容性（局部截断误差）收敛性（整体截断误差）稳定性（舍入误差）,数值方法的基本问题,61,微分方程,差分方程,真解u=u(t),真解u=un,方法分析:,相容性,收敛性,稳定性,62,数值求解微分方程过程示意,微分方程,区域剖分,离散系统的性态研究,递推计算或解线性代数方程组,微分方程离散,初始和边界条件处理,解的存在性、唯一性,解的收敛性和收敛速度,解的稳定性,得到数值解,任何模拟方法，都必须在最佳计算速度和数值精度之间寻找平衡点。,要在各种可能的求解方法中找到一种统一地适用于计算材料学领域（或其它领域）的理想方法，一般是不现实的。,由于实际问题的具体特征、复杂性以及算法自身的适用范围决定了应用中必须选择、设计适合于自己特定问题的算法，因而掌握数值方法的思想至关重要。,64,课堂练习