matlab课程设计及plc课程设计.docVIP

中南大学 MATLAB程序设计实践 班级： 姓名： 学号： 指导老师： 材料科学与工程学院 一、编程实现“四阶龙格－库塔（R-K）方法求常微分方程”，并举一例应用之。 【实例】采用龙格-库塔法求微分方程： 1、算法说明： 在龙格-库塔法中，四阶龙格-库塔法的局部截断误差约为o(h5)，被广泛应用于解微分方程的初值问题。其算法公式为： 其中： 2、流程图： 2.1、四阶龙格－库塔（R-K）方法流程图： 输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。 输入待求微分方程、求解的自变量范围、初值以及求解范围内的取点数等。 确定求解范围内的步长 k = 取点数？ 否 求解： 求解并输出： 是 结束算法 2.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 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;