东南大学_数值分析_第六章_常微分方程数值解法

东南大学《数值分析》上机练习——算法与程序设计实验报告第六章常微分方程数值解法——RK4法、AB4法******（学号）*****（姓名）上机题目要求见教材P307,23题。一、算法原理题目要求采用RK4法和AB4法求解最简单的常微分方程初值问题(,),()yfxyaxbya(1)为求解式(1)，采用离散化方法，就是寻求解)(xy在区间],[ba上的一系列点nxxxx321上的近似值,,,,21nyyy。记1(1,2,)iiihxxi表示相邻两个节点的间距，称为步长。求微分方程数值解的主要问题：(1)如何将微分方程(,)yfxy离散化,并建立求其数值解的递推公式;(2)递推公式的局部截断误差、数值数ny与精确解)(nxy的误差估计；(3)递推公式的稳定性与收敛性.a)Runge-Kutta方法基本思想：通过在1[,]iixx多预报几个点求斜率，并将其加权平均作为k的近似值，以此构造更高精度的计算公式。如果每步计算四次函数的值，完全类似的，可以导出局部截断误差为)(5hO的四阶Runge-Kutta公式(RK4)：1123412132431(22),6(,),(,),221(,),22(,).nnnnnnnnnnyykkkkkfxyhhkfxykhkfxhykkfxhyhk(2)b)Adams显式公式第六章常微分方程数值解法Runge-Kutta方法是单步法，计算1ny时，只用到ny，而已知信息1ny、2ny等没有被直接利用。可以设想如果充分利用已知信息1ny，2ny，…来计算1ny，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。我们知道常微分初值问题(1)与积分11()()(,())dnnxnnxyxyxftytt(3)是等价的。下面我们讨论最简单的多步法——Admas方法。为了导出求微分方程(1)的数值解法，我们将(3)右边积分的被积函数用插值多项式来近似替代，可以得到4阶精度的显式公式11122331[55(,)59(,)37(,)9(,)]24nnnnnnnnnnyyhfxyfxyfxyfxy(4)称为AB4法。二、流程图1.四阶Runge-Kutta公式(RK4)通用程序流程图1所示。NoYes图1关于RK4算法流程图开始确定输入参数(初值问题)计算k1,k2,k3,k4采用公式(2),计算yi+1计算完毕结束东南大学《数值分析》上机练习——算法与程序设计实验报告2.AB4法通用程序流程图2所示。NoYes图2关于AB4算法流程图三、计算代码1.RK4完整代码functiony=Runge_Kutta4(f,y0,a,b,h)%求解常微分方程初值问题%Input-fisthefunction%-aandb定义域端点%-y0初始值%-h步长%Output-y输出clcy=y0;ii=1;forx=a+h:h:bk1=feval(f,x,y(ii));k2=feval(f,x+0.5*h,y(ii)+0.5*h*k1);k3=feval(f,x+0.5*h,y(ii)+0.5*h*k2);k4=feval(f,x+h,y(ii)+h*k3);开始确定输入参数(初值问题)调用程序RK4，计算初值采用公式(4)，计算yi+1计算完毕结束第六章常微分方程数值解法temp=y(ii)+h*(k1+2*k2+2*k3+k4)/6;y=[y,temp];ii=ii+1;end2.AB4完整代码functiony=AB4(f,y0,a,b,h)%求解常微分方程初值问题%Input-fisthefunction%-aandb定义域端点%-y0初始值%-h步长%Output-y输出y=Runge_Kutta4(f,y0,a,a+3*h,h);ii=4;forx=a+3*h:h:b-htemp=y(ii)+h*(55*feval(f,x,y(ii))-59*feval(f,x-h,y(ii-1))+...37*feval(f,x-2*h,y(ii-2))-9*feval(f,x-3*h,y(ii-3)))/24;y=[y,temp];ii=ii+1;end四、计算结果及分析对于初值问题22(01.5)(0)3yxyxy(5)取步长0.1h，采用RK4法和AB4法计算，并将计算结果与精确值3()3/(1)yxx作比较，计算调用程序Runge_Kutta4和AB4，如下%求解常微分方程dy/dx=-x^2*y^2的初值问题%分别采用RK4法和AB4法求解，并与精确解作比较clc,formatlongf=inline('-x*x*y*y');%函数f(x,y)y_rk4=Runge_Kutta4(f,3,0,1.5,0.1);%RK4法disp('RK4法计算结果为：')fprintf('%8.5f',y_rk4)fprintf('\n')y_ab4=AB4(f,3,0,1.5,0.1);%AB4法disp('AB4法计算结果为：')fprintf('%8.5f',y_ab4)fprintf('\n')x=0:0.1:1.5;东南大学《数值分析》上机练习——算法与程序设计实验报告y=3./(1+x.^3);%精确解disp('精确结果为：')fprintf('%8.5f',y)fprintf('\n')error_rk4=abs(y-y_rk4);%RK4法误差error_ab4=abs(y-y_ab4);%AB4法误差disp('RK4法的误差为：')fprintf('%8.5f',error_rk4)fprintf('\n')disp('AB4法的误差为：')fprintf('%8.5f',error_ab4)fprintf('\n')在MATLAB中运行以上程序，计算结果如表1所示表1RK4法与AB4法计算结果及其误差ix()iyxRK4法AB4法iy()iiyxyiy()iiyxy03.000003.000000.000003.000000.000000.12.997002.997000.000002.997000.000000.22.976192.976190.000002.976190.000000.32.921132.921130.000002.921130.000000.42.819552.819550.000002.818390.001160.52.666672.666660.000002.664670.001990.62.467112.467100.000012.465200.001900.72.233802.233800.000012.233080.000730.81.984131.984120.000001.984950.000820.91.735111.735110.000001.737040.001941.01.500001.500010.000011.502190.002191.11.287001.287010.000011.288760.001761.21.099711.099720.000021.100720.001021.30.938380.938400.000020.938710.000331.40.801280.801300.000020.801130.000151.50.685710.685730.000020.685330.00038从表1可以看出，Rung-Kula方法和AB4法均能很好的实现常微分方程初值问题的解法，且具有较高的精度和较低的运算复杂度。同时，根据计算的误差可见，RK4法较之于AB4法具有更高的计算精度，但计算量比AB4法稍微多一些。五、结论本文介绍了常微分方程初值问题的两种解法：单步法（Rung-Kula方法）和多步法（Adams方法）。在实际问题中，如果要求解的精度高，常用的是四阶第六章常微分方程数值解法Rung-Kula方法，这是因为四阶Rung-Kula方法精度高，编程简单，易于调节步长，计算过程稳定；其缺点是计算量较大，计算函数(,)fxy的次数比同阶的线性多步法（四阶Adams方法）多几倍；另外，如果(,)fxy的光滑性较差，Rung-Kula方法的精度还不如改进的Euler方法高。一般地，若(,)fxy较复杂，用Adams方法比较合适。不论采用哪种方法，选取的步长应使落在绝对稳定区域内。一般在保证精度的条件下，步长尽可能大些，这样可以节省计算量。