详见mathworks
龙格库塔方法
写成矩阵(状态方程)的形式更简洁一点(其实这两种方法结果是一样的,如果C是[1,0,0]的话,就很明显了)
例如:求系统在0-5s内的单位阶跃相应,已知传递函数:
状态方程:
clear
%状态方程
A=[0,1
-100,-10];
B=[0;
1];
C=[100,0];
%初值
x=[0;0];y=0;t0=0;tf=5;
h=0.05;%步长
t=t0;%从t0开始迭代,储存时间
N=round((tf-t0)/h); %迭代步数
r=1;%单位阶跃输入
for i=1:N
k1=A * x+B*r;
k2=A * (x+h*k1/2)+B*r;
k3=A * (x+h*k2/2)+B*r;
k4=A * (x+h*k3)+B*r;
x=x+h*(k1+2*k2+2*k3+k4)/6; %采用四阶龙格库塔法
y=[y,C*x]; %输出值
t=[t,t(i)+h];
end
plot(t,y);
结果:
构造传函:
M=100;%分子
N=[1,10,100];%分母
G=tf(M,N);%传函
转换:
g=ss(G);%状态方程
A=g.A;%各个系数矩阵
B=g.B;
C=g.C;
ODE45
1、 直接@
A=@sin(x);
A(pi)
2、创建匿名函数
%定义函数f(x)=x^2
f = @(x) x^2;
f(2)
- ODE45使用
[t,y] = ode45(odefun,tspan,y0)
t:自变量
y(t):因变量
odefun:函数句柄(自动y’=)
tspan:t的范围,如[0,10]
y0:初值
1、一阶方程
a=@(t,y) 2*t;
[t,y] = ode45(a, [0 5], 0);
plot(t,y)
2、二阶方程
先转化为两个一阶微分:
创建m文件:
function dydt = vdp1(t,y)
dydt = [y(2);
(1-y(1)^2)*y(2)-y(1)];%括号不可省,不是变量,y是2x1列向量
调用ode45
[t,y] = ode45(@vdp1,[0 20],[2; 0]);
plot(t,y)
输出y两列,分别为y(1)和y(2)
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)