四阶龙格库塔实验报告

三、四阶Runge-Kutta法求解常微分方程一、龙格库塔法的思想根据第九章的知识可知道，Euler方法的局部截断误差是2()Oh，而当用Euler方法估计出1,()(1)nnnnyyhfxy再用梯形公式111[(,)(,)](2)2nnnnnnhyyfxyfxy进行校正，即采用改进Euler方法得出数值解的截断误差为3()Oh。由Lagrange微分中值定理'11()()()()()(,())(3)nnnnnyxyxyxxyxhfy记*(,())khfy，得到*1()()(4)nnyxyxk这样只要给出一种计算*k的算法，就能得到相应的计算公式。用这种观点的来分析Euler方法和改进Euler方法，Euler方法的迭代公式可改写为111(,)nnnnyykkhfxy改进Euler方法的预报-校正公式可改写为1121211()2(,),(,)nnnnnnyykkkhfxykhfxhykEuler方法实际上是用一个点处的值1k近似*k，而改进Euler方法是用两个点处的值1k，和2k，做算术平均值近似*k自然改进Euler方法要优于Euler方法。因此，可以想到假如在1[,]nnxx内多预报几个点值ik，并用他们的加权平均值作为*k的近似值，则有可能构造出具有更高精度的计算公式，这就是Runge-Kutta法的基本思想。二、四阶龙格库塔法由Runge-Kutta的基本思想，构造四阶Runge-Kutta法是利用1234,,kkkk和的加权平均值来近似*k，因此令1112233441211132211243312213(,)(,)(5)(,)(,)nnnnnnnnnnyywKwKwKwKKhfxyKhfxhyKKhfxhyKKKhfxhyKKK使得511()()nnyxyOh即其总体截断误差为4()Oh。采用泰勒公式展开，经过复杂的推导，得到一个具有13个参数，11个方程的线性方程组。由于方程的个数少于未知量的个数，因此方程有无穷多个解。可以根据情况得到几种常用的解，即得到相应的四阶公式。最常见的四阶公式如式（6）：，112341213243(22)6(,)1(,)(6)221(,)22(,)nnnnnnnnnnhyyKKKKKhfxyhKhfxyKhKhfxyKKhfxhyK也称为标准四阶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*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法求解，'''''''''23(0)1,(0)3,(0)2yyyyxyyy，取步长h=0.1，计算(1)y的近似值，并与解析解()21xyxxex作比较。解：首先将三阶方程改写成微分方程组的形式：令'''''12132,,yyyyyyyy得如下微分方程组'12'23'''3211'''12323(0)1,(0)3,(0)2yyyyyyyyxyyy在ydot_fun.m文件中编写待求解微分方程组，调用计算程序，保留5位小数得：表3-1三阶微分方程求解结果x00.10.20.30.40.50.60.70.80.91.01y-1.00000-0.68948-0.355720.004960.396730.824361.293271.809622.380423.013633.718272y3.000003.215693.465683.754814.088554.473084.915385.423376.005966.673237.436553y2.000002.320862.687083.104673.580384.121804.737505.437126.231507.132848.15483表3-1中第二行1y为原三阶微分方程对应的数值解，第二行2y为其一阶导数值，第三行3y为其二阶导数值。由结果可知，数值解y(1)=3.7183，其对应的精确解析解(1)1ye3.71828188的相对误差为4.8947e-6，可知四阶龙格库塔法具有很高的代数精度。