常微分方程数值解法

第九章常微分方程数值解法§1、引言00:(,)()dyfxyaxbdxyyx一阶常微分方程的初值问题'':xy-2y=4xy=f(x,y)=y(1)=-3(,)()2y4x2y4xdy2yfxy4dxxy13例方程令：且给出初值就得到一阶常微分方程的初值问题：(,)L(,)(,)fxyfxyfxyLyy只要函数适当光滑连续，且关于y满足李普希兹(Lipschitz)条件，即存在常数，使得由常微分方程理论知，初值问题的解必存在且唯一。微分方程的数值解：设方程问题的解y(x)的存在区间是[a,b]，令a=x0x1…xn=b,其中hk=xk+1-xk,如是等距节点h=(b-a)/n,h称为步长。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk上y(xk)的近似值，用yk表示，即yk≈y(xk)，这样y0,y1,...,yn称为微分方程的数值解。微分方程离散化常用方法11,A(,())nnnnnnnnyyfxyxyxdyxxxxdx用差商代替微商,2,1,0,,,,11111nhffhyyhyxyyyxyyxyxyxxnnnnnnnnnnnnnn代替，则：用11B.:(,)(0,1,)nnnnxxxxdydxfxydxndx数值积分用数值积分方法离散化1111,(),(),(,)(,)(0,1,(,))nnnnnnxnnxnnnnyyyxyxfxydxhfxyyyhfxny用代替对右端积分采用取左端点的矩形公式则有2///2//C()()()()()2()(,())()2nnnnnnnnnyxTayloryhyhyhfyxhyyxxxxhyxxxx在附近的展开：00(,)()dyfxyaxbdxyyx11()()(,)0,1,2,nnnnnnnhyyhfnTayloryxxyyyx取的线性部分，且得的近似值：展开法不仅可得到求数值解的公式，且容易估计截断误差。§2尤拉（Eular)方法000000,00001,,f,,yyxfyxyxdyyxdxyyxx，几何意义。由出发取曲线的切线（存在！），则斜率由于及已知，必有切线方程。10000,01,(,)()()nnnnEulerdyhfn,fxyaxbdxyyyyyxyxyx一、公式计算公式：),(f)x(xyyxxydxdyxyy,x00000000：由点斜式写出切线方程）（：，可由切线算出，则为等步长，yxyyyxxhfhh000110121011，，，）（：，点值在）（逐步计算出，nhfxyyyxyyxnnnnn用分段的折线逼近函数,此为“折线法”而非“切线法”，除第一个点是曲线上的切线，其它都不是。2、Euler方法的误差估计11111)~~()nnnnnnnyyEulerxTyyxyy局部截断误差。在一步中产生的误差而非累积误差：其中（精确解！）时由法求出的值，即无误差！是当xxyhxxxxxxxnnnnnnnnnyhfyhyyTaylory1//2112))(,()()()(展开：点将在112//1111~,~(),~2nnnnnnnnnnnnnyhfyyhfyyxyxxxyhyxxxTy则hhMTyMOxyxnbxa2221//22)(,)(max充分光滑，则：令3、总体方法误差(1),,,nnnnnnfxyLipschitzfyfLyyyxxxx递推方法：从任意两相邻步的总体误差关系推出总体误差与步长的关系。由微分方程解的存在唯一性自然假定（，）充分光滑，或满足条件：1111111111111n1~~~~~,nnnnnnnnnnnnnnnnnnnynyyyhfyyexyyyexxTyyyyxxxyy第步的总体截断误差记为则对步：以下估计其中总体方法误差(2)11~,,y,,1(1)nnnnnnnnnnnnnnnnnyhfyhfhfyfLipschitzhLyhLyyyxxxxyyyxxxxyxe由条件110001211211nN01111nnnNNNNNNNhLyhLhLeeTyexeeThLhLTTTT对一切成立。对取定，由，则：211111nNnNONhLhThLeTTT由局部截断误差，则1120011NNNKkkkkOhLhLhT12111NOOhhLhLh00001limlim1NhhLxxNNhxxhhLhLe与步长无关常数.0,0eNhhOEuler速率收敛：方法以4、微分方程数值解的稳定性累会淹没真值？的近似计算值，误差积是其中差：稳定性分析，对计算误1*1*111yyyynnnnn'nnkn+kn=,012hyyyhk定义：一种数值方法求解其中是复常数。对给定的步长，若计算误差在计算，，时不产生增大的误差，即，称对与这种方法是绝对稳定的。则称绝对稳定区域。定的，的允许范围内是绝对稳，对hEuler法的绝对稳定区域*n*n*nnnn/yyyyyyyhhEulery计算值：算法：的11nnnh1误差方程：是绝对稳定区域当从而1111hhnn稳定的。对稳定区域称如果整个左半平面是绝越强。可选大些，方法适应性绝对稳定区域越大，AhIm(λh)-2-10Re(λh)二、向后（后退的）Euler方法,11nnnyyhxxyx用向后差商：11100,nnnnhfyyyyxyx为隐式算法：01111E,0,1,2,1,nnnnkknnnulerhfhfknyyyxyyyx为避免解非线性方程，与法结合：11111111111,,01kkkknnnnnnkknnhffhLhLyyyyxxyy收敛条件210hL1nOOhhT局部截断误差，整体截断误差，当时：01111,0,1,2,1,nnnnkknnnhfhfknyyyxyyyx向后Euler方法收敛条件与截断误差向后Euler法的稳定性yyyynnnhEulery11/：法用向后对11nnnh误差公式：hhhnn11211111则222)(211212111hhRhhen1n()01heR只要则仍受限制。要求稳定的。但收敛法是因此向后,10hhLAEuler三、梯形公式11,nnnnxyyfxydtxxx由积分途径：yxyxyyxyxynnnnnnnnnnffhyy11111,2,,则得：：积分用梯形公式，且令,,,kffhhfnEuleryxyxyyyxyyknnnnnknnnnn210,2210111101,,，，，，对法结合，形成迭代算法同样与(,)dyfxydx（）20NEulereh梯形公式比法的局部与总体误差均高一阶，，但每次迭代均多算一次函数值—提高精度的计算代价。011111012,0122,,nnnnkknnnnnnnhfhffk,,,yyyxyyyyxx由迭代算法，对，，，11111111111,,20122kkkknnnnnnkknnhffhLhLyyyyxxyy收敛条件梯形公式的收敛性梯形公式的稳定性222241)(141)(12111121212hhRhhReehhhnnnnnn-A10)(1稳定的。梯形公式是，时当nnehR111,2,nnnnnnhffyyyyxx/112nnnhynyyyyy对用梯形公式：）（四、改进的尤拉公式梯形公式虽然提高了精度，但使算法复杂。而在实际计算中只迭代一次，这样建立的预测—校正系统称作改进的尤拉公式。1111111(,);[(,)(,)],2[(,)(,(,))],2(,);(,);()/2.nnnnnnnnnnnnnnnnnnpnnncnnpnpcyyhfxyhyyfxyfxyhyyfxyfxhyhfxyyyhfxyyyhfxyyyy预测校正这一计算公式也可表为或表为下列平均化形式'1112(01);(0)1.0.12().2();2();1().2nnnnnnpnnnncnppnpcxyyxyyhxyyhyyxyyhyyxyyhyyyyy例：用尤拉公式和改进的尤拉公式解初值问题解：取步长，尤拉公式为：改进的尤拉公式为：计算结果略。11(,);(,);()/2.pnnncnnpnpcyyhfxyyyhfxyyyy尤拉两步公式'111111()()(),22(,).nnnnnnnnnyxyxyxhyyhfxyyyn-1n上面介绍的改进尤拉公式，其预测公式（尤拉公式）的精度差，与校正公式（梯形公式）不匹配。用中心差商替代即称其为尤拉两步公式（首先要确定y,y两个值)。用尤拉两步公式与梯形公式相匹配，得到下列预测—校正系统：预测1112(,);[(,)(,)].2nnnnnnnnhfxyhyyfxyfxy校正113'11(),(),()().3nnnnnnnyyxyyxhyxyyx考虑上面给定预测－校正系统的局部截断误差：对预测公式：假定其局部截断误差：2//134///(4)1()()(,())23!2nnnnnnnnnyyhyhfyhyxxxxxxhhyyxxx''111234''''''(4)'234''''''(4)2(,)()2()()2()()()()()()2()23!4!()()()()()23!4!nnnnnnnnnnnnnnnnnnnyyhfxyyxhyxyxhhyxhhhyxhyxyxyxyxhyxhhhyxhyxyxyxyx111111111111()1.4():4()();51()().5nnnnnnnnnnnnyxyyxyyxyyyyxyyy由上两式有由此导出下列事后估计式1113'11(),()().12nnnnnnyyyxhyxyyx对校正公式：假定是准确的，即其局部截断误差为：'1112''''''''23''''''[(,)(,)]()[()2