《常微分方程数值解法》由会员分享，可在线阅读，更多相关《常微分方程数值解法（24页珍藏版）》请在金锄头文库上搜索。
1、-第八章 常微分方程的数值解法一容要点考虑一阶常微分方程初值问题:微分方程的数值解：设微分方程的解y(x)的存在区间是a,b，在a,b取一系列节点a= x0 x1 xn =b,其中hk=xk+1-xk；(一般采用等距节点，h=(b-a)/n称为步长)。在每个节点xk求解函数y(x)的近似值:yky(xk)，这样y0 , y1 ,.,yn称为微分方程的数值解。用数值方法，求得f(xk)的近似值yk，再用插值或拟合方法就求得y(x)的近似函数。(一)常微分方程处置问题解得存在唯一性定理对于常微分方程初值问题:如果：(1) 在的矩形是一个二元连续函数。(2) 对于y满足利普希茨条件，即那么在上方程的
2、解存在且唯一，这里C=min(A-x0),x0+B/L),L是利普希茨常数。定义：任何一个一步方法可以写为，其中称为算法的增量函数。收敛性定理：假设一步方法满足：(1)是p解的.(2) 增量函数对于y满足利普希茨条件.(3) 初始值y0是准确的。那么kh=x-x0,也就是有(一)、主要算法1局部截断误差局部截断误差：当y(xk)是准确解时，由y(xk)按照数值方法计算出来的的误差y(xk+1)-称为局部截断误差。注意：yk+1和的区别。因而局部截断误差与误差ek+1=y(xk+1) -yk+1不同。如果局部截断误差是O(hp+1),我们就说该数值方法具有p阶精度。1. 显式欧拉格式：k=1,2
3、,n-1.显式欧拉格式的特点：(1)、单步方法；(2)、显式格式；(3)、局部截断误差因而是一阶精度 。隐式欧拉格式2. 两步欧拉格式：k=1,2,n-1.两步欧拉格式的局部截断误差，因而是二阶精度3. 梯形方法:；3.改良的欧拉方法：预测值: 校正值:。或改写为:4、梯形方法与改良欧拉方法的截断误差是O(h3),具有二阶精度。5、龙格-库塔法的思想1).二阶龙格-库塔法计算公式：当:时，得一簇龙格-库塔公式，其局部截断误差均为O(h3)都是二阶精度。特别取,就是改良欧拉公式。取,得二阶龙格库塔法为：，称为二阶中点格式。2)、经典龙格库塔格式(也称为标准龙格库塔方法)：四阶龙格库塔方法的截断误
4、差为，具有四阶精度。一般一阶常微分方程初值问题用四阶龙格库塔方法计算，其精度均满足实际问题精度要求。3).变步长龙格库塔方法： 从节点xk出发，以步长h据四阶龙格库塔方法求出一个近似值,然后以步长h/2求出一个近似值，得误差事后估计式：根据来选取步长h。4). RKF格式：变步长龙格库塔方法，因频繁加倍或折半步长会浪费计算量。Felhberg改良了传统龙格库塔方法，得到RKF格式，较好解决了步长确实定，而且提高了精度与稳定性，为Matlab等许多数值计算软件采用。4/5阶RKF格式：由4阶龙格库塔方法与5阶龙格库塔方法结合而成。Felhberg得到的最正确步长hs，其中h为当前步长，.e为精度
5、要求，假设s1.5步长加倍。6.亚当斯方法(Adams)1).显式Adams方法：记：；二阶显式Adams方法：;三阶显式Adams方法：;四阶显式Adams方法：.2). 隐式Adams方法二阶隐式Adams方法：;三阶隐式Adams方法：;四阶隐式Adams方法：3.Adams预报-校正系统：先用显式格式作为预测值，再用隐式格式来校正。预测值：校正值：4).改良的Adams预报-校正系统：预测：改良：校正值：改良：7.收敛与稳定性对于固定的，如果数值解yk当同时时趋向于准确解y(xk)，那么称该数值方法方法是收敛的。如果一种数值方法，在节点值yk上大小为d的扰动，于以后各节点值ym(mk)
6、上产生的偏差均不超过d，那么称数值方法是稳定的. 一般，隐式格式较显式格式有较好的稳定性。8.刚性方程组：考虑n阶常微分方程组：，*假设矩阵A的特征值l1,l2,ln的实部Re(li)e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1;if kK,error(迭代发散);endend改良Euler格式functionx,y=naeulerg(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 K1=feval(dyfun,x(n),y(n); y(n+1)=y(n)+h*K1; K2=f
7、eval(dyfun,x(n+1),y(n+1); y(n+1)=y(n)+h*(K1+K2)/2;endx=x;y=y;4阶经典Runge-Kutta格式functionx,y=nak4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;for n=1:length(x)-1 K1=feval(dyfun,x(n),y(n); K2=feval(dyfun,x(n)+h/2,y(n)+h/2*K1); K3=feval(dyfun,x(n)+h/2,y(n)+h/2*K2); K4=feval(dyfun,x(n+1),y(n)+h*K3); y(
8、n+1)=y(n)+h*(K1+2*K2+2*K3+K4)/6;endx=x;y=y;变步长4阶经典Runge-Kutta格式functionx,y=nark4v(dyfun,xspan,y0,e,h)if nargin5,h=(xspan(2)-xspan(1)/10;endn=1;x(n)=xspan(1);y(n)=y0;y1,y2=put(dyfun,x(n),y(n),h);while x(n)e while abs(y2-y1)/10e h=h/2; y1,y2=put(dyfun,x(n),y(n),h); endelse while abs(y2-y1)/10=e h=2*h;
9、 y1,y2=put(dyfun,x(n),y(n),h); end h=h/2;h=min(h,xspan(2)-x(n);y1,y2=put(dyfun,x(n),y(n),h);end n=n+1;x(n)=x(n-1)+h;y(n)=y2; y1,y2=put(dyfun,x(n),y(n),h); endx=x;y=y;function y1,y2=put(dyfun,x,y,h); y1=rk4(dyfun,x,y,h);y21=rk4(dyfun,x,y,h/2);y2=rk4(dyfun,x+h/2,y21,h/2);function y=rk4(dyfun,x,y,h); K
10、1=feval(dyfun,x,y); K2=feval(dyfun,x+h/2,y+h/2*K1); K3=feval(dyfun,x+h/2,y+h/2*K2); K4=feval(dyfun,x+h,y+h*K3); y=y+h*(K1+2*K2+2*K3+K4)/6;五.试验例题例1 不同方法的精度比拟用多种方法解下述初值问题，并与其准确解进展比拟。解：取步长h=0.1, xk=kh (k=0,1,6)。各方法计算结果及绝对误差见表(1)、(2)、(3)。表1xn欧拉公式改良的欧拉公式四阶标准龙格库塔公式yn误差yn误差yn误差0.01.000000010100.11.0000001.00500001.004837500.21.0100001.0190251.018730900.31.0290001.0412181.040818420.41.0561001.0708021.070320290.51.0904901.1070761.106530930.61.1314411.1494041.14881193表2四阶阿当姆斯公式nxn显示公式隐式公式yn误差yn误差30.3取自准确解1.0408180140.41.070322921.0703196650.51.10
