第十章.常微分方程数值解法

第十章常微分方程初值问题的数值解法最为普遍的模型，最为复杂的问题建模框架：明确地叙述问题→分析其中的物理、化学、生物定律→寻找数量关系：变化率=输入率-输出率求数值解及定性分析建立微分方程世界本质上是非线性的：混沌与现实一阶常系数微分方程的初值问题——基本模型1)(),(00yxyyxfdxdy微分方程的解形成曲线族，或解流形，初值条件定出其中一条曲线。增长常数人口模型例rrNdtdNNN(t)No0tN=N0exp(rt)已知每一点的导数及起始点，重构这条曲线—数值解数值解的基本方法,,,上存在且唯一。在一定条件有定理保证上充分光滑，满足的显式且在某区间已知baxyybayxf取定步长h，χi=χi-1+h(i=1,…)求近似值yi≈y(χi)yy=y(χ)·（χi，yi)hχ0χiχyxyxyxyxnnpnpnnnnn;,,,p,,,111121021,，多步法：单步法：数值方法分两大类：§1简单方法——解法的基础,10,100,nhfyEuleryxyyyxnnnn计算公式：法一、已知，必有切线方程。及由于斜率的切线（存在！），则出发取解曲线由，几何意义。,,f,,10000000,000yxyxyxdxdyyxfxyyyx),(f)x(xyyxxydxdyxyy,x00000000：由点斜式写出切线方程）（：，可由切线算出，则为等步长，yxyyyxxhfhh000110121011，，，）（：，点值在）（逐步计算出，nhfxyyyxyyxnnnnn注意：这是“折线法”而非“切线法”除第一个点是曲线切线外，其他点不是！yχ0χ1χ2χ3χ2，Euler法的分析与解释)(2))(,()()(2)()()()(//2//2/xyhxxxxyhxyxxxnnnnnnnnnyhfyhyhyTaylorxyA展开：的附近在,2,1,0),()()(11nhfyyhyxyyxxynnnnnnn的近似值：得且的线性部分，取xxxxdxdynnnnnnyyyx11,B用差商代替微商,2,1,0,,,,11111nhffhyyhyxyyyxyyxyxyxxnnnnnnnnnnnnnn代替，则：用C.数值积分dt)t(y,tfdtdtdy:kxxkxxkxx积分该微分方程：的途径。三种解释给出各自推广则：，取且210111,C,B,A,,,nhfh),hy,y,x,hkdt)t(y,tf)x(y)kx(yy,xyyxxyx(fyyxyxyxnnnnnnnnnnnnnnnkxx3,Euler方法的误差估计无误差！法求出的值，即由（精确等！）时其中非累积误差：在一步中产生的误差而局部截断误差。是当)()1~~1111yxyyyxTnnnnnnnEuleryyxxyhxxxxxxxnnnnnnnnnyhfyhyyTaylory1//2112))(,()()()(展开：点将在xxyhyxTxxxynnnnnnnnnyyhfy1//211112则,)(~~hhMTyMOxyxnbxa2221//22)(,)(max充分光滑，则：令2),总体方法误差(1)yxyxxxnnnnnnyL,fy,fLipschitzyxf,条件：充分光滑，或满足），（自然假定一性由微分方程解的存在唯关系。推出总体误差与步长的邻步的总体误差关系递推方法：从任意两相yyyyTyyyxyxeyxen~nn~nnn~n~nnnnnnnnyyny1111111111111n以下估计步：则对步的总体截断误差记为第2),总体方法误差(2)eyxyxxxyxyxyxxxyynnnnnnnnnnnnnnnn~n)hL(yhLLipschitz,fy,fh,hfy,hfy11y11条件由ThLThLTTeTeyxeeTeNNNNNNNnnnhLhLyhL1122110001111110Nn1则：，由，对取定成立。对一切ThLTTehTNnNNnhLO111211则，由局部截断误差1021011NkkKNNkkhLhThLOhOOhLhhLN21111无关常数与步长001100hehLlimhLlimxxLhxxhNhNN.0,0eNhhOEuler速率收敛：方法以微分方程数值解的稳定性累会淹没真值？的近似计算值，误差积是其中差：稳定性分析，对计算误1*1*111yyyynnnnn。这种方法是绝对稳定的与称对时不产生增大的误差，，，在计算若计算误差，对给定的步长是复常数。其中：定义一种数值方法求解h210,yknn/khyy则称绝对稳定区域。定的，的允许范围内是绝对稳，对hEuler法的绝对稳定区域*n*n*nnnn/yyyyyyyhhEulery计算值：算法：的11nnnh1误差方程：是绝对稳定区域当从而1111hhnn稳定的。对稳定区域称如果整个左半平面是绝越强。可选大些，方法适应性绝对稳定区域越大，AhIm(λh)-2-10Re(λh)（二）向后Euler方法hyyxxxynnn11,用向前差商：yxyxyyyhfnnnn00111,则隐式算法：hOhLOhTn1021时：当整体截断误差，，局部截断误差,,,khf,hfEulery,xyyyxyynknnknnnnn210111101法结合：与为避免解非线性方程，（二）向后Euler法的稳定性yyyynnnhEulery11/：法用向后对11nnnh误差公式：hhhnn11211111则222)(211212111hhRhhe10)(n1n则只要hRe仍受限制。要求稳定的。但收敛法是因此向后,10hhLAEuler（三）梯形公式dttytfyyxxnnnnxx11)(,由积分途径：yxyxyyxyxynnnnnnnnnnffhyy11111,2,,则得：：积分用梯形公式，且令,,,kffhhfnEuleryxyxyyyxyyknnnnnknnnnn210,2210111101,,，，，，对法结合，形成迭代算法同样与1202,,211111111111hLhLffhyyyxyxyyknknknnknnknkn收敛条件收敛性由提高精度的计算代价—函数值但每次迭代均多算一次，均高一阶，方法的局部与总体误差梯形公式比heNEuler20梯形公式的稳定性222241)(141)(12111121212hhRhhReehhhnnnnnn-A10)(1稳定的。梯形公式是，时当nnehR.Euler112ThehNNOhOO一般但多计算一次函数值。高一阶，法比，梯形法的精度§2.Runge-Kutta法A.一步法分析：Euler法与梯形法),(2合。在某些点上值的线性组，两种计算公式全是yxf同的越多，精度越高。级数，比较低阶项，相成展开及点将法：在，判断局部截断误差方)3Taylorx(yyxnnn),(合构造求解公式。在某些点上值的线性组想法：用yxf1在哪些点上的值？，计算问题：)y,x(f，线性组合系数？2用待定系数法来求这些参数，使上述3，中相同项更多KCyyiNiinn11构造B.一般公式,2,1,01,2,i111nKbyh,axhfK,yxhfKjijijnininn其中：点的使计算公式在取参数Tayloryxbacnnijii,,,,hdhddh33221n1n!31!21yyxyhxyhxyxxnnnnn!!hyhy///3//2/3121及越好。的低次幂项相同的越多时当hyyxnnKbyaxKyxKKcKcyynnnnnnhhfhfN12122122111,2,，则算法：一般公式中理论的简单例——梯形法yxTnnny~111局部截断误差：hxyhxyxxxxOhyyTaylorynnnnnn3//2/11!2展开式：点的在先求ffyfffyyxyxf////)!(~1xyxynnnnyTaylor展开，由公式点的在求xyxxKn/nnh),y(hf1KbxaxKnnyhhf12122)(,hfKbfaxyxfOhhyxnn31212)()(,hfabfhacxyccxyxxffOfhyyyxn/nnnnyx3221222211~,,取值，从而均在yxyxyxyyhTnnnnnnnnnhfhhfhfO,,21,21Euler131合的一步迭代算法：法结即梯形法和，达到1,211,21,1212212212221311baccabaccchTxOynn取则要的展开式，对比付出多少得到多少！hfffffffhffhxyxxOfffhyyyyyxyxyxxnnn4222232/12!3!2计算两次函数值还能提高误差阶数吗？多展开一项：hfKbfKbafhafbfaxyxKOhKh)(,fhyxxyyxnnn32212212122222212212!21hffabfabfhafabfhaxyxOff)(,hfyxyxyxnnn4222222122123222212222项！掉无论怎样取参数无法去2fffyyxf四阶Rung-Kutta算法KKKKyynn432112261yxKnnhf,10,1,2,n,2,22,2342312KyxKKyxKKyxKnnnnnnhhfhhfhhfhehTOONn451,则R-K方法的绝对稳定区域)()()(()()(()(4323432232121/4121)4121)212121,),(hhhyKyKhhyKyKhyKyhyKyhhhhhhhKRyyxfnnnnnnn公式：代入将