欧拉法与matlab,龙格-库塔法与欧拉法的比较-Matlab仿真

%数值积分方法求解连续系统输出(龙格-库塔法与欧拉法的比较)

clear all; close all;

h=0.3; L=15/h; %h计算步长、L为仿真步数

z=[-1 -2]; p=[-4 -0.5+j -0.5-j]; k=2.5; %对象的零极点型

[A,B,C,D]=zp2ss(z,p,k); %转化为状态空间型

u=1*ones(L,1); u0=0; %输入及初值

n=length(p); %对象阶次

x0=zeros(n,1); xe0=zeros(n,1); %x0、xe0分别为龙格-库塔法、欧拉法的状态初值

for i=1:L

time(i)=i*h;

%欧拉法

xe=xe0+h*(A*xe0+B*u0);

ye(i)=C*xe;

%龙格-库塔法

k1=A*x0+B*u0;

k2=A*(x0+h*k1/2)+B*u0;

k3=A*(x0+h*k2/2)+B*u0;

k4=A*(x0+h*k3)+B*u(i);

x=x0+h*(k1+2*k2+2*k3+k4)/6;

y(i)=C*x;

%更新数据

u0=u(i);x0=x;xe0=xe;

end

plot(time,u,'k-.',time,ye,':',time,y,'r');

xlabel('t'); ylabel('y_r(t)、y(t)');

legend('y_r(t)','Euler:y(t)','Runge-Kutta:y(t)');