Mạch phi tuyến động: Các phương pháp số tối ưu

(cid:37)(cid:172)(cid:44)(cid:3)(cid:42)(cid:44)(cid:1190)(cid:49)(cid:42)(cid:3)(cid:48)(cid:208)(cid:49)(cid:3)(cid:48)(cid:1188)(cid:38)(cid:43)(cid:3)(cid:264)(cid:44)(cid:1226)(cid:49)(cid:3)(cid:21)(cid:37)

9.5 Maïch phi tuyeán ñoäng vaø caùc PP soá

,...,

x t , ) n x t , ) n

f x x ( 1 1 2 f x x ( 2 1 2

(cid:132) Caùc baøi toaùn coù 2 phaàn : 1. Vieát phöông trình traïng thaùi. 2. Söû duïng thuaät toaùn laëp. (cid:132) Caùc thuaät toaùn laëp :

Coù hai nhoùm :

,...,

, ) x t n

( f x x n 1 2

' x (cid:173) (cid:32) 1 (cid:176) ' x (cid:32)(cid:176) 2 (cid:174) ... (cid:176) (cid:176) (cid:32)(cid:175) ' x n

1. Phöông phaùp khai trieån Taylor : tìm

chuoãi gaàn ñuùng haøm taïi x(tn).

(cid:132) Caùc bieán traïng thaùi :

2. Phöông phaùp ña thöùc hoùa : tìm ña thöùc gaàn ñuùng giaù trò x(tn) theo p giaù trò tröôùc ñoù .

9.5.1 Phöông trình traïng thaùi : (cid:132) Phöông trình traïng thaùi coù daïng : x’ = f(x,t) , x0 = x(t0) , hay neáu coù n bieán :

Tuï ñieän: Choïn uC hay qC. Cuoän daây: Choïn iL hay (cid:92)L. (cid:132) Phöông trình traïng thaùi thöôøng ñöôïc thieát laäp theo caùc phöông phaùp : doøng nhaùnh , theá nuùt hay maéc löôùi.

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

1 9.5.2 Phöông phaùp khai trieån Taylor

1. Giaûi thuaät Euler thuaän : x , n

x n

(cid:14)

(cid:14) (cid:32) 1

h f x nh . (cid:11) n

(cid:12)

2. Giaûi thuaät Euler ngöôïc : n ,(

(cid:132) Khai trieån giaù trò taïi böôùc tính thöù (n+1) : x(tn+1) theo giaù trò taïi böôùc tính thöù (n) : x(tn) baèng chuoãi Taylor vaø choïn ñeán baäc p thích hôïp (sai soá h(p+1)) , vôùi h laø böôùc tính :

x (cid:32) (cid:14) n

(cid:14)

x 1 n (cid:14)

h 1) (cid:12)

)

( x t

h f x . (cid:11) 1 n (cid:14) Giaûi ra x(n+1) vaø laëp.

( ) x t n

(cid:14)

(cid:32)

(cid:16)

1 (cid:14)

(cid:11)

(cid:12)

x t '( ) n 1!

...

(cid:14)

(cid:16)

(cid:14)

1 (cid:14)

(cid:11)

(cid:12)

''( ) x t n 2!

3. Giaûi thuaät Runge-Kutta : xaáp xæ tieáp tuïc caùc ñaïo haøm baäc cao vaø thieát laäp coâng thöùc laëp cho giaù trò : x(tn+1) .

(cid:132) Thuaät toaùn khai trieån Taylor laø thuaät toaùn ñôn böôùc : giaù trò ñöôïc tìm khi chæ caàn bieát moät giaù trò cuûa böôùc tröôùc ñoù .

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

2

(cid:131) Giaûi thuaät Euler thuaän : ví duï 1

function vidu1=pplap()% Euler thuan % e(t) = 200i((cid:92)) + d(cid:92)/dt t0=0;tf=0.002;delta = 110^(-6); y0 = 0;% bien tu thong (cid:92) y=linspace(0,1,round(tf/delta)+1); y(1) = y0; for k=2:round(tf/delta)+1

et = 19sin(10^4(k-1)delta); iL = 50y(k-1)+ 410^7y(k-1)^3; y(k) = y(k-1)+delta(et-200iL);

end % Ve dang dong dien t=0:delta:deltaround(tf/delta); plot(t1000,y*1000); ylabel('Dang tu thong (mWb)'); % ==== End of Program ========

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

3

(cid:131) Giaûi thuaät Euler thuaän : ví duï 2

Xaùc ñònh daïng töø thoâng vaø doøng ñieän treân caûm phi

i

5 (cid:58)

e(t)

1 2 0

2 s in (1 0 0

9 0 )o

tuyeán . e t ( )

Bieát :

(cid:92)

(cid:32)

t(cid:83)

(cid:14)

_ +

Giaûi

1 V

(cid:132) Phöông trình traïng thaùi :

(cid:92) (Wb)

e t ( )

i 5 .

(cid:32)

1 (cid:14) (cid:16)

(cid:11) (cid:92)

(cid:12)

0,5625

d (cid:92) d t

0,45

-0,45

-0,225

i(A)

0,45

0,225 -0,45 -0,5625

(cid:132) Chöông trình MATLAB: Wb_dt=[-5.3375 -0.5625 -0.45 0 0.45 0.5625 5.3375] ; dong_dt=[-10 -0.45 -0.225 0 0.225 0.45 10];% A chuky = 1/50; t0 = 0; % gia tri ban dau tf = 75*chuky ;%gia tri cuoi cua thoi gian toi da delta = chuky/100; % buoc thoi gian tinh Wb=linspace(0,1,round(tf/delta)+1); % Tao mang tu thong iL=linspace(0,1,round(tf/delta)+1); % Tao mang luu dong % Gan gia tri so kien vao mang Wb(1) = 0; iL(1) = 0;

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

4

(cid:131) Giaûi thuaät Euler thuaän : ví duï 2 (tt)

for k=2:round(tf/delta)+1

% Tinh toan cac gia tri dong theo tu thong

gt_hdt= interp1(Wb_dt,dong_dt,Wb(k-1));

% Noi suy theo MATLAB

et=120sqrt(2)

sin(100pi(k-1)*delta+pi/2);

fn = et + 1 -5gt_hdt; Wb(k)=Wb(k-1)+deltafn; % Euler thuan iL(k)=interp1(Wb_dt,dong_dt,Wb(k));

end % Ve dang dong dien : chi khao sat

50 ms cuoi cung N = round(tf/delta); t=0:delta:delta*N; subplot(211); plot(t(N-250:N),Wb(N-250:N)); ylabel('Tu thong (Wb)'); subplot(212); plot(t(N-250:N),iL(N-250:N)); ylabel('Dong dien (A)');

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

5

(cid:131) So saùnh keát quaû do Simulink tính

(cid:132) Hình beân laø keát quaû do Simulink cuûa MATLAB tính toaùn vaø veõ .Coù theå tham khaûo taïi trang web cuûa MathWorks- >Documentation-> SimpowerSystems ; hay http.see.deis.unical.it/ita/ Kizilcay/pdf/chap2.pdf.

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

6 3. Giaûi thuaät Runge-Kutta

k 2

x n

k 1

k 2 2

k 2 3

k 4

x n

(cid:14)

(cid:14) (cid:32) 1

(cid:132) Gæai Thuaät R-K baäc 2 : (cid:12)

(cid:132) GT R-K baäc 4 (k1 tính nhö treân) : k (cid:12) 1

(cid:32)

k t , 1

(cid:32)

(cid:14)

h (cid:11) (cid:14) 2 x t ,n x

(cid:12) hk t , 1

(cid:32)

(cid:14)

(cid:11) (cid:11)

(cid:12)

(cid:183) (cid:184) (cid:185) (cid:183) h (cid:184) (cid:185) h

k 4

hk t 3, n

h (cid:11) 6 h 2 h 2 (cid:14)

(cid:32)

1 2 1 2 (cid:14)

k 33

(cid:14)

(cid:14) (cid:32) 1

(cid:167) f x (cid:168) n (cid:169) (cid:167) f x (cid:168) n (cid:169) f x (cid:11) n

k 3 k t , 2 n (cid:32) (cid:14) (cid:14)

(cid:12)

(cid:132) GT R-K baäc 3 (k1 nhö treân) : h 4

k t , 1

(cid:32)

(cid:14)

(cid:132) Giaûi thuaät R-K tính nhieàu caùc keát quaû trung gian vaø löu chuùng -> toán boä nhôù .

(cid:11) h 3 h 2 3

(cid:12) 1 3 2 3

(cid:167) (cid:168) (cid:169) (cid:167) f x (cid:168) n (cid:169)

(cid:183) (cid:184) (cid:185) (cid:183) (cid:184) (cid:185)

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

h k 3 k t , 2 (cid:32) (cid:14) (cid:14)

7

(cid:131) GT Runge-Kutta baäc 2 : Ví duï 1

10 (cid:80)F

+

Tìm aùp ra cuûa maïch keïp döông (ñænh aâm) bieát e(t) = 2sin(2(cid:83).103t) V.

+

100 K(cid:58)

-uC

Giaûi

u(t)

e(t)

-

)

i e u ( (cid:16) C

(cid:32)

(cid:14)

i(mA)

a) Thieát laäp pt traïng thaùi : Vieát KCL e u (cid:16) C R

du C dt

u(V)

0,4 0,6

= [-10 -5 -0.6 -0.4 0 5 10 ];% V

Löu yù: Ñaëc tuyeán i(e-uC) nhaän ñöôïc khi laáy ñoái xöùng qua goác toïa ñoä ñaëc tuyeán u-i cuûa diode ñaõ cho . b) Chöông trình MATLAB: % Thuat toan Runge-Kutta bac 2 % x(n+1) = x(n) + h/2(k1 + k2) % k1 = f(xn,tn) ; k2 = f(xn+hk1,tn+1) % Cho phuong trình phi tuyen - Dac tuyen dang bang-do thi % duc/dt = (e -uc) + 10^5hdt(e-uc) ; e = 2sin(2pift) % u = e - uc % Ve dang dac tuyen : dong = hdt(u) ap_dt dong_dt = [-40 -20 -0.004 0 0 0 0 ];% A

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

