四阶龙格-库塔R-K方法求常微分方程

实用标准文案中南大学MATLAB 程序设计实践材料科学与工程学院2013 年 3 月 26 日一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。【实例】采用龙格-库塔法求微分方程：文档实用标准文案y' y 1y(x0)  0 , x0  01、算法说明：在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为 o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为：hyn1  yn (k1  2k2  2k3)6其中：k1 k2 k 3k 4f (xn, yn )11h, yn hk1)2211f (xn h, yn hk2)22f (xn  h, yn  hk3)f (xn 2、流程图：2.1、四阶龙格－库塔（R-K）方法流程图：输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。确定求解范围内的步长k = 取点数？是否k1 求解：k 2k 3k 4f (xn, yn )11h, yn hk1)2211f (xn h, yn hk2)22f (xn  h, yn  hk3)f (xn 结束算法求解并输出：hyn1  yn (k1  2k2  2k3)62.2、实例求解流程图：文档实用标准文案开始输入求解的自变量范围求出待求简单微分方程的真值解用 MATLAB 自带函数 ode23求解待求微分方程用自编函数四阶龙格－库塔（R-K）方法求解待求微分方程结束3、源程序代码3.1、四阶龙格－库塔（R-K）方法源程序：function [x,y] = MyRunge_Kutta(fun,x0,xt,y0,PointNum,varargin)%Runge-Kutta 方法解微分方程形为 y'(t)=f(x,y(x))%此程序可解高阶的微分方程。只要将其形式写为上述微分方程的向量形式%函数 f(x,y): fun%自变量的初值和终值：x0, xt%y0表示函数在x0处的值，输入初值为列向量形式%自变量在[x0,xt]上取的点数：PointNum%varargin为可输入项，可传适当参数给函数f(x,y)%x：所取的点的x值%y：对应点上的函数值if nargin<4 | PointNum<=0 PointNum=100;endif nargin<3 y0=0;endy(1,:)=y0(:)';%初值存为行向量形式h=(xt-x0)/(PointNum-1);%计算步长x=x0+[0:(PointNum-1)]'*h;%得x向量值for k=1:(PointNum)%迭代计算文档实用标准文案 f1=h*feval(fun,x(k),y(k,:),varargin{:}); f1=f1(:)';%得公式k1 f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin{:}); f2=f2(:)';%得公式k2 f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin{:}); f3=f3(:)';%得公式k3 f4=h*feval(fun,x(k)+h,y(k,:)+f3,var...