Hd Bt9 Gui Matlab

HÖÔÙNG DAÃN BT9 MATLAB GUI

(caùi nì daønh cho Trinh)

1. Thieát keá giao dieän: - Taïo giao dieän nhö sau:

Panel

Tag: NumDim String: 1 2 3 FontSize: 10

9) Baøi toaùn: Moâ phoûng chuyeån ñoäng nhieät vaø phaân boá Bolzmann. Tham soá: soá chieàu chuyeån ñoäng, soá laàn moâ phoûng. Output leân ñoà thò: Trình baøy ñoà thò keát quaû moâ phoûng.

Tag: NmbrFois String: FontSize: 10 Tag: sim String: Simulate

Trinh coù thaáy String cuûa NumDim coù 3 haøng khoâng? Ñeå taïo moät String nhieàu haøng thì taïi thuoäc tính String cuûa NumDim, thay vì nhaäp thaúng chöõ vaøo oâ keá beân, Trinh haõy baám vaøo nuùt keá beân vaø nhaäp giaù trò thaønh nhieàu doøng:

Ñoái vôùi caùc Static Text khaùc trong GUI, neáu muoán phaân thaønh nhieàu doøng, Trinh cuõng laøm nhö vaäy.

2. Laäp trình Nhaéc laïi giaûi thuaät moâ phoûng chuyeån ñoäng nhieät ñaõ giôùi thieäu trong taøi lieäu MonteCarlo (maø mình ñaõ ñöa cho Trinh): i) Cho moãi phaân töû moät vaän toác ban ñaàu baát kyø; ii) Taïi moãi böôùc laëp choïn moät phaân töû khí xaùc ñònh, cho vaän toác cuûa noù thay ñoåi ngaãu nhieân moät löôïng nhoû; iii) Tính ñoä bieán thieân ñoäng naêng cuûa heä thoáng do thay ñoåi ngaãu

e nhieân vöøa roài; Neáu thì traïng thaùi môùi noùi treân ñöôïc giöõ laïi laøm traïng thaùi ban

e  0 ñaàu cuûa böôùc laëp keá tieáp; Neáu thì traïng thaùi môùi chæ ñöôïc chaáp thuaän vôùi xaùc suaát laø

e  0 . Ñeå thöïc hieän ñieàu naøy treân maùy tính ta choïn moät soá ngaãu

e



e kT

