rungeandkutta2.m

内燃机转子仿真

clear all
clc
close all
mm=[5.98000 1.08000 1.04000 2.91300 2.91300 2.91300 2.91300 2.91300 51.46300 0.60000 1.11500 3.94400];
cc=[0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000];
kk=[0.00000 8.20000 392.20000 150.00000 112.78000 112.78000 112.78000 112.78000 112.78000 169.66000 0.50000 0.50000 50.29000 0.00000];
FF=[0.00000 20.00000 0.00000 40.00000 0.00000 30.00000 0.00000 50.00000 0.00000 0.00000 0.00000 0.00000];
%n=1;
n=length(FF);
m=zeros(n);
c=zeros(n);
k=zeros(n);
for i=1:n
    m(i,i)=mm(i);
end
kk=kk*1.0e5
for i=1:n
    k(i,i)=kk(i)+kk(i+1);
    c(i,i)=cc(i)+cc(i+1);
end
for i=1:n-1
            k(i,i+1)=-kk(i+1);
            c(i,i+1)=-cc(i+1);
            k(i+1,i)=-kk(i+1);
            c(i+1,i)=-cc(i+1);
end
Is=eye(n);
o=zeros(n,1);
os=zeros(n);
ms=inv(m);
A=[os Is;-ms*k -ms*c];
B=[os;ms];
h=0.001;
n1=5000;
w=30;
ft=zeros(n,1);
t=0;
i=1;
x=zeros(2*n,1);
xy(1,i)=i;
xy((2:2*n+1),i)=x;
ft=FF';
for i=2:n1
    
    
    f=A*x+B*ft*sin(w*t);
    f1=h*f;
    f2=h*(A*(x+0.5*f)+B*ft*sin(w*(t+0.5*h)));

    f3=h*(A*(x+0.5*f2)+B*ft*sin(w*(t+0.5+h)));

     f4=h*(A*(x+f3)+B*ft*sin(w*t));
    
    y=x+(f1+2*f2+2*f3+f4)/6;
    x=y;
    t=t+h;
    xy(1,i)=i;
    xy((2:2*n+1),i)=y;
end
figure(1);
plot(xy(1,:),xy(2,:));
figure(2);
plot(xy(1,:),xy(3,:));