richardson.m

内燃机转子仿真

%clear all;
%clc;
%close all;
function Tmax=richardson(n,t0,N,step,T,T1,f0,K,dK,zhuansu,M,C)
T11=zeros(2*n,1);
T21=zeros(2*n,1);
T22=zeros(2*n,1);
Tmax=zeros(2*n,1);
T11=step*(T1+T)*f0/2;
for k=1:N;
    dH=zeros(2*n,1);
    for j=1:2^k-1
        K=K+dK*sin(pi*zhuansu/30*(t0+j*step/2^k));
        H=[-inv(M)*C/2 inv(M);0.25*C*inv(M)*C-K -0.5*C*inv(M)];
        I=eye(size(H));
        Tk=H*(t0-j*step/2^k)+(H*(t0-j*step/2^k))^2*(I+(H*(t0-j*step/2^k))/3+(H*(t0-j*step/2^k))^2/12+(H*(t0-j*step/2^k))^3/60)/2; 
        dH=dH+0.5*step/2^k*Tk*f0;
    end
T21=0.5*T11+dH;
for i=2:k+1
T22=4^i*T21/(4^i-1)-T11/(4^i-1);
end 
T11=T21;
cp=T22-T11;
if cp<=1e-2
    Tmax=T22;
    break;
end
end