搜索
常微分方程组求解matlab程序代码clc,clear,closeallCT=[1;1;1];h=0.25;[k,X,Y,wucha,P]=RK4z(@dydx,0.1,60,CT,h),plot(X,Y(:,1),'g--',X,Y(:,2),'b*--',X,Y(:,3),'mp--')xlabel('轴\itx');ylabel('轴\ity')gridonlegend('方程解z1的曲线','方程解z2的曲线','方程解z3的曲线')functiondY=dydx(X,Y)dY(1)=Y(2);dY(2)=Y(3);dY(3)=Y(3)*X^(-1)-3*X^(-2)*Y(2)+2*X^(-3)*Y(1)+9*X^3*sin(X);dY=[dY(1);dY(2);dY(3)];endfunction[k,X,Y,wucha,P]=RK4z(dydx,a,b,CT,h)n=fix((b-a)/h);X=zeros(n+1,1);Y=zeros(n+1,length(CT));X=a:h:b;Y(1,:)=CT';fork=1:nk1=feval(dydx,X(k),Y(k,:))x2=X(k)+h/2;y2=Y(k,:)'+k1*h/2;k2=feval(dydx,x2,y2);k3=feval(dydx,x2,Y(k,:)'+k2*h/2);k4=feval(dydx,X(k)+h,Y(k,:)'+k3*h);Y(k+1,:)=Y(k,:)+h*(k1'+2*k2'+2*k3'+k4')/6;k=k+1;endfork=2:n+1wucha(k)=norm(Y(k)-Y(k-1));k=k+1;endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:);P=[k',X',Y,wucha'];