SOS! PTVP trong matlab. Bất kì ai trả lời cũng được vote

SOS! PTVP trong matlab. Bất kì ai trả lời cũng được vote

E đang có chút thắc mắc về giải hệ ptvp phi tuyến trong matlab.
Các bác giúp e nhé.
Bác nào giúp e, e sẽ gửi mail mfile đến.
Xin cảm ơn và hậu tạ vote5*

Bạn có thắc mắc gì thì cứ post lên đây mới có nhiều chuyên gia để ý chứ!
WJT.

E gửi các bác 2 file mfile. Các bác xem giúp nhé.
Khi e chạy file pl1 thi nó báo lỗi như sau:
??? Input argument ''t'' is undefined.
Error in ==> D:MATLAB6p5workQTQD.m
On line 7 ==> uA=sqrt(2)*Um*sin(2*pi*fv*t);
Error in ==> D:MATLAB6p5workPL1.M
On line 37 ==> [t,si] = ode45(QTQD,[0 0.7],[0;0;0;0;0]);

Đây là file chính
%clc
%clear
global r1 F1 x10 x20 Xo ktg Dx2s Omeo r20 r21 is p G ...
a1 a2 b1 b2 s X1 X2 r2 x2s Mc Um fv;
a1 = 0; a2 = 0; b1 = 0; b2 = 0;
is = 0; f = 50;
x10 = 0.98; %so lieu dong co
x20 = 0.92; %so lieu dong co
x2s = 0.81; %so lieu dong co
x2bh = 0.755; %so lieu dong co
r20 = 0.235; %so lieu dong co
r21 = 0.376; %so lieu dong co
r1 = 0.4; %so lieu dong co
Xo = 22; %so lieu dong co
p = 2; %so lieu dong co,so doi cuc
J = 0.125; %so lieu dong co, momen quan tinh
ktg = 1/115;
Dx2s = x20 - x2s;
Omeo = 314; %Toc do goc dong bo omega_0
F1 = p/J;
fv = 50; %Tan so can thay doi
%----------------------------------------------------------------------------------
%U/f khong doi, Mc = const
ku = 220/f;
Um = ku*fv;
%----------------------------------------------------------------------------------
%U/f^2 khong doi, Mc ~ Omega^2
%ku = 220/f^2;
%Um = ku*fv^2;
%----------------------------------------------------------------------------------
%U/sqrt(f^3) khong doi, Mc ~ Omega
%ku = 220/sqrt(f^3);
%Um = ku*sqrt(fv^3);
%----------------------------------------------------------------------------------
[t,si] = ode45(QTQD,[0 0.7],[0;0;0;0;0]);
uA=sqrt(2)*Um*sin(2*pi*fv*t'');
uB=sqrt(2)*Um*sin(2*pi*fv*t''+2*pi/3);
uC=sqrt(2)*Um*sin(2*pi*fv*t''-2*pi/3);
u1_anpha = uA;
u1_beta = (sqrt(3)/3)*(uC-uB);
i1_anpha = si,1)*a1 - si,3)*a2;
i2_anpha = si,3)*b1 - si,1)*b2;
i1_beta = si,2)*a1 - si,4)*a2;
i2_beta = si,4)*b1 - si,2)*b2;
M = G*(si,3).*i1_beta - si,4).*i1_anpha);
Di1_anpha = a1*u1_anpha - a1*r1*i1_anpha + a2*r2*i2_anpha + a2*si,5).*si,4);
Di1_beta = a1*u1_beta - a1*r1*i1_beta + a2*r2*i2_beta - a2*si,5).*si,3);
Di2_anpha = -b2*(u1_anpha - r1*i1_anpha) - b1*(r2*i2_anpha + si,5).*si,4));
Di2_beta = -b1*(u1_beta-r1*i1_beta) + b1*(-r2*i2_beta + si,5).*si,3));
pd1 = r1*(i1_anpha.^2 + i1_beta.^2);
pd2 = r2*(i2_anpha.^2 + i2_beta.^2);
pd = pd1 + pd2;
im_anpha = i1_anpha + i2_anpha;
im_beta = i1_beta + i2_beta;
Dim_anpha = Di1_anpha + Di2_anpha;
Dim_beta = Di1_beta + Di2_beta;
pdg1 = (X1/Omeo)*(i1_anpha.*Di1_anpha + i1_beta.*Di1_beta);
pdg2 = (X2/Omeo)*(i2_anpha.*Di2_anpha + i2_beta.*Di2_beta);
pdg = pdg1 + pdg2;
pm = (Xo/Omeo)*(im_anpha.*Dim_anpha + im_beta.*Dim_beta);
P1 = u1_anpha.*i1_anpha+u1_beta.*i1_beta;
Tp = pd + pdg + pm;
temp4 = size(t);
kt = temp4(1,1);
%Tinh thoi gian tich phan
Count = 0;
for i=1:kt
if si(i,5) < 0.95*si(kt,5);
Count = Count + 1;
else
break
end
end