8

(cid:131) GT Runge-Kutta baäc 2 : Ví duï 1 (tt)

for k=2:round(tf/delta)+1 % Tinh toan cac gia tri ap

e = 2sin(2pif(k-1)*delta);% tn gt_hdt = interp1(ap_dt,dong_dt,

e - uc(k-1)); % Noi suy

k1 = e - uc(k-1) + 10^5gt_hdt; e = 2sin(2pif(k)delta);% tn+1 gt_hdt = interp1(ap_dt,dong_dt,

e-(uc(k-1)+ delta*k1)); % Noi suy

k2 = e - (uc(k-1) + deltak1) + 10^5gt_hdt; uc(k) = uc(k-1) + (delta/2)*(k1 + k2); % bien trang thai ap(k) = e - uc(k);% tin hieu ra

end % Ve dang dien ap ra t=0:delta:deltaround(tf/delta); plot(t1000,ap); % don vi ms

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

9

(cid:131) GT Runge-Kutta baäc 4 : Ví duï 1 )

i (cid:92) (

0, 5 (cid:92)(cid:32)

Khoùa K môû ra taïi t = 0, bieát ñaëc tuyeán cuûa caûm phi tuyeán : a) Vieát phöông trình traïng thaùi theo töø thoâng treân cuoän daây ? Xaùc ñònh giaù trò töø

