四阶龙格库塔法的基本思想_四阶龙格库塔实验报告.docx

四阶龙格库塔实验报告

三、四阶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] =Runge(ydot_fun,x0,y0,h,N )x=zeros(1,N+1);y=zeros(length(y0),N+1);x(1)=x0;y(:,1)=y0;for n=1:N x(n+1)=x(n)+h; k1=h*feval(ydot_fun,x(n),y(:,n)); k2=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k1); k3=h*feval(ydot_fun,x(n)+1/2*h,y(:,n)+1/2*k2); k4=h*feval(ydot_fun,x(n)+h,y(:,n)+k3); y(:,n+1)=y(:,n)+1/6*(k1+2*k2+2*k3+k4); end3.2 程序解释及使用 该算法可以对一阶微分方程，一阶微分方程组进行有效的求解。ydot_fun为一阶微分方程的函数，x0为初始点，y0为初始向量，h为步长，N为区间的等分数，x为Xn构成的向量，y为Yn构成的矩阵。 程序调用方法：1，先编写要求解的一阶微分方程或方程组的函数文件ydot_fun.m文件，将该文件和Runge文件放到同一个目录下。2. 调用求解程序，[x,y]=Runge(@dot_fun,x0,y0,h,N)，运行后即可得出结果。或者用内部函数调用：输入：ydot_fun=(x,y)[] [x,y]= Runge(ydot_fun,x0,y0,h,N)3.3实例求解课本304页9.2题目：用标准4级4阶R-K法求解，，取步长h=0.1，计算的近似值，并与解析解作比较。解：首先将三阶方程改写成微分方程组的形式： 令得如下微分方程组在ydot_fun.m文件中编写待求解微分方程组，调用计算程序，保留5位小数得：表3-1 三阶微分方程求解结果00.10.20.30.40.50.60.70.80.91.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中第二行为原三阶微分方程对应的数值解，第二行为其一阶导数值，第三行为其二阶导数值。由结果可知，数值解y(1)=3.7183，其对应的精确解析解3相对误差为4.8947e-6，可知四阶龙格库塔法具有很高的代数精度。