nhieân r trong khoaûng ]0,1[, neáu

thì traïng thaùi môùi ñöôïc giöõ laïi,

re



e kT

neáu khoâng thì traïng thaùi cuõ seõ ñöôïc duøng ñeå thöïc hieän böôùc laëp keá tieáp. Sau nhieàu böôùc laëp nhö treân, heä seõ tieán daàn ñeán phaân boá Maxwell.

Sau khi löu GUI, Trinh haõy tìm ñeán haøm sim_Callback trong M-file. Nhôù laïi trong worksheet Maple, ngöôøi ta coù söû duïng pheùp tính modulo (chia laáy phaàn dö) k:=(k+1) mod NMol;

nhöng trong Matlab laïi khoâng coù toaùn töû thöïc hieän pheùp tính naøy neân ta seõ vieát theâm moät function cho noù. Trinh haõy cheøn ñoaïn maõ sau vaøo ngay tröôùc haøm sim_Callback: % --- Function for calculating modulo operating function r = modulo(p,q) if p == q r = 0; else r = p - floor(p/q); % floor: lấy phần nguyên về phía giảm end

Ví duï neáu x = modulo(17, 3) thì x = 2 vì 17 = 3.5 + 2 Ngoaøi ra, trong whorksheet Maple, ngöôøi ta coù ñònh nghóa moät haøm ñeå tính toång naêng löôïng cuûa caùc haït nhö sau: #--------------------------------#Calculates system kinetic energy #--------------------------------energy:=proc(v,NMol) local e,i; e:=0;

v: array chöùa caùc giaù trò vaän toác caùc haït trong heä. NMol: Soá haït trong heä (Number of Molecule)

Caùch tính ôû treân laø

nhöng vì ngöôøi ta ñaõ giaû thieát mi =

e

1 mi vi2  2

1 (nhö phaàn ñaàu taøi lieäu MonteCarlo ñaõ noùi) neân chæ coøn .

e

1 vi2  2

Vieát laïi haøm naøy beân Matlab laø (ñoaïn maõ naøy cuõng ñöôïc cheøn vaøo phía tröôùc haøm sim_Callback trong M-file): % --- Function for calculating system kinetic energy. function e = Energy(V, m) % V: The array of velocities % m: The mass of each partical

vôùi V laø ma traän 3 coät chöùa caùc giaù trò vaän toác cuûa caùc haït trong heä (moãi haøng ma traän laø daønh cho moät haït), m laø khoái löôïng moãi haït. Maëc duø coù maët m nhöng trong sim_Callback ta cuõng seõ cho m = 1 ñeå vieäc moâ phoûng, tính toaùn ñöôïc ñôn giaûn (nghó coi, phaân töû maø cho m = 1kg thì quaù ñaùng thaät, nhöng chöa heát ñaâu, vaøo chöông trình chính roài môùi thaáy coøn maáy caùi kinh khuûng nöõa: haèng soá Boltzmann kB = 1,38065.10-23 thì ñaåy leân baèng 1, coøn nhieät ñoä T thì chæ coù 1oK (–272oC) ). Haõy xem laïi haøm Energy(V, m) ôû treân, taïi sao laïi laø W = V.^2 chöù khoâng phaûi V^2? Haõy luoân nhôù raèng döõ lieäu trong Matlab luoân ôû daïng ma traän (hoaëc taäp hôïp caùc ma traän). Moãi phaàn töû cuûa ma traän V laø moät giaù trò vaän toác vaø ta caàn tính bình phöông vaän toác, töùc laø bình phöông cuûa töøng phaàn töû ma traän (chöù khoâng phaûi bình phöông ma traän) neân ta phaûi vieát V.^2 (V^2 laø pheùp nhaân hai ma traän V*V). Khaùc vôùi chöông trình moâ phoûng trong worksheet Maple chæ xeùt chuyeån ñoäng 1 chieàu (theo truïc x), baøi taäp cuûa Trinh laïi phaûi moâ phoûng chuyeån ñoäng 1 – 3 chieàu neân vôùi moãi haït (moãi haøng cuûa ma traän), ta phaûi löu 3 giaù trò vaän toác laø vx, vy, vz. Ñieàu ñoù giaûi thích taïi sao V coù 3 coät. Doøng maõ e = m*sum(W(:))/2; nghóa laø gì? Ta haõy xem:

e  m vi2 / 2  m (vix2  viy2  viz2 ) / 2 i

i

cho neân sum(W(:)) tính toång caùc phaàn töû cuûa ma traän W. Chuù yù raèng sum(W) chæ traû veà moät ma traän haøng maø phaàn töû thöù i laø toång caùc phaàn töû treân coät thöù i cuûa ma traän W. Vaäy laø ñaõ xong nhöõng “haøm ngoaøi leà”, giôø ta ñi vaøo haøm chính, haøm sim_Callback. Ñaàu tieân laø böôùc chuaån bò cho tính toaùn: % --- Executes on button press in sim. function sim_Callback(hObject, eventdata, handles) % hObject

handle to sim (see GCBO)

% eventdata

reserved - to be defined in a future version of MATLAB

% handles

structure with handles and user data (see GUIDATA)

clc NMol = 10; T = 1; kB = 1; Theo worksheet Maple, mình cuõng cho T = kB = m = 1 ñeå tính toaùn ñôn

giaûn (maø vaãn khoâng aûnh höôûng ñeán quy luaät phaân boá Boltzmann m = 1; cuûa chyeån ñoäng nhieät). Cho soá haït laø Nmol = 10. Doøng NmbrFois = str2double(get(handles.NmbrFois, 'String')); chaéc cuõng quen thuoäc vôùi Trinh roài, phaûi khoâng? Noù laáy soá laàn moâ phoûng maø ngöôøi duøng nhaäp vaøo oâ input text NmbrFois treân GUI roài gaùn cho ma traän NmbrFois. Coøn doøng NumDim = get(handles.NumDim, 'Value'); laáy soá chieàu chuyeån ñoäng töø popup menu NumDim. Popup menu naøy coù 3 löïa choïn, cöù moãi caùi ñöôïc choïn thì Value seõ ñöôïc gaùn giaù trò laø soá thöù töï cuûa löïa choïn ñoù. Ví duï neáu popup menu coù caùc löïa choïn ‘heo’, ‘boø’, ‘gaø’ thì Value = 2 khi ‘boø’ ñöôïc choïn. Theo yù cuûa giaûi thuaät ôû treân, ta seõ duøng ma traän V0 ñeå moâ taû (baèng vaän toác) traïng thaùi ban ñaàu vaø V1 moâ taû traïng thaùi ñaït ñöôïc sau khi thay ñoåi ngaãu nhieân vaän toác cuûa 1 haït naøo ñoù. Tuøy theo ta muoán moâ phoûng vôùi bao nhieâu chieàu chuyeån ñoäng maø V0, V1 seõ coù baáy nhieâu coät (thaønh phaàn x, y, z cuûa vaän toác). zeros(NMol, NumDim) taïo ra ma traän NMol haøng, NumDim coät chöùa toaøn soá 0. Cöù moãi laàn ñaït ñeán traïng thaùi môùi, traïng thaùi ñoù chæ ñöôïc chaáp nhaän neáu noù thoûa ñk ôû muïc iii) neâu trong giaûi thuaät, ta seõ

