rungeandkutta.m

内燃机转子仿真

clear all
clc
m=4;
k=2000;
x0=0;
y0=0;
w=sqrt(k/m);
t=2*pi/w;
dt=0.02;
z0=(100-k*x0)/m;
f=[100 100 100 100 80 60 40 20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
for i=1:27
    xz=x0+y0*dt/2;
    yz=y0+z0*dt/2;
    zz=(f(1,i)-k*xz)/m;
    xz1=x0+yz*dt/2;
    yz1=y0+zz*dt/2;
    zz1=(f(1,i)-k*xz1)/m;
    x11=x0+yz1*dt;
    y11=y0+zz1*dt;
    z1=(f(1,i)-k*x11)/m;
    x1=x0+dt*(y0+2*yz+2*yz1+y11)/6;
    y1=y0+dt*(z0+2*zz+2*zz1+z1)/6;
    x0=x1
    y0=y1;
    z0=(100-k*x0)/m;
end