thoâng taïi t = 0-.

b) Vieát chöông trình tính töø thoâng khi 0 < t < 50(ms) , vôùi böôùc tính h = 1(ms) theo thuaät toaùn Runge-Kutta baäc 4. Töø ñoù xaùc ñònh thôøi gian quaù ñoä cuûa maïch (thôøi gian khi ñoä thay ñoåi töø thoâng moùc voøng sau moät böôùc tính beù hôn 0,01% ).

Giaûi

250

a) Xaùc ñònh : Khi t < 0 : cuoän daây xem nhö ngaén maïch :

i (cid:32)

(cid:32)

5, 2 625

120

(cid:14)

0, 4(

)Wb

0, 08(

A )

i 1

(cid:32)

(cid:92)(cid:159) (cid:32) 0

0, 2 5, 2

1 5, 2 250 62,5 (cid:16)

2 (cid:92)

i 250 125 (cid:16)

(cid:32)

d (cid:92) dt

120 (cid:58)

5 (cid:58)

K

i 1

(cid:92)

t = 0

250 V

Khi t > 0 : phöông trình moâ taû maïch : b) Chöông trình MATLAB: function vidu11=pplap() t0 = 0 ; tf = 0.05; delta = 0.001; x0 = 0.4; x = linspace(0,1,round(tf/delta)+1); x(1) = x0;

