四阶龙格库塔实验报告.docx

《四阶龙格库塔实验报告.docx》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、三、四阶Runge-Kutta法求解常微分方程一、龙格库塔法的思想根据第九章的知识可知道，Euler方法的局部截断误差是，而当用Euler方法估计出再用梯形公式进行校正，即采用改进Euler方法得出数值解的截断误差为。由Lagrange微分中值定理记，得到这样只要给出一种计算的算法，就能得到相应的计算公式。用这种观点的来分析Euler方法和改进Euler方法，Euler方法的迭代公式可改写为改进Euler方法的预报-校正公式可改写为Euler方法实际上是用一个点处的值近似，而改进Euler方法是用

2、两个点处的值，和，做算术平均值近似自然改进Euler方法要优于Euler方法。因此，可以想到假如在内多预报几个点值，并用他们的加权平均值作为的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta法的基本思想。二、四阶龙格库塔法由Runge-Kutta的基本思想，构造四阶Runge-Kutta法是利用的加权平均值来近似，因此令使得即其总体截断误差为。采用泰勒公式展开，经过复杂的推导，得到一个具有13个参数，11个方程的线性方程组。由于方程的个数少于未知量的个数，因此方程有无穷多

3、个解。可以根据情况得到几种常用的解，即得到相应的四阶公式。最常见的四阶公式如式（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;forn=1:Nx(n+1)=x(n)+h;k1=h*feval(ydot_fun,x(n),y(:,n));k2=h*fev

4、al(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为

5、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，计算的近似值，并与解析解作比

6、较。解：首先将三阶方程改写成微分方程组的形式：令得如下微分方程组在ydot_fun.m文件中编写待求解微分方程组，调用计算程序，保留5位小数得：表3-1三阶微分方程求解结果00.10.20.30.40.50.60.70.80.91.0-1.00000-0.68948-0.355720.004960.396730.824361.293271.809622.380423.013633.718273.000003.215693.465683.754814.088554.473084.915385.423

7、376.005966.673237.436552.000002.320862.687083.104673.580384.121804.737505.437126.231507.132848.15483表3-1中第二行为原三阶微分方程对应的数值解，第二行为其一阶导数值，第三行为其二阶导数值。由结果可知，数值解y(1)=3.7183，其对应的精确解析解3.71828188的相对误差为4.8947e-6，可知四阶龙格库塔法具有很高的代数精度。