四阶Runge-Kutta法求解常微分方程电气与自动化工程学院1011203007秦.doc

四阶Runge-Kutta法求解常微分方程 电气与自动化工程学院 1011203007 秦超 算法原理 龙格-库塔(Runge-Kutta)法是一类用于求解常微分方程重要的隐式或显式迭代法，由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。Runge-Kutta法具有精度高，收敛稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，在常微分方程的求解中得到了广泛的应用。 对于常微分方程，n阶Runge-Kutta法的一般形式为： ， 其中，，，， 龙格库塔法家族中的 其中，，，， 下一个值由现在的值加上时间间隔和估算的斜率的乘积决定。 是时间段开始时的斜率； 是时间段中点的斜率，通过欧拉法采用斜率来决定在点的值； 也是中点的斜率，用斜率决定值；是时间段终点的斜率，其值用决定。 当四个斜率取平均时，中点的斜率有更大的权值： 编程实现 四阶Runge-kutta法的matlab源程序如下所示： % % ts：待求时间区间 % u0：初值 % h：步长 % 输出参数说明： % t：离散化时间序列 % u：微分方程数值解 function [t, u] = RK4(f, ts, h, u0) t = ts(1):h:ts(2); u(1) = u0; for k = 1:length(t)-1 k1 = feval(f, t(k), u(k)); k2 = feval(f, t(k) + h/2, u(k) + h*k1/2); k3 = feval(f, t(k) + h/2, u(k) + h*k2/2); k4 = feval(f, t(k) + h, u(k) + h*k3); u(k+1) = u(k) + h*(k1 + 2*k2 + 2*k3 + k4)/6; end 算例验证 通过以下常微分方程来验证所编算法程序的正确性： 该方程的精确解为，在采用四阶Runge-Kutta法进行求解时，步长h取为0.05，程序调用代码如下： ts = [0 1]; h = 0.05; f = @(t,u)(u – 2*t/u); u0 = 1; [t, u] = RK4(f, ts, h, u0); 解析解的计算： us = sqrt(1+2*t); 所得计算结果如下表所示： t u/数值解 us/解析解 t u/数值解 us/解析解 0.05 1.0488 1.0488 0.55 1.4491 1.4491 0.1 1.0954 1.0954 0.6 1.4832 1.4832 0.15 1.1402 1.1402 0.65 1.5166 1.5166 0.2 1.1832 1.1832 0.7 1.5492 1.5492 0.25 1.2247 1.2247 0.75 1.5811 1.5811 0.3 1.2649 1.2649 0.8 1.6125 1.6125 0.35 1.3038 1.3038 0.85 1.6432 1.6432 0.4 1.3416 1.3416 0.9 1.6733 1.6733 0.45 1.3784 1.3784 0.95 1.7029 1.7029 0.5 1.4142 1.4142 1.0 1.7321 1.7321 由上表结果可知，通过四阶Runge-Kutta法所求得的结果与解析解保持很好的一致性，因此验证了该方法在求解常微分方程中的有效性和所编写程序的正确性。