0,2 (cid:58)

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

10 saiso = 1; flag = 0; k = 1; while (saiso > 0.0001) & (flag == 0)

(cid:131) GT Runge-Kutta baäc 4 : Ví duï 1 (tt) % So sanh va bao ket qua if flag == 1

k = k + 1; % Tinh gia tri bien theo Runge-

Kutta bac 4

disp('Vout qua tmax cho tinh toan'); t=0:delta:(k-1)*delta; plot(t,x(1:k)); % ve cac so lieu da tinh

else

k1 = 250-62.5x(k-1)^2; k2 = 250-62.5(x(k-

1)+(delta/2)*k1)^2;

k3 = 250-62.5*(x(k-

1)+(delta/2)*k2)^2;

k4 = 250-62.5(x(k-1)+(delta)k3)^2; x(k) = x(k-

1)+(delta/6)(k1+2k2+2*k3+k4);

disp('Thoi gian qua do la :'); disp((k-1)delta); disp('Gia tri tu thong khi do la :'); disp(x(k)); t=0:delta:(k-1)delta; plot(t,x(1:k)); % ve cac so lieu da tinh

% Tinh sai so saiso = abs((x(k)-x(k-1))/x(k)); % Kiem tra thoi gian tinh if k > round(tf/delta) flag = 1; end

end %----------------------------------------- Thoi gian qua do la : 0.0330 Gia tri tu thong khi do la : 1.9993

end % == end of main ======

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

11 9.5.3 Phöông phaùp ña thöùc hoùa

a x 0 n

(cid:132) Xaáp xæ haøm vaø ñaïo haøm baèng ña thöùc baäc k, coù ñaïi löôïng ôû böôùc tính (n+1) theo p giaù trò ôû p böôùc tröôùc ñoù (goïi laø giaûi thuaät ña böôùc) . Nghieäm toång quaùt cuûa giaûi thuaät coù daïng : a x p n p

(cid:32)

... (cid:14) (cid:14)

(cid:14)

x 1 n (cid:14)

t , n p n p (cid:16)

t 1 n (cid:14)

(cid:16)

h b f x , (cid:11) 1 1 n (cid:14) (cid:16)

b f x t , (cid:11) n n 0

(cid:12)

(cid:11) b f x p

(cid:170) (cid:172)

1. Giaûi thuaät hình thang :

(cid:12) (cid:186) (cid:188) (cid:132) Ñeå tính ñöôïc x(n+1) , phaûi coù caùc giaù trò ôû p böôùc tröôùc ñoù. Duøng giaûi thuaät ñôn böôùc (kieåu Taylor) ñeå khôûi ñoäng , sau ñoù duøng caùc giaûi thuaät ña böôùc ñeå deã kieåm soaùt sai soá. x (cid:32) (cid:14) n

(cid:14)

t 1 n (cid:14)

x 1 n (cid:14)

f x (cid:11) 1 n (cid:14)

(cid:12)

f x t , (cid:11) n n

(cid:12)

(cid:186) (cid:188)

(cid:170) (cid:172)

h 2

2. GT Adams-Bashforth (hieån do b-1 = 0) baäc 2 vaø 3 nhö sau :

(cid:12)

h , x (cid:32) (cid:14) n (cid:16) x 1 n (cid:14) t 1 n (cid:16) f x t , (cid:11) n n f x (cid:11) 1 n (cid:16) 1 2 (cid:186) (cid:187) (cid:188)

1 (cid:14)

1 (cid:16)

(cid:16)

(cid:12)

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

h t t , , x n x n (cid:32) (cid:14) (cid:16) (cid:14) f x t , (cid:11) n f x (cid:11) n f x (cid:11) n 3 2 23 12 16 12 5 12 (cid:170) (cid:171) (cid:172) (cid:170) (cid:171) (cid:172) (cid:186) (cid:187) (cid:188)

12

(cid:14)

(cid:32)

(cid:14)

(cid:11)

1 2

(cid:131) Giaûi Thuaät Adams-Moulton 3. GT Adams-Moulton (aån do b-1 (cid:122) 0) baäc 2 vaø 3 nhö sau : (cid:186) x (cid:187) (cid:188)

1 (cid:14)