%Tinh dien nang tieu thu
for i=1:kt
if t(i) < t(count)
temp(i) = t(i);
end
end
temp2 = size(temp);
Dtp = temp2(1,2);
A1 = trapz(t(1tp),P1(1tp));
Adg = trapz(t(1tp),pdg(1tp));
ATp = trapz(t(1tp),Tp(1tp));
Ad = trapz(t(1tp),pd(1tp));
Phan_tram_Adg_A1 = Adg/A1;
Phan_tram_Ad_A1 = Ad/A1;
Phan_tram_Adg_ATp = Adg/ATp;
%----------------------------------------------------------------------------------
figure(1);
plot(t,si,5),t,M);
hold on;
title([''Thoi diem dong co o che do xac lap: '' num2str(t(Dtp))]);
xlabel(''thoi gian'');
ylabel(''MOMEN va TOC DO'');
%legend(''Tocdo'',''Momen'');
grid on;
figure(2);
subplot(2,1,1);
plot(t,u1_anpha);
hold on;
title(''Dien ap pha thuc'');
xlabel(''Thoi gian (giay)'');
ylabel(''Dien ap (V)'');
grid on;
subplot(2,1,2);
plot(t,i1_anpha);
hold on;
title(''Dong dien pha A stato'');
xlabel(''Thoi gian (giay)'');
ylabel(''Dong dien (A)'');
grid on;
figure(3);
subplot(3,1,1);
plot(t, pdg);
hold on;
title(''Ton hao qua do'');
xlabel(''Thoi gian (giay)'');
ylabel(''Cong suat (W)'');
grid on;
subplot(3,1,2);
plot(t, pd);
hold on;
title(''Ton hao dong trong qua trinh qua do'');
xlabel(''Thoi gian (giay)'');
ylabel(''Cong suat (W)'');
grid on;
subplot(3,1,3);
plot(t, Tp);
hold on;
title(''Tong ton hao trong qua trinh qua do'');
xlabel(''Thoi gian (giay)'');
ylabel(''Cong suat (W)'');
grid on;
Đây là chương trình con
unction dsi = QTQD(t,si);
global r1 F1 i1_anpha i1_beta x10 Xo ktg x20 Dx2s Omeo r20 r21 is p G ...
a1 a2 b1 b2 X1 X2 r2 s x2s X1 X2 Mc Um fv;

