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

第六章 常微分方程数值解法 ——RK4法、AB4法 ******（学号） *****（姓名） 上机题目要求见教材P307,23题。 一、算法原理 题目要求采用RK4法和AB4法求解最简单的常微分方程初值问题 为求解式(1)，采用离散化方法，就是寻求解在区间上的一系列点 上的近似值。记表示相邻两个节点的间距，称为步长。 求微分方程数值解的主要问题： (1) 如何将微分方程离散化,并建立求其数值解的递推公式; (2) 递推公式的局部截断误差、数值数与精确解的误差估计； (3) 递推公式的稳定性与收敛性. Runge-Kutta方法 基本思想：通过在多预报几个点求斜率，并将其加权平均作为的近似值，以此构造更高精度的计算公式。如果每步计算四次函数 的值，完全类似的，可以导出局部截断误差为的四阶Runge-Kutta公式(RK4)： Adams显式公式 Runge-Kutta方法是单步法，计算时，只用到， 而已知信息、等没有被直接利用。可以设想如果充分利用已知信息，，…来计算，那么不但有可能提高精度，而且大大减少了计算量，这就是构造所谓线性多步法的基本思想。我们知道常微分初值问题(1)与积分 是等价的。下面我们讨论最简单的多步法——Admas方法。为了导出求微分方程(1)的数值解法，我们将(3)右边积分的被积函数用插值多项式来近似替代，可以得到4阶精度的显式公式 称为AB4法。 二、流程图 四阶Runge-Kutta公式(RK4)通用程序流程图1所示。 No Yes 图1关于RK4算法流程图 AB4法通用程序流程图2所示。 No Yes 图2关于AB4算法流程图 三、计算代码 RK4完整代码 function y=Runge_Kutta4(f,y0,a,b,h) %求解常微分方程初值问题 %Input - f is the function % - a and b 定义域端点 % - y0 初始值 % - h 步长 %Output - y 输出 clc y=y0; ii=1; for x=a+h:h:b k1=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); temp=y(ii)+h*(k1+2*k2+2*k3+k4)/6; y=[y,temp]; ii=ii+1; end AB4完整代码 function y=AB4(f,y0,a,b,h) %求解常微分方程初值问题 %Input - f is the function % - a and b 定义域端点 % - y0 初始值 % - h 步长 %Output - y 输出 y=Runge_Kutta4(f,y0,a,a+3*h,h); ii=4; for x=a+3*h:h:b-h temp=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 四、计算结果及分析 对于初值问题 取步长，采用RK4法和AB4法计算，并将计算结果与精确值作比较，计算调用程序Runge_Kutta4和AB4，如下 %求解常微分方程dy/dx=-x^2*y^2的初值问题 %分别采用RK4法和AB4法求解，并与精确解作比较 clc,format long f=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法误差