rungekutax.m

内燃机转子仿真

clear all
clc
close all
m=[2 0;0 1];
c=[0 0;0 0];
k=[6 -2;-2 4];
F=[0;10];

%n=1;
n=length(F);
%w=0;
i=1;
x0=zeros(n,1);
y0=zeros(n,1);
w=1200;
dt=60/w/360*2;
t=2;
t1=0.5;
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);
z0=inv(m)*(f0-c*y0-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);
xyz=zeros(3*n,n1);
xyz(1,i)=i;
xyz((2:n+1),i)=x0;
xyz((n+2):(2*n+1),i)=y0;
xyz((2*n+2):(3*n+1),i)=z0;
for i=2:n1
    
    f=F;
    xz=x0+y0*dt/2+z0*dt^2/8;
    yz=y0+z0*dt/2;
    zz=inv(m)*(f-c*yz-k*xz);
    for j=1:5
    xz1=x0+y0*dt/2+(2*z0+zz)*dt^2/24;
    yz1=y0+(z0+zz)*dt/4;
    zz1=inv(m)*(f-c*xz1-k*xz1);
    yz=yz1;
    zz=zz1;
    end
    x11=x0+yz1*dt;
    y11=y0+zz1*dt;
    z11=inv(m)*(f-c*x11-k*x11);
    x21=x0+dt*(y0+2*yz+2*yz1+y11)/6;
    y21=y0+dt*(z0+2*zz+2*zz1+z11)/6;
    z21=inv(m)*(f-c*y21-k*x21);
    t=(i-1)*dt;
    x1=x0+y11*dt;
    y1=y0+z11*dt;
    z1=inv(m)*(f-c*y1-k*x1);
    x0=x1;
    y0=y1;
    z0=z1;
    xyz(1,i)=i;
    xyz(2:n+1,i)=x0;
    xyz((n+2):2*n+1,i)=y0;
    xyz((2*n+2):3*n+1,i)=z0;
end
figure(1);
plot(xyz(1,:),xyz(2,:));
figure(2);
plot(xyz(1,:),xyz(3,:));