% Tinh cac tham so
uA=sqrt(2)*Um*sin(2*pi*fv*t);
uB=sqrt(2)*Um*sin(2*pi*fv*t+2*pi/3);
uC=sqrt(2)*Um*sin(2*pi*fv*t-2*pi/3);
u1_anpha = uA;
u1_beta = (sqrt(3)/3)*(uC-uB);
% Tinh toan phan tham so dong
s = (1/Omeo)*(Omeo - si(5));
is = sqrt((si(1)*a1 - si(3)*a2)^2 + (si(2)*a1 - si(4)*a2)^2);
X1 = x10 - 1.15*(1e-2)*ktg*is - 0.25;
X2 = x20 - 0.11*(1e-2)*ktg*is - s*Dx2s;
r2 = r20 + s*(r21 - r20);
%--------------------------------------------------------------------------
Xs = Xo + X1;
Xr = Xo + X2;
a1 = (Omeo*Xr) / (Xs*Xr - Xo^2);
a2 = a1*(Xo/Xr);
b1 = (Omeo*Xs) / (Xs*Xr - Xo^2);
b2 = b1*(Xo/Xs);
kr = Xo/Xr;
G = (3*p*kr)/2;
% Xac dinh phu tai can tinh
Mc = 0;
%Mc = 92.2;
%Mc = 0.00093*si(5)^2;
%Mc = 0.2936*si(5);
%--------------------------------------------------------------------------
i1_anpha = si(1)*a1 - si(3)*a2;
i2_anpha = si(3)*b1 - si(1)*b2;
i1_beta = si(2)*a1 - si(4)*a2;
i2_beta = si(4)*b1 - si(2)*b2;
M = G*(si(3).*(si(2)*a1 - si(4)*a2) - si(4).*(si(1)*a1 - si(3)*a2));
dsi = zeros(5,1);
% Xac dinh he ham vi phan
dsi(1) = u1_anpha - r1*i1_anpha;
dsi(2) = u1_beta - r1*i1_beta;
dsi(3) = -r2*i2_anpha - si(5)*si(4);
dsi(4) = -r2*i2_beta - si(5)*si(3);
dsi(5) = F1*(M-Mc);
Được chutmayman sửa chữa / chuyển vào 21:23 ngày 26/10/2005

Các chuyên gia giúp em với.