duøng bieán j ñeå ñeám soá traïng thaùi chaáp nhaän ñc trong suoát quaù trình moâ phoûng. Theo muïc ii) cuûa giaûi thuaät thì ta seõ choïn 1 haït, nhöng laøm vaäy thì sôï (Nmol – 1) haït coøn laïi ghen tò, thoâi thì cöù cho haït naøo cuõng ñöôïc choïn ñeå coâng baèng. Nhö vaäy, laàn moâ phoûng thöù nhaát ta choïn haït 1, laàn 2 choïn haït 2, laàn 3 laø haït 3 vaø cöù theá ñeán heát roài quay voøng. Ñieàu naøy gôïi cho ta söû duïng pheùp chia laáy dö (modulo) ñeå dieãn taû böôùc ñi vöøa roài. Ñieàu naøy lyù giaûi söï coù maët cuûa bieán k vaø haøm modulo() maø ta ñaõ xaây döïng (thaày Chieán noùi roài, caùi gì toàn taïi cuõng coù lyù do cuûa noù ). Baây giôø ta ñi tieáp coâng ñoaïn tính toaùn: if ~isnan(NmbrFois) for i = 1:NmbrFois; k = modulo(k, NMol) + 1; % Phep chia lay phan du V1(k, :) = V0(k, :) + randn(1, NumDim); e0 = Energy(V0, m); % Động năng của hệ ở trạng thái ban đầu e1 = Energy(V1, m); % Động năng của hệ ở trạng thái mới de = e1 - e0; aRan = rand(1); ExP = exp(-de/(kB*T));

