常微分方程数值.

,3,2,1k第七章常微分方程数值解2第七章常微分方程数值解7.1引言(基本求解公式)§7.2Runge-Kutta法§7.3微分方程组和高阶方程解法简介§3本章要点：本章主要研究基于微积分数值解法的常微分方程数值解，主要方法有线性单步法中的Euler方法、Euler方法的改进方法、Runge-Kutta方法高阶微分方程和微分方程组的数值解法47.1引言(基本求解公式)§在工程和科学技术的实际问题中，常需要求解微分方程只有简单的和典型的微分方程可以求出解析解而在实际问题中的微分方程往往无法求出解析解在高等数学中我们见过以下常微分方程：0)(),(yaybxayxfy)(,)(),,(0ayyaybxayyxfy-----------(1)-----------(2)5nybyyaybxayyxfy)(,)(),,(0-----------(3)(1),(2)式称为初值问题,(3)式称为边值问题本课程主要研究问题(1)的数值解法我们首先介绍初值问题(1)的解存在的条件6定理1.条件，即满足如果Lipschitzyxf),(，均有使得正数],[,baxL|||),(),(|2121yyLyxfyxf的解存在且唯一则初值问题)1(对于问题(1),要求它的数值解)(],[)(节点上的一系列离散点在区间就是求未知函数baxybxxxxan210),,2,1()(nkyxykk的近似值上函数值的数值解就是问题而)1(),,2,1(nkyk70)(),(yaybxayxfy-----------(1)从(1)的表达式可以看出,求它的数值解的关键在于数值计算问题)(xy中程或者它的等价的积分方xxdttytfyxy0))(,()(0的数值计算问题积分xxdttytf0))(,(而数值微分或数值积分问题我们都已经学习过8一、基于数值微分的常微分方程数值解法0)(),(yaybxayxfy-----------(1)对于初值问题(1)bxxxxan210为了讨论方便,假设以下节点为等距节点khaxnabhk,9(一)Euler公式将在按泰勒公式展开，有1xy10,xx0201010012yxxxxxyxyxy01xxh)()()(001xyhxyxy)(202yh10--------(2)),(1jjjjyxhfyy)(2)(21jjyhhe)(22jxyh--------(3)记)(jjxyy1,,1,0nj其中11)(jjyxy)()()(1jjjxyhxyxy)(22jyh1,,1,0nj),()(jjjyxfxy(2)和(3)式称为求解初值问题(1)的(前进)Euler公式和误差项在按泰勒公式展开1,jjxx11同理，将在按泰勒公式展开，有0xy10,xx0210101102yxxxxxyxyxy01xxh)()()(101xyhxyxy)(202yh12)()()(11jjjxyhxyxy)(22jyh记),(111jjjjyxhfyy)(2)(21jjyhhe)(212jxyh)(jjxyy其中11)(jjyxy--------(4)--------(5)(4)和(5)式称为求解初值问题(1)的后退Euler公式和误差项式形公此类公式称为隐式右端含有注意)(,)4(1jy1,,1,0nj),()(111jjjyxfxy式形公此类公式称为显公式右端不含而前进)(,1jyEuler13从(2)或(4)式不难看出,jjyy时只要用到前一个值在计算1这种类型的方法称为单步格式或单步法Euler方法的几何体现:)(jjxyhy),(1jjjjyxhfyy前进Euler公式)(1jjxyhy),(111jjjjyxhfyy后退Euler公式0123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.50123456-1.5-1-0.500.511.514Euler1.m例1.公式求解初值问题用前进Euler1)0(102yxyxyy解:yxyyxf2),(显然1,1,10,000ybnax由前进Euler公式),(1jjjjyxhfyy0,1,2,,jn)2(jjjjyxyhy1.0h取1500.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8得)2(00001yxyhyy1.1)1021(1.01)2(11112yxyhyy1918.1)1.11.021.1(1.01.1依此类推,有],[yx00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前前Euler前前前前前01.00000.10001.10000.20001.19180.30001.27740.40001.35820.50001.43510.60001.50900.70001.58030.80001.64980.90001.71781.00001.784816由于后退Euler公式是隐形公式,计算例1将很麻烦事实上大多数情况下用后退Euler公式都较困难),(111jjjjyxhfyy),(1jjjjyxhfyy11jjyyEuler的预测值公式得到如果用前进11jjyEulery公式计算代入后退然后将就可得到新的Euler公式0)(yay--------(6)1,,1,0nj此方法称为预测—校正系统17用Euler公式的预测——校正系统求解例1.例2.解:由(6)式,有)2(00001yxyhyy1.1)1021(1.01)2(11101yxyhyy0918.1)1.11.021.1(1.01)2(11112yxyhyy1827.1)0918.11.020918.1(1.00918.1)2(22212yxyhyy1763.1)1827.12.021827.1(1.00918.1Euler1.m18依此类推,得01.00000.10001.09180.20001.17630.30001.25460.40001.32780.50001.39640.60001.46090.70001.52160.80001.57860.90001.63211.00001.6819],[yx比较不同的结果00.10.20.30.40.50.60.70.80.9111.11.21.31.41.51.61.71.8前前Euler前前前前Euler前前前前前校正系统预测19(二)常微分方程数值解的截断误差评价一个微分方程求解公式的标准当然是其精度:)(11的差与计算值也就是精确值jjyxy)()(111heyxyjjj而在求解公式中),(1jjjjyxhfyy误差项1,,1,0nj)(jjjxyyy都是近似值，即一般步的误差只能表示求解公式第1)()(111jyxyhejjj20定义1.部截断误差步的局的求解公式第为计算称jyyxyhejjjj)()(因为一般情况下,求解公式的每一步都存在误差,因此有定义2.且步的截断误差的求解公式第为计算设,)(jyhejjkjjkhehE1)()(步的累计截断误差为该求解公式第则称khEk)(点上的总体截断误差即该求解公式在kx定义3.)()(1pjhOhe误差为若求解公式的局部截断阶精度则称该求积公式具有p21-2-1.8-1.6-1.4-1.2-1-0.8-0.6-0.4-0.2000.511.522.53210xxxx)()(1122xyyxyy)(1he)(2he)(3he)(1hE)(2hE)(3hE22Euler公式的局部截断误差为)(2)(21jjyhhe具有1阶精度后退Euler公式的局部截断误差为)(2)(21jjyhhe也具有1阶精度显然一个求解公式的精度越高,计算解的精确性也就越好从前面的分析可知,Euler法的精度并不算高因此有必要找寻精度更高的求解公式23二、基于数值积分的常微分方程数值解法0)(),(yaybxayxfy-----------(1)对于初值问题上积分对上式在区间],[1kkxxkkxxkkdxyxfxyxy1),()()(1-----------(7)kkxxkkdxyxfxyxy1),()(,)(1就要求计算积分则计算已知假设kkkkxxxxdxyxfdxy11),(24:),(1的计算kkxxdxyxf),(),(111kkxxyxhfdxyxfkk矩形求积公式))](,(),([2),(111kkkkxxxyxfyxfhdxyxfkk已知假设11)(kkyxy梯形求积公式,误差为将以上求积公式代入(7)式,并加以处理就可得到相对应的求解公式))(,(12)(3kkkyfhTR))](,())2(,2(4),([6),(11111kkkkkkxxxyxfhxyhxfyxfhdxyxfkkSimpson求积公式,误差为))(,(180)()4(5kkkyfhSR25(一)矩形求解公式kkxxkkdxyxfxyxy1),()()(1),(),(111kkxxyxhfdxyxfkk由可得),()(111kkkkyxhfyxy令),(111kkkkyxhfyy-----------(8)(8)式称为矩形公式(矩形法)实际上就是Euler求解公式26(二)梯形求解公式kkxxkkdxyxfxyxy1),()()(1由))](,(),([2),(111kkkkxxxyxfyxfhdxyxfkk可得)])(,(),([2)(111kkkkkkxyxfyxfhyxy)],(),([2111kkkkkkyxfyxfhyy令------(9)称(9)式为梯形求解公式(梯形法)注意:(9)式是隐形公式27))(,(12)()(3kkkkyfhTRhe)()9(11kkxyy式中如果在则梯形公式第k步的截断误差为)(123kyh],[1kkkxx)(3hO显然梯形法具有二阶精度由于梯形公式为隐形公式,一般情况下不易显化kkyyEuler的预测值求出矩形法公式可以先使用)8)((28即进行校正代入梯形公式然后将,)9(ky),(111kkkkyxhfyy)],(),([2111kkkkkkyxfyxfhyy------(10)以上公式称为改进的Euler求解公式(改进Euler法),即)],(),([2111kkkkkkyxfyxfhyy))],(,(),([2111111kkkkkkkyxhfyxfyxfhy------(11)29例3.用Euler公式、梯形公式和改进Euler公式求解初值问题,并比较结果的精度1)0(]5.0,0[,yxydxdy解:1.0h取步长(1)Euler公式