内燃机转子仿真

				clear all
				clc
				fni=input('the name of the film:','s');
				fid=fopen(fni,'r');
				n=fscanf(fid,'%d',1);
				F0=fscanf(fid,'%f',[n,1]);
				x10=fscanf(fid,'%f',[n,1]);
				x20=fscanf(fid,'%f',[n,1]);
				w=fscanf(fid,'%f',1);
				b=fscanf(fid,'%f',1);
				K=fscanf(fid,'%f',[n,n]);
				M=fscanf(fid,'%f',[n,n]);
				C=fscanf(fid,'%f',[n,n]);
				pn=fscanf(fid,'%d',1);
				panduan=fscanf(fid,'%d',1);
				fclose(fid);
				K1=zeros(n);
				dt=0.0005*pi/40;
				endt=pi/40;%
				K1=K+(6/b/dt^2)*M+(3/b/dt)*C;
				%对K1作ldlt分解
				K2=K1;
				k1=zeros(n-1,1);
				k=zeros(n,1);
				for j=1:n;
				    for i=1:j-1 
				    k1(i)=K1(j,i)*K1(i,i);
				    end
				    k(j)=K1(j,j)-K1(j,1:j-1)*k1(1:j-1);
				    K1(j,j)=k(j);
				    K1(j+1:n,j)=(K1(j+1:n,j)-K1(j+1:n,1:j-1)*k(1:j-1))/k(j);
				end
				l=tril(K1,0)-diag(diag(K1))*diag(ones(n,1))+diag(ones(n,1));
				D=diag(diag(K1))*diag(ones(n,1));
				fni=input('the name of the rusl:','s');
				fid=fopen(fni,'wt');
				%开始循环      
				for i=0:dt:endt    
				 f0=F0*sin(w*i);
				x30=inv(M)*(f0-C*x20-K*x10);
				 f=F0*sin(w*i)+b*F0*[sin(w*(i+dt))-sin(w*i)];   
				f1=f+M*((6/b/dt^2)*x10+(6/b/dt)*x20+2*x30)+C*(3/b/dt)*x10+2*x20+(b*dt/2)*x30;
				x11=inv(l')*inv(D)*inv(l)*f1;
				x3=(6/b^3/dt^2)*(x11-x10-x20*b*dt)-(6/b^2/dt)*x20+(1-3/b)*x30;
				x2=x20+(dt/2)*(x3+x30);
				x1=x10+dt*x20+(dt^2/6)*(x3+2*x30);
				x10=x1;
				x20=x2;
				x30=x3;
				for j=1:3  
				    if j				fprintf(fid,'%5.4f      ',i,x1(j,1));  
				    else
				   fprintf(fid,'%5.4f\n',x1(j,1));   
				    end
				end
				end
				fclose(fid);