Doøng if ~isnan(NmbrFois) noùi raèng vieäc tính toaùn seõ ñöôïc thöïc if de < 0 || aRan < ExP hieän neáu NmbrFois khoâng phaûi NaN (Not a Number). Haõy xem laïi doøng NmbrFois = str2double(get(handles.NmbrFois, 'String'));. Neáâu ai ñoù lôõ queân, thay vì nhaäp moät soá vaøo input text treân GUI haén ta laïi nhaäp moät chuoãi kyù töï nhö Doraemon chaúng haïn, thì haøm str2double seõ traû veà giaù trò NaN (Not a Number), khoâng phuø hôïp cho tính toaùn neân ta phaûi duøng leänh if treân ñeå lôõ gaëp thì khoûi laøm nöõa. Thöïc ra thì neáu khoâng coù if, nhöõng leänh beân döôùi vaãn cöù tieán haønh cho duø coù gaëp NaN ñi nöõa (ñaây laø lí do khieán paø Haûi löôøi bieáng töø choái tieáp thu khi mình chæ caùi naøy cho “paû”, thieät laø…, neáu “paû” chòu nghe ñeå truyeàn laïi cho Trinh thì mình ñaâu caàn phaûi noùi laïi ôû ñaây nöõa ). ÖØ thì maáy baøi khaùc coù theå boû qua ñöôïc böôùc kieåm tra naøy, nhöng vôùi baøi taäp naøy cuûa Trinh thì khaùc, neáu cöù ñeå caùi NaN “nhôûn nhô ngoaøi voøng phaùp luaät” thì noù seõ laøm voøng laëp beân döôùi chaïy maõi chaïy maõi khoâng döøng ñöôïc, keát quaû laø Crash!

Doøng k = modulo(k, NMol) + 1; keát hôïp vôùi voøng laëp for i = 1:NmbrFois; seõ laøm cho k chaïy töø tuaàn hoaøn 1  NMol, töùc laø ñi töø haït 1 ñeán haït thöù NMol. Doøng V1(k, :) = V0(k, :) + randn(1, NumDim); dieãn taû: Töø traïng thaùi ban ñaàu V0, cho caùc thaønh phaàn vaän toác cuûa haït thöù k thay ñoåi ngaãu nhieân ñöôïc traïng thaùi môùi V1. [ : ] nghóa laø laáy moïi phaàn töû trong daõy. V(k, :) coù nghóa: laáy moïi coät treân haøng thöù k. Haøm randn(m, n) taïo ma traän m x n chöùa caùc soá ñöôïc phaùt sinh ngaãu nhieân trong khoaûng (-1; 1), khaùc vôùi haøm rand(m, n) laáy soá ngaãu nhieân trong khoaûng (0; 1). Ñk de < 0 || aRan < exp(-de/(kB*T)) laø ñieàu kieän ñöôïc noùi trong muïc iii) cuûa giaûi thuaät. e(j) = e1; löu laïi naêng löôïng e1 cuûa traïng thaùi môùi naøy ñeå tí nöõa

veõ ñoà thò keát quaû moâ phoûng. Nhaéc laïi: j duøng ñeå ñeám soá traïng thaùi ñöôïc chaáp nhaän. Vaäy laø xong phaàn moâ phoûng – tính toaùn roài ñoù. Baây giôø ta seõ xuaát keát quaû moâ phoûng leân ñoà thò. Ta seõ veõ bieåu ñoà coät phaân boá naêng löôïng vaø chöùng minh vaän toác trung bình moãi haït khi heä tieán daàn ñeán phaân boá Boltzmann laø NumDim*kBT/2. Ta chia truïc hoaønh (E ~ kBT) thaønh nhöõng khoaûng naêng löôïng (E, E+dE) vaø tính soá haït coù naêng löôïng naèm trong vuøng ñoù. e1

e0

e

C B A D

nints = 20; de = kB*T/nints; for i = 1:3*nints n(i, 1) = 0;

