(完整版)偏微分方程数值解

1偏微分方程数值解(NumericalSolutionofPartialDifferentialEquations)主讲：王曰朋[e-mail]2参考数目1.GeorgeJ.Haltiner,RogerTerryWilliams,NumericalPredictionandDynamicMeteorology(2ndEdition),theUnitedStatesofAmerica,1979.2.CurtisF.GeraldandPatrickO.,AppliedNumericalAnalysis,PersonEducation,Inc.,2004.3.EugeniaKalnay,AtmosphericModeling,DataAssimilationandPredictability,thepressSyndicateoftheUniversityofCambridge,2003.4.AriehIserles,AFirstCourseintheNumericalAnalysisofDifferentialEquations,CambridgeUniversityPress,1996.5.李荣华，冯国忱.微分方程数值解.北京：人民教育出版社，1980.6.徐长发，李红.实用偏微分方程数值解法.华中科技大学出版社，2003.7.沈桐立，田永祥等.数值天气预报.北京：气象出版社，2007.3数值天气预报—PDE数值解1.挪威气象学家V.Bjerknes（1904）提出数值预报的思想：通过求解一组方程的初值问题可以预报将来某个时刻的天气—思想；2.L.F.Richardson(1922)：开创了利用数值积分进行预报天气的先例，由于一些原因（如，计算稳定性问题“Courant,1928”）并没有取得预期的效果—尝试；3.Charney,Fjortoft,andVonNeumann(1950),借助于Princeton大学的的计算机(ENIAC)，利用一个简单的正压涡度方程（C.G.Rossby,1940）对500mb的天气形式作了24小时预报---成功；4TheElectronicNumericalIntegratorandComputer(ENIAC).5常微分方程的数值解大气科学中常微分方程和偏微分方程的关系1.大气行星边界层(近地面具有湍流运动特性的大气薄层,1—1.5km),埃克曼（V.W.Ekman）（瑞典）螺线的导出；2.1963年，美国气象学家Lorenz在研究热对流的不稳定问题时，使用高截断的谱方法，由Boussinesq流体的闭合方程组得到了一个完全确定的三阶常微分方程组，即著名的Lorenz系统。6Lorenz系统dx/dt=a(y-x)dy/dt=x(b-z)-ydz/dt=xy-cz其中，a=10,(Prandtlnumber);b=28(Rayleighnumber);c=8/3;(x,y,z)_0=(0.01;0.01;1e-10)705101520253035404550-30-20-10010203040508-20020-30-20-1001020300102030405091011Franceshini将Navier-Stokes方程截断为五维的截谱模型如下：11234522133312441655142449357Re53xxxxxxxxxxxxxxxxxxxxxxì¢=-++ïïïï¢=-+ïïïïï¢=--+íïï¢ï=--ïïï¢=--ïïîï12欧拉法—折线法•常微分方程能直接进行积分的是少数，而多数是借助于计算机来求常微分方程的近似解；•有限差分法是常微分方程中数值解法中通常有效的方法；•建立差分算法的两个基本的步骤：1.建立差分格式，包括：a.对解的存在域剖分；b.采用不同的算法可得到不同的逼近误差—截断误差（相容性）；c.数值解对真解的精度—整体截断误差（收敛性）；d.数值解收敛于真解的速度；e.差分算法—舍人误差（稳定性）.132.差分格式求解将积分方程通过差分方程转化为代数方程求解，一般常用递推算法。在常微分方程差分法中最简单的方法是Euler方法，尽管在计算中不会使用，但从中可领悟到建立差分格式的技术路线，下面将对其作详细介绍:14差分方法的基本思想“就是以差商代替微商”23()1111()()()()()()()2!3!!nnnuthututhuthuthuthOhn23()1111()()()()()()()2!3!!nnnuthututhuthuthuthOhn考虑如下两个Taylor公式：（1）（2）从（1）得到：1()()()()iiiutututOhh15从（2）得到：1()()()()iiiutututOhh211()()()()2iiiutututOhh从（1）-（2）得到：从（1）+（2）得到：2112()2()()()()iiiiututututOhh16对经典的初值问题0(,)(0)duftudtuu(0,)tT1212(,)(,)ftuftuLuu满足Lipschitz条件保证了方程组的初值问题有唯一解。17一、算法构造：0tuTnt0tit1iitth1.在求解域上等距离分割：2.在有：1[,]iitt1()()(,)()iiututftuOhh1(,)iiiiuuftuh1(,)iiiiuuhftu微分方程的精确解差分方程的精确解183.应用时采用如下递推方式计算：0u1u2u3u0t1t2t4.例题对初值问题2(0)1yxyy(0,1)x5,N0.2h用Euler法求解，用即，001111223334445(,)1.000,1.200;(,)1.600,1.520;(,)2.320,1.984;(,)3.184,2.621;(,)4.221,3.465fxyyfxyyfxyyfxyyfxyy195.Euler法的几何意义0t0tit1iitth在递推的每一步，设定()iiutu1()iut1iu2()Oh()iut(),iiutu过点(,)iitu作的切线，该切线的()ut方程为：11()()iiiiiuuuttt即：1(,)iiiiuuhftu20二、误差分析构造算法后，这一算法在实际中是否可行呢？也就是说是否使计算机仿真而不失真，这还需要进一步分析。1.局部截断误差--相容性为了分析分析数值方法的精确度，常常在()iiuut成立的假定下，估计误差111()iiieutu这种误差称为“局部截断误差”，如图。21()()()()2iiihututhutu2()(,())()2iiihuthftutu22111()()()2iiiheutuuOh局部截断误差是以点it的精确解()iut为出发值，用数值方法推进到下一个点1it而产生的误差。21整体截断误差是以点0t的初始值0u为出发值，用数值方法推进i+1步到点1it，所得的近似值与精确值的偏差：2.整体截断误差—收敛性1iu1()iut111()iiiutu称为整体截断误差。1(,)iiiiuuhftu11()()(,())iiiiiututhftute21[(,()(,)]iiiiiihftutftuDh21iiihLDh12210(1)[1(1)(1)(1)]iiihLDhhLhLhL2110(1)[(1)1]iiDhhLhLhL220(1)TLTLDheeL特例，00若不计初始误差，即则1(1)TLiDheL1()iOh即3.舍入误差—稳定性假设一个计算机仅表示4个数字（小数点后面），那么10.3333,310.11119计算20.373410.1112011119110.33340.33333xxx20.33341190.6667133xxxx23我们的要求是：最初产生的小误差在以后的计算中虽然会传递下去，但不会无限制的扩大，这就是稳定性所描述的问题。下面引进稳定性的概念：设由初值0u0v得到精确解，nu由初值得到精确解，nv若存在常数C和充分小的步长0h使得00nnuvCuv则称数值方法是稳定的。tu00u0vntnunv242(0)1dyxydxyy计算例题(0,1)x其解析解为：12yxx=00.20000.40000.60000.80001.0000y=1.00001.20001.37331.53151.68111.82692500.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.81.9xy26三、改进的Euler法将微分方程()(,())utftut在区间1[,]iitt上积分，得到11()()(,())iitiitututftutdt用梯形法计算积分的近似值，有1111(,())[(,())(,())]2iitiiiitftutdtftutftut于是1111()()[(,())(,())]2iiiiiiututhftutftut1111[(,)(,)]2iiiiiiuuhftuftu这是一个隐式格式，一般需要用迭代法来求，而用显式的Euler法提供初值。1iu27[0]1(,)iiiiuuhftu[1][]111/2[(,)(,)]kkiiiiiiuuhftuftu为了简化计算的过程，在此基础上进一步变为如下算法：111/2[(,)(,)]iiiiiiuuhftuftu1(,)iiiiuuhftu此式称为“改进的Euler法。接下来讨论其几何意义预估校正其局部截断误差为3()Oh这个问题将在下节讨论。28tu0it1it0(,)iimftu()utiu1iu111(,)iimftu122averagemmm1iu1()iut2900.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30四、（龙格-库塔）Runge-Kutta方法简单的Euler法是建立在Taylor级数的一项展开；改进的Euler法是以两项Taylor级数为基础建立的，如：23()()()()()2iiiihuthuthututOh23()()()[(()())/]/2()iiiiiuthuthuthuthuthOh()()()2iiiuthututh11(,)(,)2iiiiiftuftuuh如果我们截取Taylor级数的更多项会得到什么样的求解方法呢？两个德国数学家（C.Runge&M.kutta）以这种思想为基础建立了求解微分方程的龙格-库塔方法。它是常微分方程数值解法中使用最为广泛的方法之一。31一般地，一个K阶的Runge-Kutta方法可用下面的公式表示：11111(,)(,)kiijjjiijjijijllluuwKKhftuKhftchuaK2,3,,jk其中，jw是待定的加权系数，2321,1,,,,,,kkkcccaa是待定的系数。Euler法就是11,1kw的R-K法。其系数的确定如下：将1iu展开成h的幂级数，并与微分方程的精确解1()iut在点it的Taylor展开式相比较，使两者的前1p项相同，这样确定的R-K法，其