1 (cid:16)

(cid:12)

(cid:170) (cid:171) (cid:172) 5 12

(cid:12) 8 12

(cid:12) 1 12

h t t , , x n x n (cid:32) (cid:14) (cid:14) (cid:16) f x (cid:11) n f x t , (cid:11) n f x (cid:11) n

500 (cid:58)

(cid:170) (cid:171) (cid:172) (cid:186) (cid:187) (cid:188)

(cid:132) Ví duï : Veõ daïng doøng vaø töø thoâng ? Giaûi

t=0

i 5 0 0 (

e t

( )

) (cid:92) (cid:14)

(cid:32)

i = 1,5(cid:92)5

d (cid:92) d t

(cid:92)

e(t)+ _

e(t) =1200sin(800(cid:83)t) V

function vidu16=pplap() % PP Adams-Bashforth bac 2 % Khoi dong theo Euler thuan % d(cid:92)/dt = 1200sin(800(cid:83)t)-750(cid:92)5 = f((cid:92),t) % x(n+1) = x(n) +h((3/2)fn - (1/2)fn_1) T = 1/400;t0 = 0; tf = 4*T; delta = T/50;

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

13

(cid:131) GT Adams-Bashforth baäc 2 : Ví duï 1

Wb = linspace(0,1,round(tf/delta)+1); iL = linspace(0,1,round(tf/delta)+1); Wb(1) = Wb0; iL(1) = 1.5Wb(1)^5; % gia tri tiep theo Euler thuan Wb(2) = Wb(1) + delta(-750*Wb(1)^5

+1200sin(0)); iL(2) = 1.5Wb(2)^5; for k=3:round(tf/delta)+1

fn = -750*Wb(k-1)^5

+1200sin((2pi/T)(k-1)delta);

fn_1 = -750*Wb(k-2)^5

+1200sin((2pi/T)(k-2)delta); Wb(k) = Wb(k-1) + delta((3/2)fn -

(1/2)*fn_1);

iL(k) = 1.5*Wb(k)^5;

end %=== end of main ========= t=0:delta:deltaround(tf/delta); subplot(211); plot(t1000,Wb); ylabel('Tu thong (Wb)'); subplot(212); plot(t*1000,iL); ylabel('Dong dien (A)'); % == End of Program ======

http://www.khvt.com

Chapter-94

(cid:11)(cid:38)(cid:12)(cid:47)(cid:114)(cid:3)(cid:48)(cid:76)(cid:81)(cid:75)(cid:3)(cid:38)(cid:1133)(cid:1249)(cid:81)(cid:74)

Mạch phi tuyến động và các PP số

Một mạch gọi là tuyến tính khi tuân theo định luật: Nếu y1(t) và y2(t) lần lượt là đáp ứng của hai nguồn kích thích độc lập với nhau x1(t) và x2(t), mạch là tuyến tính nếu và chỉ nếu đáp ứng đối với x(t)= k1x1(t) + k2x2(t)

Chủ đề:

9.5 Maïch phi tuyeán ñoäng vaø caùc PP soá

(cid:132) Caùc baøi toaùn coù 2 phaàn : 1. Vieát phöông trình traïng thaùi. 2. Söû duïng thuaät toaùn laëp. (cid:132) Caùc thuaät toaùn laëp :

Coù hai nhoùm :

1. Phöông phaùp khai trieån Taylor : tìm

chuoãi gaàn ñuùng haøm taïi x(tn).

2. Phöông phaùp ña thöùc hoùa : tìm ña thöùc gaàn ñuùng giaù trò x(tn) theo p giaù trò tröôùc ñoù .

9.5.1 Phöông trình traïng thaùi : (cid:132) Phöông trình traïng thaùi coù daïng : x’ = f(x,t) , x0 = x(t0) , hay neáu coù n bieán :

Tuï ñieän: Choïn uC hay qC. Cuoän daây: Choïn iL hay (cid:92)L. (cid:132) Phöông trình traïng thaùi thöôøng ñöôïc thieát laäp theo caùc phöông phaùp : doøng nhaùnh , theá nuùt hay maéc löôùi.

1

9.5.2 Phöông phaùp khai trieån Taylor

1. Giaûi thuaät Euler thuaän : x , n

(cid:12)

