M_1.zip_定步长积分

a=input('n+1维分母系数向量：'); b=input('n+1维分子系数向量：'); X0=input('初值：'); V=input('反馈系数：'); n=input('系统阶数：'); T0=input('起始时间：'); Tf=input('终止时间；'); h=input('计算步长；'); R=input('阶跃输入函数的幅值:'); b=b/a(1); a=a/a(1); A=a(2: n+1); A=[rot90(rot90(eye(n-1,n)));-fliplr(A)]; B=[zeros(1,n-1),1]'; m1=length(b); C=[fliplr(b),zeros(1,n-m1)]; Ab=A-B*C*V; X=X0'; y=0; t=T0; N=round(Tf-T0)/h for i=1:N Kg=Ab*X+B*R; Kh=Ab*(X+h*Kg/2)+B*R; Ki=Ab*(X+h*Kh/2)+B*R; Kj=Ab*(X+h*Ki)+B*R; X=X+h*(Kg+2*Kh+2*Ki+Kj)/6; y=[y,C*X]; t=[t,t(i)+h]; end [t',y'] plot(t,y)