四阶龙格库塔RK方法求常微分方程.docxVIP

四阶龙格库塔RK方法求常微分方程 四阶龙格库塔RK方法求常微分方程 PAGE / NUMPAGES 四阶龙格库塔RK方法求常微分方程 标准文档 中南大学 MATLAB 程序设计实践 资料科学与工程学院 2013年3月 26日 一、编程实现“四阶龙格－库塔（ R-K）方法求常微分方程” ，并举一 例应用之。 【实例】采用龙格 -库塔法求微分方程： 合用文案 标准文档 y y 1 y( x0 ) 0 , x0 0 1、算法说明： 在龙格 -库塔法中，四阶龙格 -库塔法的局部截断误差约为 o(h5)，被广泛应用于解微分 方程的初值问题。其算法公式为： yn 1 yn h 2k 2 2k3 ) (k1 6 其中： k1 f ( xn , yn ) k2 f (xn 1 h, yn 1 hk1 ) 2 2 k3 f ( xn 1 h, yn 1 hk2 ) 2 2 k4 f (xn h, yn hk3 ) 2、流程图： 、四阶龙格－库塔（ R-K ）方法流程图： 输入待求微分方程、 求解的自变量范围、 初值以及求解范围内的取点数等。 确定求解范围内的步长 k = 取点数？ 结束算法 是 否 k1 f ( xn , yn ) 求解： k2 f ( xn 1 h, yn 1 hk1 ) 2 2 k3 f ( xn 1 h, yn 1 hk 2 ) 2 2 k4 f ( xn h , y n hk3 ) 求解并输出： y n 1 y n h ( k 2k 2 2k 3 ) 6 、实例求解流程图： 合用文案 标准文档 开始 输入求解的自变量范围 求出待求简单微分 方程的真值解 用 MA TLAB 自带函数 ode23 求解待求微分方程 用自编函数四阶龙格－库塔 R-K ）方法求解待求微分方程 结束 3、源程序代码 、四阶龙格－库塔（ 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 nargin4 | PointNum=0 PointNum=100; end if nargin3 y0=0; end y(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,varargin{:}); f4=f4(:); %得公式 k4 y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6; %得 y(n+1) end 、实例求解源程序： %运行四阶 R-K法 clear, clc %除去内存中的变量 x0=0; xt=2; Num=100; h=(xt-x0)/(Num-1); x=x0+[0:Num]*h; a=1; yt=1-exp(-a*x); %真值解 fun=inline( -y+1 , x , y ); %用 inline 构造函数 f(x,y) y0=0; %设定函数初值 PointNum=5; %设定取点数 [x1,y1]=ode23(fun,[0,2],0); [xr,yr]=MyRunge_Kutta(fun,x0,xt,y0,PointNum); MyRunge_Kutta_x=xr MyRunge_Kutta_y=yr plot(x,yt, k ,x1,y1, b-- ,xr,yr, r- ) legend( 真值 , ode23 , Rung-Kutta 法解 ,2) hold on plot(x1,y1, bo ,