Scilab圖顯示不正確的ode值

Question

我試圖在scilab中創建一個異步電動機模型，並顯示轉速，電流和扭矩如何隨時間變化的圖形。它看起來很長，但你不需要全部閱讀。

fHz = 50; 
Um = 230; 
p = 3; 
we = 2*%pi*fHz/p; 
wb = 2*%pi*50; 
Rs = 0.435; 
Rr = 0.64; 
Ls = 0.0477; 
Xls = wb*Ls; // [Ohm] 
Lr = 0.0577; 
Xlr = wb*Lr; // [Ohm] 
Lm = 0.012; 
Xm = wb*Lm; // [Ohm] 
Xml = 1/(1/Xls + 1/Xm + 1/Xlr) // [Ohm]; 
D = 0.0002; 
J = 0.28; 
Mt = 0.0; 

function [xdot]=AszinkronGep(t, x, Um, fHz) 

xdot = zeros(12, 1); 

Fsq = x(1); 
Fsd = x(2); 
Frq = x(3); 
Frd = x(4); 
wr = x(5); 
isabc(1) = x(6); 
isabc(2) = x(7); 
isabc(3) = x(8); 
irabc(1) = x(9); 
irabc(2) = x(10); 
irabc(3) = x(11); 

Ua = Um*sin(2*%pi*fHz*t); 
Ub = Um*sin(2*%pi*fHz*t - 2*%pi/3); 
Uc = Um*sin(2*%pi*fHz*t + 2*%pi/3); 

Uab = 2/3*[1, -0.5, -0.5; 0, sqrt(3)/2, -sqrt(3)/2]*[Ua;Ub;Uc]; 

phi = 2*%pi*fHz*t; 
Udq = [cos(phi), sin(phi); -sin(phi), cos(phi)]*Uab; 
Usd = Udq(1); 
Usq = Udq(2); 
Urd = 0; 
Urq = 0; 

isd = (Fsd-Xml*(Fsd/Xls + Frd/Xlr))/Xls; 
isq = (Fsq-Xml*(Fsq/Xls + Frq/Xlr))/Xls; 
ird = (Frd-Xml*(Fsd/Xls + Frd/Xlr))/Xlr; 
irq = (Frq-Xml*(Fsq/Xls + Frq/Xlr))/Xlr; 

isdq = [isd; isq]; 
isalphabeta = [cos(phi), -sin(phi); sin(phi), cos(phi)]*isdq; 
isabc = [1, 0; -0.5, sqrt(3)/2; -0.5, -sqrt(3)/2]*isalphabeta; 

irdq = [ird; irq]; 
iralphabeta = [cos(phi), -sin(phi); sin(phi), cos(phi)]*irdq; 
irabc = [1, 0; -0.5, sqrt(3)/2; -0.5, -sqrt(3)/2]*iralphabeta; 

//TORQUE 
Me = (3/2)*p*(Fsd*isq - Fsq*isd)/wb 

Fmq = Xml*(Fsq/Xls + Frq /Xlr); 
Fmd = Xml*(Fsd/Xls + Frd /Xlr); 

//Differential equations 
xdot(1) = wb*(Usq - we/wb*Fsd + Rs/Xls*(Fmq - Fsq)); 
xdot(2) = wb*(Usd + we/wb*Fsq + Rs/Xls*(Fmd - Fsd)); 
xdot(3) = wb*(Urq - (we - wr)/wb*Frd + Rr/Xlr *(Fmq - Frq)); 
xdot(4) = wb*(Urd + (we - wr)/wb*Frq + Rr/Xlr *(Fmd - Frd)); 
xdot(5) = p*(Me - D*wr - Mt)/J; 
xdot(6) = isabc(1); 
xdot(7) = isabc(2); 
xdot(8) = isabc(3); 
xdot(9) = irabc(1); 
xdot(10) = irabc(2); 
xdot(11) = irabc(3); 
xdot(12) = Me; 

if t <= 5 then 
disp(Me); 
end 

endfunction 


//Simulation parameter 
t = 0:0.001:5; 
t0 = 0; 

//Starting parameters 
y0 = [0;0;0;0;0;0;0;0;0;0;0;0] 
y = ode(y0,t0,t,list(AszinkronGep,Um,fHz)); 

//Graphs 
figure(1) 
plot(t,y(5,:), "linewidth", 3); 
xlabel("time [s]", "fontsize", 3, "color", "blue"); 
ylabel("rpm [rpm]", "fontsize", 3, "color", "blue"); 

figure(4) 
plot(t,y(12,:), "linewidth", 3); 
xlabel("time [s]", "fontsize", 3, "color", "blue"); 
ylabel("torque [Nm]", "fontsize", 3, "color", "blue"); 

我想要一張顯示'我'作爲時間函數的圖。所以我寫：xdot（12）= Me，然後繪製它，但它看起來不像它應該如何。爲了檢查，我在函數的末尾添加了「disp（Me）」，以查看計算是否正確。是的，這些是正確的價值。爲什麼我繪製它時會給我不同的價值？

Answer 1

正如在評論中指出的，y(12)是Me的積分超過t。 如果你想Me，你只需要區分開來：

//after you run ode() and have y values 
h = t(2) - t(1); 
Me = diff(y(12,:)) ./ h; 

//plotting 
scf(); clf(); 

subplot(2,1,1); 
xtitle("y(12,:)"); 
plot2d(t,y(12,:)); 

subplot(2,1,2); 
xtitle("Me"); 
plot2d(t(1:$-1),Me); 

這裏是輸出：