四阶龙格库塔实验报告

三、四阶Runge-Kutta法求解常微分方程 一、龙格库塔法的思想 依据第九章的学问可知道，Euler方法的局部截断误差是，而当用Euler方法估计出 再用梯形公式 进行校正，即采纳改进Euler方法得出数值解的截断误差为。 由Lagrange微分中值定理 记，得到 这样只要给出一种计算的算法，就能得到相应的计算公式。 用这种观点的来分析Euler方法和改进Euler方法，Euler方法的迭代公式可改写为 改进Euler方法的预报-校正公式可改写为 Euler方法事实上是用一个点处的值近似，而改进Euler方法是用两个点处的值，和，做算术平均值近似自然改进Euler方法要优于Euler方法。 因此，可以想到假如在内多预报几个点值，并用他们的加权平均值作为的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta法的基本思想。 二、四阶龙格库塔法 由Runge-Kutta的基本思想，构造四阶Runge-Kutta法是利用的加权平均值来近似，因此令 使得 即其总体截断误差为。 采纳泰勒公式绽开，经过困难的推导，得到一个具有13个参数，11个方程的线性方程组。由于方程的个数少于未知量的个数，因此方程有无穷多个解。可以依据状况得到几种常用的解，即得到相应的四阶公式。最常见的四阶公式如式（6） ， 也称为标准四阶Runge-Kutta法。 三、四阶龙格库塔法程序说明及应用 3.1龙格库塔的计算程序 function [x,y] Rungeydot_fun,x0,y0,h,N xzeros1,N1;yzeroslengthy0,N1; x1x0;y,1y0; for n1N xn1xnh; k1h*fydot_fun,xn,y,n; k2h*fydot_fun,xn1/2*h,y,n1/2*k1; k3h*fydot_fun,xn1/2*h,y,n1/2*k2; k4h*fydot_fun,xnh,y,nk3; y,n1y,n1/6*k12*k22*k3k4; end 3.2 程序说明及运用 该算法可以对一阶微分方程，一阶微分方程组进行有效的求解。ydot_fun为一阶微分方程的函数，x0为初始点，y0为初始向量，h为步长，N为区间的等分数，x为Xn构成的向量，y为Yn构成的矩阵。 程序调用方法 1，先编写要求解的一阶微分方程或方程组的函数文件ydot_fun.m文件，将该文件和Runge文件放到同一个书目下。 2. 调用求解程序，[x,y]Rungedot_fun,x0,y0,h,N，运行后即可得出结果。 或者用内部函数调用 输入ydot_funx,y[] [x,y] Rungeydot_fun,x0,y0,h,N 3.3实例求解 课本304页9.2题目 用标准4级4阶R-K法求解 ，， 取步长h0.1，计算的近似值，并与解析解作比较。 解首先将三阶方程改写成微分方程组的形式 令得如下微分方程组 在ydot_fun.m文件中编写待求解微分方程组，调用计算程序，保留5位小数得 表3-1 三阶微分方程求解结果 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 -1.00000 -0.68948 -0.35572 0.00496 0.39673 0.82436 1.29327 1.80962 2.38042 3.01363 3.71827 3.00000 3.21569 3.46568 3.75481 4.08855 4.47308 4.91538 5.42337 6.00596 6.67323 7.43655 2.00000 2.32086 2.68708 3.10467 3.58038 4.12180 4.73750 5.43712 6.23150 7.13284 8.15483 表3-1中其次行为原三阶微分方程对应的数值解，其次行为其一阶导数值，第三行为其二阶导数值。 由结果可知，数值解y13.7183，其对应的精确解析解3.71828188的相对误差为4.8947e-6，可知四阶龙格库塔法具有很高的代数精度。