CT của bạn bị lặp va số loi dùng biến thế nên không chạy được là phải.
VD:
- Lời gọi của hàm ode45 tham chiếu tới hàm QTQD, trong hàm này lại dùng biến u1_alpha, biến được tính phía sau lời gọi hàm ode45.
- Trong lời gọi của QTQD cần có tham số cho nó (cái này do bạn đã định nghĩa cho QTQD).
Tớ copy lại và trình bầy qua CT của bạn. Tuy nhiên nó vẫn chưa chạy được đâu, bạn xem lại nhé.
function[] = un()
global r1 F1 x10 x20 Xo ktg Dx2s Omeo r20 r21 is p G ...
a1 a2 b1 b2 s X1 X2 r2 x2s Mc Um fv;
a1 = 0; a2 = 0; b1 = 0; b2 = 0;
is = 0; f = 50;
x10 = 0.98; %so lieu dong co
x20 = 0.92; %so lieu dong co
x2s = 0.81; %so lieu dong co
x2bh = 0.755; %so lieu dong co
r20 = 0.235; %so lieu dong co
r21 = 0.376; %so lieu dong co
r1 = 0.4; %so lieu dong co
Xo = 22; %so lieu dong co
p = 2; %so lieu dong co,so doi cuc
J = 0.125; %so lieu dong co, momen quan tinh
ktg = 1/115;
Dx2s = x20 - x2s;
Omeo = 314; %Toc do goc dong bo omega_0
F1 = p/J;
fv = 50; %Tan so can thay doi
%----------------------------------------------------------------------------------
%U/f khong doi, Mc = const
ku = 220/f;
Um = ku*fv;
%----------------------------------------------------------------------------------
%U/f^2 khong doi, Mc ~ Omega^2
%ku = 220/f^2;
%Um = ku*fv^2;
%----------------------------------------------------------------------------------
%U/sqrt(f^3) khong doi, Mc ~ Omega
%ku = 220/sqrt(f^3);
%Um = ku*sqrt(fv^3);
%----------------------------------------------------------------------------------
[t,si] = ode45(QTQD,[0, 0.7],[0,0,0,0,0]); %%%%% THAM SO CUA QTQD
uA=sqrt(2)*Um*sin(2*pi*fv*t'''');
uB=sqrt(2)*Um*sin(2*pi*fv*t''''+2*pi/3);
uC=sqrt(2)*Um*sin(2*pi*fv*t''''-2*pi/3);
u1_anpha = uA;
u1_beta = (sqrt(3)/3)*(uC-uB);
i1_anpha = si,1)*a1 - si,3)*a2;
i2_anpha = si,3)*b1 - si,1)*b2;
i1_beta = si,2)*a1 - si,4)*a2;
i2_beta = si,4)*b1 - si,2)*b2;
M = G*(si,3).*i1_beta - si,4).*i1_anpha);
Di1_anpha = a1*u1_anpha - a1*r1*i1_anpha + a2*r2*i2_anpha + a2*si,5).*si,4);
Di1_beta = a1*u1_beta - a1*r1*i1_beta + a2*r2*i2_beta - a2*si,5).*si,3);
Di2_anpha = -b2*(u1_anpha - r1*i1_anpha) - b1*(r2*i2_anpha + si,5).*si,4));
Di2_beta = -b1*(u1_beta-r1*i1_beta) + b1*(-r2*i2_beta + si,5).*si,3));
pd1 = r1*(i1_anpha.^2 + i1_beta.^2);
pd2 = r2*(i2_anpha.^2 + i2_beta.^2);
pd = pd1 + pd2;
im_anpha = i1_anpha + i2_anpha;
im_beta = i1_beta + i2_beta;
Dim_anpha = Di1_anpha + Di2_anpha;
Dim_beta = Di1_beta + Di2_beta;
pdg1 = (X1/Omeo)*(i1_anpha.*Di1_anpha + i1_beta.*Di1_beta);
pdg2 = (X2/Omeo)*(i2_anpha.*Di2_anpha + i2_beta.*Di2_beta);
pdg = pdg1 + pdg2;
pm = (Xo/Omeo)*(im_anpha.*Dim_anpha + im_beta.*Dim_beta);
P1 = u1_anpha.*i1_anpha+u1_beta.*i1_beta;
Tp = pd + pdg + pm;
temp4 = size(t);
kt = temp4(1,1);
%Tinh thoi gian tich phan
Count = 0;
for i=1:kt
if si(i,5) < 0.95*si(kt,5);
Count = Count + 1;
else
break
end
end
%Tinh dien nang tieu thu
for i=1:kt
if t(i) < t(count)
temp(i) = t(i);
end
end
temp2 = size(temp);
Dtp = temp2(1,2);
A1 = trapz(t(1tp),P1(1tp));
Adg = trapz(t(1tp),pdg(1tp));
ATp = trapz(t(1tp),Tp(1tp));
Ad = trapz(t(1tp),pd(1tp));
Phan_tram_Adg_A1 = Adg/A1;
Phan_tram_Ad_A1 = Ad/A1;
Phan_tram_Adg_ATp = Adg/ATp;
%----------------------------------------------------------------------------------
figure(1);
plot(t,si,5),t,M);
hold on;
title([''''Thoi diem dong co o che do xac lap: '''' num2str(t(Dtp))]);
xlabel(''''thoi gian'''');
ylabel(''''MOMEN va TOC DO'''');
%legend(''''Tocdo'''',''''Momen'''');
grid on;
figure(2);
subplot(2,1,1);
plot(t,u1_anpha);
hold on;
title(''''Dien ap pha thuc'''');
xlabel(''''Thoi gian (giay)'''');
ylabel(''''Dien ap (V)'''');
grid on;
subplot(2,1,2);
plot(t,i1_anpha);
hold on;
title(''''Dong dien pha A stato'''');
xlabel(''''Thoi gian (giay)'''');
ylabel(''''Dong dien (A)'''');
grid on;
figure(3);
subplot(3,1,1);
plot(t, pdg);
hold on;
title(''''Ton hao qua do'''');
xlabel(''''Thoi gian (giay)'''');
ylabel(''''Cong suat (W)'''');
grid on;
subplot(3,1,2);
plot(t, pd);
hold on;
title(''''Ton hao dong trong qua trinh qua do'''');
xlabel(''''Thoi gian (giay)'''');
ylabel(''''Cong suat (W)'''');
grid on;
subplot(3,1,3);
plot(t, Tp);
hold on;
title(''''Tong ton hao trong qua trinh qua do'''');
xlabel(''''Thoi gian (giay)'''');
ylabel(''''Cong suat (W)'''');
grid on;
%??y l? chuong tr?nh con
function[dsi] = QTQD(t,si);
global r1 F1 i1_anpha i1_beta x10 Xo ktg x20 Dx2s Omeo r20 r21 is p G a1 a2 b1 b2 X1 X2 r2 s x2s X1 X2 Mc Um fv;
% Tinh cac tham so
uA=sqrt(2)*Um*sin(2*pi*fv*t);
uB=sqrt(2)*Um*sin(2*pi*fv*t+2*pi/3);
uC=sqrt(2)*Um*sin(2*pi*fv*t-2*pi/3);
u1_anpha = uA;
u1_beta = (sqrt(3)/3)*(uC-uB);
% Tinh toan phan tham so dong
s = (1/Omeo)*(Omeo - si(5));
is = sqrt((si(1)*a1 - si(3)*a2)^2 + (si(2)*a1 - si(4)*a2)^2);
X1 = x10 - 1.15*(1e-2)*ktg*is - 0.25;
X2 = x20 - 0.11*(1e-2)*ktg*is - s*Dx2s;
r2 = r20 + s*(r21 - r20);
%--------------------------------------------------------------------------
Xs = Xo + X1;
Xr = Xo + X2;
a1 = (Omeo*Xr) / (Xs*Xr - Xo^2);
a2 = a1*(Xo/Xr);
b1 = (Omeo*Xs) / (Xs*Xr - Xo^2);
b2 = b1*(Xo/Xs);
kr = Xo/Xr;
G = (3*p*kr)/2;
% Xac dinh phu tai can tinh
Mc = 0;
%Mc = 92.2;
%Mc = 0.00093*si(5)^2;
%Mc = 0.2936*si(5);
%--------------------------------------------------------------------------
i1_anpha = si(1)*a1 - si(3)*a2;
i2_anpha = si(3)*b1 - si(1)*b2;
i1_beta = si(2)*a1 - si(4)*a2;
i2_beta = si(4)*b1 - si(2)*b2;
M = G*(si(3).*(si(2)*a1 - si(4)*a2) - si(4).*(si(1)*a1 - si(3)*a2));
dsi = zeros(5,1);
% Xac dinh he ham vi phan
dsi(1) = u1_anpha - r1*i1_anpha;
dsi(2) = u1_beta - r1*i1_beta;
dsi(3) = -r2*i2_anpha - si(5)*si(4);
dsi(4) = -r2*i2_beta - si(5)*si(3);
dsi(5) = F1*(M-Mc);
Được le_rob sửa chữa / chuyển vào 01:41 ngày 02/11/2005

Cảm ơn nhiều nhé.
Bạn có thể cho địa chỉ liên hệ để mình trao đổi trực tiếp k?
Mà bạn ở đâu nhỉ?

Hi, tớ không ở VN có gì bạn cứ PS trong TTVNOL cho tớ, thỉnh thoảng tớ vào đây chơi mà.

[/quote]
Hi, tớ không ở VN có gì bạn cứ PS trong TTVNOL cho tớ, thỉnh thoảng tớ vào đây chơi mà.
[/quote]
Hoặc bạn gửi cho tớ tới hộp thư: let_le_rob@yahoo.com