2. Giaûi thuaät Euler ngöôïc : n ,(

(cid:132) Khai trieån giaù trò taïi böôùc tính thöù (n+1) : x(tn+1) theo giaù trò taïi böôùc tính thöù (n) : x(tn) baèng chuoãi Taylor vaø choïn ñeán baäc p thích hôïp (sai soá h(p+1)) , vôùi h laø böôùc tính :

h f x . (cid:11) 1 n (cid:14) Giaûi ra x(n+1) vaø laëp.

3. Giaûi thuaät Runge-Kutta : xaáp xæ tieáp tuïc caùc ñaïo haøm baäc cao vaø thieát laäp coâng thöùc laëp cho giaù trò : x(tn+1) .

(cid:132) Thuaät toaùn khai trieån Taylor laø thuaät toaùn ñôn böôùc : giaù trò ñöôïc tìm khi chæ caàn bieát moät giaù trò cuûa böôùc tröôùc ñoù .

2

function vidu1=pplap()% Euler thuan % e(t) = 200*i((cid:92)) + d(cid:92)/dt t0=0;tf=0.002;delta = 1*10^(-6); y0 = 0;% bien tu thong (cid:92) y=linspace(0,1,round(tf/delta)+1); y(1) = y0; for k=2:round(tf/delta)+1

et = 19*sin(10^4*(k-1)*delta); iL = 50*y(k-1)+ 4*10^7*y(k-1)^3; y(k) = y(k-1)+delta*(et-200*iL);

end % Ve dang dong dien t=0:delta:delta*round(tf/delta); plot(t*1000,y*1000); ylabel('Dang tu thong (mWb)'); % ==== End of Program ========

3

Xaùc ñònh daïng töø thoâng vaø doøng ñieän treân caûm phi

i

tuyeán . e t ( )

Bieát :

Giaûi

4

for k=2:round(tf/delta)+1

% Tinh toan cac gia tri dong theo tu thong

gt_hdt= interp1(Wb_dt,dong_dt,Wb(k-1));

% Noi suy theo MATLAB

et=120*sqrt(2)*

sin(100*pi*(k-1)*delta+pi/2);

fn = et + 1 -5*gt_hdt; Wb(k)=Wb(k-1)+delta*fn; % Euler thuan iL(k)=interp1(Wb_dt,dong_dt,Wb(k));

end % Ve dang dong dien : chi khao sat

50 ms cuoi cung N = round(tf/delta); t=0:delta:delta*N; subplot(211); plot(t(N-250:N),Wb(N-250:N)); ylabel('Tu thong (Wb)'); subplot(212); plot(t(N-250:N),iL(N-250:N)); ylabel('Dong dien (A)');

5

(cid:132) Hình beân laø keát quaû do Simulink cuûa MATLAB tính toaùn vaø veõ .Coù theå tham khaûo taïi trang web cuûa MathWorks- >Documentation-> SimpowerSystems ; hay http.see.deis.unical.it/ita/ Kizilcay/pdf/chap2.pdf.

6

3. Giaûi thuaät Runge-Kutta

(cid:12)

(cid:132) Giaûi thuaät R-K tính nhieàu caùc keát quaû trung gian vaø löu chuùng -> toán boä nhôù .

7

+

Tìm aùp ra cuûa maïch keïp döông (ñænh aâm) bieát e(t) = 2sin(2(cid:83).103t) V.

+

-uC

Giaûi

-

a) Thieát laäp pt traïng thaùi : Vieát KCL e u (cid:16) C R

= [-10 -5 -0.6 -0.4 0 5 10 ];% V

8

for k=2:round(tf/delta)+1 % Tinh toan cac gia tri ap

e = 2*sin(2*pi*f*(k-1)*delta);% tn gt_hdt = interp1(ap_dt,dong_dt,

e - uc(k-1)); % Noi suy

k1 = e - uc(k-1) + 10^5*gt_hdt; e = 2*sin(2*pi*f*(k)*delta);% tn+1 gt_hdt = interp1(ap_dt,dong_dt,

e-(uc(k-1)+ delta*k1)); % Noi suy

k2 = e - (uc(k-1) + delta*k1) + 10^5*gt_hdt; uc(k) = uc(k-1) + (delta/2)*(k1 + k2); % bien trang thai ap(k) = e - uc(k);% tin hieu ra

end % Ve dang dien ap ra t=0:delta:delta*round(tf/delta); plot(t*1000,ap); % don vi ms

