MATLAB程序设计实践课程设计.docVIP

中南大学 MATLAB程序设计实践 班级：材料0901 姓名：杨力通 学号：0604090102 指导老师：汪冰峰 材料科学与工程学院 2012年4月10日 一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。 【实例】采用龙格-库塔法求微分方程： 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为： 其中： 2、流程图： 2.1、四阶龙格－库塔（R-K）方法流程图： 2.2、实例求解流程图： 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 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 3.2、实例求解源程序： %运行四阶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,