rungeandkutta1.m

内燃机转子仿真

clear all
clc
close all
m=[2 0;0 1];
c=[0 0;0 0];
k=[6 -2;-2 4];
o=zeros(2);
F=[0;10;o(:,1)];
M=[c m;m o];
K=[k o;o -m];

%n=1;
n=length(F);
%w=0;
i=1;
x0=zeros(n,1);
y0=zeros(n,1);
w=1200;
%dt=60/w/360*2;
dt=0.0281;
t=2.81;
t1=1;
%dn=rem(t,dt);
%dn1=t-dn;
%n1=dn1/dt;
%k1=0.5;
%dn1=rem(t1,dt);
%dn2=t1-dn;
%n2=dn2/dt;
f0=zeros(n,1);
f=zeros(n,1);
f1=inv(M)*(f0-K*x0);
xz=zeros(n,1);
yz=zeros(n,1);
zz=zeros(n,1);
xz1=zeros(n,1);
yz1=zeros(n,1);
zz1=zeros(n,1);
x11=zeros(n,1);
y11=zeros(n,1);
z11=zeros(n,1);
x1=zeros(n,1);
y1=zeros(n,1);
n1=1000;
%p=[0.0002;0.0002];
xyz=zeros(2*n+1,n1);
xyz(1,i)=i;
xyz((2:n+1),i)=x0;
xyz((n+2):(2*n+1),i)=y0;
t=0;
for i=2:n1
    t=t+dt;
    %if i<=n2
       % f=F*(i-1)*dt;
    %else
        %f=F;
    %end
    %dk=k1*sin(w*(i-1)*dt);
    %k(2,2)=k(2,2)+dk;
    f=F;
    xz=x0+y0*dt/2;%预测中值
    yz=y0+f1*dt/2;
    f2=inv(M)*(f-K*xz);
    for j=1:5                                     %校正中值
    xz1=x0+yz*dt/2;
    yz1=y0+f1*dt/2;
    f3=inv(M)*(f-K*xz);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*dt;
    y11=y0+f3*dt;
    f4=inv(M)*(f-K*y11);
    x21=x0+dt*(y0+2*yz+2*yz1+y11)/6;        %二阶龙格库塔
    y21=y0+dt*(f1+2*f2+2*f3+f4)/6;
    x1=x0+x21*dt
    y1=y0+y21*dt;
    x0=x1;
    y0=y1;
    inv(M)*(f0-K*x0);
    xyz(1,i)=i;
    xyz(2:5,i)=x0;    
end
figure(1);
plot(xyz(1,:),xyz(2,:));
figure(2);
plot(xyz(1,:),xyz(3,:));