4阶龙格库塔1(2)

实验中所完成的主要内容

所遇到的问题及解决方案

尚未解决的问题

一、用给定的系统模型或者自己找个系统进行仿真，要求：

1）二阶系统

2）迭代方法：定步长的4阶龙格库塔积分方法

G(s) s 1

s 3s 2 2

(0) 1, y(0) 1, u(0) 0 y

根据上次实验一求初值计算得X0=[4 1]

程序

clc; clear all; %清空%

num=[1 1]; %n维分子系数向量%

den=[1 3 2]; %n维分母系数向量%

X0=[4 1]; %状态向量初值%

V=1; %反馈系数%

n=2; %系统阶次%

T0=0; %起始时间%

Tf=50; %终止时间%

h=0.1; %计算步长%

R=1; %阶跃输入函数的幅值%

num=num/den(1);den=den/den(1);A=den(2:n+1);%首一化处理%

A=[rot90(rot90(eye(n-1:n)));-fliplr(A)]; %形成能控标准型A阵% B=[zeros(1,n-1),1]'; %形成输入阵B%

m1=length(num); %分子系数向量维数%

C=[fliplr(num),zeros(1,n-m1)]; %形成输出阵C（n维行向量）%

Ab=A-B*C*V; %形成闭环系数阵Ab%

X=X0';y=0;t=T0; %设初值，准备开始递推运算%

N=round(Tf-T0)/h; %确定输出点数%

for i=1:N %四阶龙格——库塔法%

k1=Ab*X+B*R;

k2=Ab*(X+h*k1)+B*R;

k3=Ab*(X+h*k2/2)+B*R;

k4=Ab*(X+h*k3)+B*R; %求各次斜率K%

X=X+h*(k1+2*k2+2*k3+k4)/6; %求状态%

y=[y,C*X]; %求输出并以向量形式保存%

t=[t,t(i)+h]; %输出对应时刻以向量形式保存%

end

[t',y'] %输出数据形式结果%

plot(t,y); %输出曲线形式结果%

输出结果：