ode45求解微分方程

用MATLAB实现的四阶龙格库塔，并用来求解微分方程组function[x,y]=Runge_kutta(ufunc,y0,h,a,b)%%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数）n=floor((b-a)/h);%求步数,迭代次数%x=zeros(n+1,1);%n+1行1列的全零矩阵x(1)=a;%时间起点y(:,1)=y0;%赋初值，可以是向量，但是要注意维数forii=1:nx(ii+1)=x(ii)+h;k1=ufunc(x(ii),y(:,ii));k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);k4=ufunc(x(ii)+h,y(:,ii)+h*k3);y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;%按照龙格库塔方法进行数值求解%y(:,ii+1)=roundn(y(:,ii+1),-8);endendclc;clearall;closeall;%-------初值和参数----------a=16.923;b=20.653;r=0.51;G=1.41;k1=1/1.355;%k1=1/2.5;k2=1.25;k3=0;x0=[0.001;0.1;0;0];%-------定义方程----------fun=@(t,y)[a*(y(2)-y(1)+G*y(1)-y(1)*(k1*((1-k2*(y(4)-k3))^(-1/2))));y(1)-y(2)+y(3);-b*y(2)-r*y(3);y(1);];%--------龙格库塔解-----------[t,x]=Runge_Kutta(fun,x0,0.01,0,200);%测试时改变test_fun的函数维数，别忘记改变初始值的维数%--------龙格库塔解-----------v=x(1,3000:end);i=((k1.*((1-k2.*(x(4,3000:end)-k3)).^(-1/2)))).*x(1,3000:end);x(1,:)=round((x(1,:)+13)*10+1);x(2,:)=round((x(2,:)+13)*10+1);x(3,:)=round((x(3,:)+13)*10+1);x(4,:)=round((x(4,:)+13)*10+1);v=fix((v+15)*10);i=fix((i+15)*10);figure(1)plot(t,x(1,:))%自编的龙格库塔函数效果xlabel('x')ylabel('y');使用方法[T,Y]=ode45(odefun,tspan,y0)odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名tspan是区间[t0tf]或者一系列散点[t0,t1,...,tf]y0是初始值向量T返回列向量的时间点Y返回对应T的求解列向量[T,Y]=ode45(odefun,tspan,y0,options)options是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等[T,Y,TE,YE,IE]=ode45(odefun,tspan,y0,options)每组(t,Y)之产生称为事件函数。每次均会检查是否函数等于零。并决定是否在零时终止运算。这可以在函数中之特性上设定。例如以events或@events产生一函数。[value,isterminal,direction]=events(t,y)其中，value(i)为函数之值，isterminal(i)=1时运算在等于零时停止，=0时继续；direction(i)=0时所有零时均需计算（默认值），+1在事件函数增加时等于零，-1在事件函数减少时等于零等状况。此外，TE,YE,IE则分别为事件发生之时间，事件发生时之答案及事件函数消失时之指针i。sol=ode45(odefun,[t0tf],y0...)sol结构体输出结果应用举例1求解一阶常微分方程程序:odefun=@(t,y)(y+3*t)/t^2;%定义函数tspan=[14];%求解区间y0=-2;%初值[t,y]=ode45(odefun,tspan,y0);plot(t,y)%作图title('t^2y''=y+3t,y(1)=-2,1t4')legend('t^2y''=y+3t')xlabel('t')ylabel('y')%精确解%dsolve('t^2*Dy=y+3*t','y(1)=-2')%ans=%(3*Ei(1)-2*exp(1))/exp(1/t)-(3*Ei(1/t))/exp(1/t)2求解高阶常微分方程关键是将高阶转为一阶,odefun的书写.F(y,y',y''...y(n-1),t)=0用变量替换,y1=y,y2=y'...注意odefun方程定义为列向量dxdy=[y(1),y(2)....]程序:functionTestode45tspan=[3.94.0];%求解区间y0=[28];%初值[t,x]=ode45(@odefun,tspan,y0);plot(t,x(:,1),'-o',t,x(:,2),'-*')legend('y1','y2')title('y''''=-t*y+e^t*y''+3sin2t')xlabel('t')ylabel('y')functiony=odefun(t,x)y=zeros(2,1);%列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);endend