soá haït Haøm fill(X, Y,Ñeám ‘colour’) ñeå coù naêng löôïng veõ ña giaùc naèm trong (e, Vôùi X = [xA xB xCkhoaûng xD] Y = [yA yB yC yD]

n(i)/j: Soá haït moãi traïng thaùi (~ E(i, 2) = E(i, 1) + de; r : red (ñoû), b : blue, : green dieängtích coät i). n(i)/j/de: Chieàu cao for k = 1 : j coät. E(i, 1) = (i - 1)*de; colour: maøu

t = e(k)/NMol;

if t >= E(i, 1) && t < E(i, 2) n(i, 1) = n(i) + 1;

Cöù moãi voøng laëp laø veõ moät ñoà thò, bình thöôøng thì ñoà thò môùi seõ thay theá ñoà thò cuõ. Ñeå giöõ laïi ñoà thò, ta duøng leänh hold on. Tieáp theo, ta ñaët teân cho truïc hoaønh vaø truïc tung. Ngoaøi ra, ta vieát theâm thoâng baùo loãi neáu ngöôøi duøng nhaäp soá laàn moâ phoûng khoâng hôïp leä: xlabel('E') ylabel('f(E)') else errordlg('Oops, you''ve entered wrong value!')

errordlg nghóa laø “error dialogue”. Trong caâu thoâng baùo “Oops, you've

entered wrong value!”, ta coù duøng daáu nhaùy maø daáu naøy laïi ñöôïc Matlab duøng ñeå bao 2 ñaàu chuoãi kyù töï neân ñeå traùnh Matlab nhaàm töôûng chuoãi cuûa mình keát thuùc taïi “you”, ta haõy theâm 1 daáu nhaùy

nöõa laøm kyù töï thoaùt (töông töï daáu “\” trong C). Luùc naøy '' töông ñöông vôùi ' thoâng thöôøng. Vaäy laø chöông trình cuûa mình ñaõ xong roài ñoù. Cho Trinh xem tröôùc keát quaû moâ phoûng vôùi NmbrFoise = 5000:

Soá chieàu chuyeån ñoäng : 1 => ñoäng naêng trung bình: kT/2

Soá chieàu chuyeån ñoäng : 2 => ñoäng naêng trung bình: 2kT/2

Soá chieàu chuyeån ñoäng : 3 => ñoäng naêng trung bình: 3kT/2

Cuøng xem laïi phaàn laäp trình ta ñaõ laøm treân M-file:

% --- Function for calculating modulo operating function r = modulo(p,q) if p == q r = 0; else r = p - floor(p/q); end

% --- Function for calculating system kinetic energy. function e = Energy(V, m) % V: The array of velocities % m: The mass of each partical W = V.^2; e = m*sum(W(:))/2;

% --- Executes on button press in sim. function sim_Callback(hObject, eventdata, handles) % hObject

handle to sim (see GCBO)

% eventdata

reserved - to be defined in a future version of MATLAB

% handles

structure with handles and user data (see GUIDATA)

clc NMol = 10; T = 1; kB = 1; m = 1; NmbrFois = str2double(get(handles.NmbrFois, 'String')); NumDim = get(handles.NumDim, 'Value'); V0 = zeros(NMol, NumDim); V1 = zeros(NMol, NumDim); k = 0; j = 0; if ~isnan(NmbrFois) for i = 1:NmbrFois; k = modulo(k, NMol) + 1; % Phep chia lay phan du V1(k, :) = V0(k, :) + randn(1, NumDim); e0 = Energy(V0, m); % Động năng của hệ ở trạng thái ban đầu e1 = Energy(V1, m); % Động năng của hệ ở trạng thái mới de = e1 - e0; aRan = rand(1); ExP = exp(-de/(kB*T)); if de < 0 || aRan < ExP j = j + 1; V0(k, :) = V1(k, :);

e(j) = e1; else V1(k, :) = V0(k, :); end end nints = 20; de = kB*T/nints; for i = 1:3*nints n(i, 1) = 0; E(i, 1) = (i - 1)*de; E(i, 2) = E(i, 1) + de; for k = 1 : j t = e(k)/NMol; if t >= E(i, 1) && t < E(i, 2) n(i, 1) = n(i) + 1; end end n(i) = n(i)/j/de; fill([E(i, 1) E(i, 2) E(i, 2) E(i, 1)], [0 0 n(i) n(i)], 'c') hold on end hold off xlabel('E') ylabel('f(E)') else errordlg('Oops, you''ve entered wrong value!') end

Theá laø xong, phuø… N.H.Q Khoâng hieåu choã naøo, cöù maïnh daïn hoûi nheù. Neáu ngaïi thì nhôø ñöùa naøo ñoù chuyeån lôøi ñeán cho mình, OK? (ÔØ, taïi sao laïi ngaïi cô chöù? Ngöôøi ta coù aên thòt ñaâu…

)