[問題] 使用runge-kutta解微分方程式已回收
(3*a*Jt)*s^2+(Jt+3*a*Bm)s+Bm
有一轉移函數G(s) = ------------------------------
a^3*s^3+3*a^2*s^2+3*a*s+1
不想使用simulink以及ode45或是ode23等matlab的解題器得到系統的輸出
我的方法是將此轉移函數化成狀態方程式
然後自己寫runge-kutta，之後利用lsim指令驗證是否正確
可是兩者結果差異很大
但是我用自己寫的runge-kutta解一階的轉移函數，結果是相近的
想請問板上的各位，問題是出在哪呢???
我將轉移函數化成狀態方程式之後，寫成function
==========================函式=========================================
function xg1dot=QGpn_inv(t,ub1,xg1,AQGpn,BQGpn);
xg1dot=AQGpn*xg1+BQGpn*ub1;
==========================函式=========================================
==========================以下是主程式===============================
clc
clear
close all
%================================data======================================
Jt=0.027;%total rotating inertial---(kg*m^2)
Bm=0.016*2*pi/60;%total damping coefficient---(N*m*min/rev)
a=0.001;
AQGpn=[0 1 0;0 0 1;-1/(a^3) -3/(a^2) -3/(a)];
BQGpn=[0;0;1];
CQGpn=[Bm/(a^3) (Jt+3*a*Bm)/(a^3) 3*Jt/(a^2)];
QGpn_inv='QGpn_inv';
h=0.001;
tf=3.5-h;%模擬終止時間
t0=0;ub10=0;xg10=[0 0 0];%宣告初始值
t=t0;ub1=ub10;xg1=xg10';%先將初始值代入actuator和motor中解
i=0;%迭代因子
n=fix(tf/h)+1;%t=0算第一次(即n=1)
tout=zeros(n,1);
Tdgout=zeros(n,1);
%================================data======================================
%=================解Q/Gpn by Runge-Kutta===================================
while t<=tf+h%到t=tf停止模擬
ub1=100*sin(t);
q1=feval(QGpn_inv,t ,ub1,xg1 ,AQGpn,BQGpn);
q2=feval(QGpn_inv,t+h/2,ub1,xg1+h*q1/2,AQGpn,BQGpn);
q3=feval(QGpn_inv,t+h/2,ub1,xg1+h*q2/2,AQGpn,BQGpn);
q4=feval(QGpn_inv,t+h ,ub1,xg1+h*q3 ,AQGpn,BQGpn);
xg1=xg1+h*(q1+2*q2+2*q3+q4)/6;
Tdg=CQGpn*xg1;
i=i+1;
tout(i)=t;
Tdgout(i)=Tdg;
t=t+h;
end
%=================解Q/Gpn by Runge-Kutta===================================
%=================解Q/Gpn by lsim===========================================
t=t0:h:3.5;
ub=100*sin(t);
[yg1,xg1]=lsim(AQGpn,BQGpn,CQGpn,0,ub,t,xg10);
%=================解Q/Gpn by lsim===========================================
figure(1)
plot(t,yg1)
figure(2)
plot(tout,Tdgout)
請板上的各位不吝指教~~
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 124.8.138.176
→ fishfree:宅神原po真是帥!!! 09/16 06:23
推 lovetone:樓上真是專業!!! 09/16 12:19
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.171.123.86
※ 編輯: lovetone 來自: 118.171.123.86 (01/25 17:10)
→
01/26 11:25, , 1F
01/26 11:25, 1F