9

Khoùa K môû ra taïi t = 0, bieát ñaëc tuyeán cuûa caûm phi tuyeán : a) Vieát phöông trình traïng thaùi theo töø thoâng treân cuoän daây ? Xaùc ñònh giaù trò töø

thoâng taïi t = 0-.

Giaûi

a) Xaùc ñònh : Khi t < 0 : cuoän daây xem nhö ngaén maïch :

K

Khi t > 0 : phöông trình moâ taû maïch : b) Chöông trình MATLAB: function vidu11=pplap() t0 = 0 ; tf = 0.05; delta = 0.001; x0 = 0.4; x = linspace(0,1,round(tf/delta)+1); x(1) = x0;

10

saiso = 1; flag = 0; k = 1; while (saiso > 0.0001) & (flag == 0)

(cid:131) GT Runge-Kutta baäc 4 : Ví duï 1 (tt) % So sanh va bao ket qua if flag == 1

k = k + 1; % Tinh gia tri bien theo Runge-

Kutta bac 4

disp('Vout qua tmax cho tinh toan'); t=0:delta:(k-1)*delta; plot(t,x(1:k)); % ve cac so lieu da tinh

else

k1 = 250-62.5*x(k-1)^2; k2 = 250-62.5*(x(k-

1)+(delta/2)*k1)^2;

k3 = 250-62.5*(x(k-

function vidu1=pplap()% Euler thuan % e(t) = 200i((cid:92)) + d(cid:92)/dt t0=0;tf=0.002;delta = 110^(-6); y0 = 0;% bien tu thong (cid:92) y=linspace(0,1,round(tf/delta)+1); y(1) = y0; for k=2:round(tf/delta)+1

et = 19sin(10^4(k-1)delta); iL = 50y(k-1)+ 410^7y(k-1)^3; y(k) = y(k-1)+delta(et-200iL);

end % Ve dang dong dien t=0:delta:deltaround(tf/delta); plot(t1000,y*1000); ylabel('Dang tu thong (mWb)'); % ==== End of Program ========

et=120sqrt(2)

sin(100pi(k-1)*delta+pi/2);

fn = et + 1 -5gt_hdt; Wb(k)=Wb(k-1)+deltafn; % Euler thuan iL(k)=interp1(Wb_dt,dong_dt,Wb(k));

e = 2sin(2pif(k-1)*delta);% tn gt_hdt = interp1(ap_dt,dong_dt,

k1 = e - uc(k-1) + 10^5gt_hdt; e = 2sin(2pif(k)delta);% tn+1 gt_hdt = interp1(ap_dt,dong_dt,

k2 = e - (uc(k-1) + deltak1) + 10^5gt_hdt; uc(k) = uc(k-1) + (delta/2)*(k1 + k2); % bien trang thai ap(k) = e - uc(k);% tin hieu ra

end % Ve dang dien ap ra t=0:delta:deltaround(tf/delta); plot(t1000,ap); % don vi ms

k1 = 250-62.5x(k-1)^2; k2 = 250-62.5(x(k-

k4 = 250-62.5(x(k-1)+(delta)k3)^2; x(k) = x(k-

1)+(delta/6)(k1+2k2+2*k3+k4);

disp('Thoi gian qua do la :'); disp((k-1)delta); disp('Gia tri tu thong khi do la :'); disp(x(k)); t=0:delta:(k-1)delta; plot(t,x(1:k)); % ve cac so lieu da tinh

Wb = linspace(0,1,round(tf/delta)+1); iL = linspace(0,1,round(tf/delta)+1); Wb(1) = Wb0; iL(1) = 1.5Wb(1)^5; % gia tri tiep theo Euler thuan Wb(2) = Wb(1) + delta(-750*Wb(1)^5

+1200sin(0)); iL(2) = 1.5Wb(2)^5; for k=3:round(tf/delta)+1

+1200sin((2pi/T)(k-1)delta);

+1200sin((2pi/T)(k-2)delta); Wb(k) = Wb(k-1) + delta((3/2)fn -

end %=== end of main ========= t=0:delta:deltaround(tf/delta); subplot(211); plot(t1000,Wb); ylabel('Tu thong (Wb)'); subplot(212); plot(t*1000,iL); ylabel('Dong dien (A)'); % == End of Program ======