2. 11. 2011
KTE / PPEL Počítačová podpora v elektrotechnice Ing. Lenka Šroubová, Ph.D. email:
[email protected] http://home.zcu.cz/~lsroubov Polynomiální regrese polyfit(x, y, st) – proloží množinu bodů o souřadnicích obsažených ve vektorech x a y polynomem stupně st. Výstupem jsou koeficienty polynomu zvoleného stupně st. Např.: Při měření byla získána data uvedená v tabulce. t[s] 1 2 3 4 5 u[V] 5.5 43.1 100.2 190.7 218.4 t = [1:5]; u = [5.5, 43.1, 100.2, 190.7, 218.4] plot(t,u,'o') xlabel('t[s]') ylabel('u[V]') title('Data z mereni') show Body na příslušných souřadnicích zobrazeny jako kolečka.
Proložíme naměřené hodnoty polynomem 3. řádu a vykreslíme jeho průběh společně s naměřenými hodnotami do grafu. Nejprve vypočteme regresní polynom: regr = polyfit(t,u,3) regr = -6.8583 62.6964 -110.3452 61.5800 tj.: regr(x) = – 6.8583 x3 + 62.6964 x2 -110.3452 x1 + 61.5800 x0 = – 6.8583 x3 + 62.6964 x2 -110.3452 x + 61.5800 Regresní křivku (zjištěný polynom) přidáme do grafu t t = 1 2 3 4 5 Pro zobrazení polynomu musíme zjemnit x-ovou osu (na ní zobrazeno t), vytvoříme vektor rt: rt = [0:0.05:6]; Vektor rt je z intervalu od 0 do 6 => dali jsme kousek i okolo. Pak vypočteme hodnoty polynomu pro všechna rt: hp = polyval(regr,rt); plot(t,u,'o',rt, hp, 'r') show
Prokládané body i křivka (polynom) jsou v jednom grafu, polynom zobrazen červeně, aby byl dobře vidět ;-)
Křivka není úplně v pořádku, zvláště na krajích intervalu. Zkusíme proložit naměřené hodnoty polynomem 4. řádu a vykreslíme jeho průběh společně s naměřenými hodnotami do grafu. regr_moc = polyfit(t,u,4) regr_moc = -4.5875 48.1917 -164.7125 263.2083 -136.6000 hp_moc = polyval(regr_moc,rt); plot(rt, hp_moc, 'g') show
Prokládané body i křivky (polynomy) jsou v jednom grafu, polynom 3. st. zobrazen červeně, polynom 4. st. zobrazen zeleně. Dostali jsme nevhodnou regresní křivku, nesmyslnou. Ani ta červená není úplně v pořádku…Zkusíme tedy ještě polynom ještě vyššího stupně. regr_jeste_vic = polyfit(t,u,6); hp_jeste_vic = polyval(regr_jeste_vic,rt); plot(t,u,'o',rt, hp, 'r',rt, hp_moc, 'g',... rt, hp_jeste_vic, 'k') show
Polynom 6. st. zobrazen černě, je to ještě horší, proto zkusíme proložit naměřené hodnoty polynomem 2. stupně: regr_mene = polyfit(t,u,2) regr_mene = 0.9714 51.5114
-53.6400
hp_mene = polyval(regr_mene,rt); plot(t, u,'o',... rt, hp, 'r',... rt, hp_moc, 'g',... rt, hp_jeste_vic, 'k',... rt, hp_mene, 'm')
show
Nejlépe by zde vyhověl polynom 2. stupně. Zkusíme ještě přímku (polynom 1. stupně): regr_min = polyfit(t,u,1) regr_min = 57.340 -60.440 hp_min = polyval(regr_min,rt); regr_min = polyfit(t,u,1) hp_min = polyval(regr_min,rt); plot(t,u,'o', rt,hp_mene,'m', rt,hp_min,'c') legend('body','polynom 2.st.','polynom 1.st.','Location','SouthEast') show Polynom 1. stupně y(x) = 57.34 x – 60.44 (rovnice přímky) Polynom 2. stupně y(x) = 0.9714 x2 + 51.5114 x – 53.6400 (rovnice paraboly) Polynom 1. stupně a polynom 2. stupně na intervalu od 0 do 6 mají velmi podobný průběh a jsou vhodné pro proložení těchto naměřených dat z tabulky.
Příklady – viz http://edison.fel.zcu.cz – publikované sešity Lze použít i další typy regresních křivek, nejen výše uvedené polynomy Např.: spline Použitelné funkce: interp1 – 1Dimenzionální interpolace interp2 – 2D interpolace interp3 – 3D interpolace interpn – nDimenzionální interpolace interpft – interpolace s využitím FFT Např. 1D interpolace pomocí spline fukcí spl_int = interp1(t,u,rt,'spline'); plot(t,u,'o',rt,spl_int, 'g')
Příklad: vlastní funkce pro polynomiální interpolaci Vstupní data: x-ové a y-ové souřadnice bodů (2 vektory), stupeň polynomu, kterým budou body proloženy. Nejprve provedeme jemnější dělení osy x pro zobrazení křivky (od minima z vektoru x do maxima z vektoru x s krokem vypočteným jako (max(x)-min(x))/100), poté získáme koeficienty polynomu ntého stupně pro polynomiální interpolaci příkazem polyfit a vyčíslíme hodnoty polynomu pro všechny body z jemnějšího dělení osy x příkazem polyval. Nakonec vykreslíme graf. function prolozeni_bodu(x,y,n) xp=[min(x):((max(x)-min(x))/100):max(x)]; yn=polyval(polyfit(x,y,n),xp); plot(x,y,'*',xp,yn,'m') show endfunction Volání funkce: prolozeni_bodu([1:7],[1,2,2,4,3,4,5],3)
Výstup – textový (bez možnosti formátování textu) disp – výstup bez názvu proměnné, nemá možnost formátovat text promenna=7; disp(promenna); 7 disp([1:5]) 1 2
3
4
5
disp('textovy retezec 1'); textovy retezec 1 disp(['textovy ','retezec ', '2']) textovy retezec 2 Pozn.: num2str – převod čísel na řetězec, např. num2str(rand) ans = 0.81472 whos Name Size Bytes ans 1x7 14
Class char
Využití disp a num2str: x=5; disp(['Zvetsime-li cislo ', num2str(x), ' o jednicku, dostaneme ', num2str(x+1)]); Výstup na obrazovku: Zvetsime-li cislo 5 o jednicku, dostaneme 6 Vlastní uživatelské funkce – opakování a pokračování Příklad – výpočet obsahu a obvodu kruhu function [S,o]=obsah_obvod(r) S=pi.*(r.^2); o=2.*pi.*r; endfunction Volání funkce [obs,obv]=obsah_obvod(10) obs = 314.16 obv = 62.832 Některé z uvedených funkcí v MATLABu – m-file File –> New –> M-file (editor) function [S,o]=obsah_obvod(r) S=pi.*(r.^2); o=2.*pi.*r; Soubor je nutné uložit. Save as: obsah_obvod.m Volání funkce v příkazovém okně Command Window [obs,obv]=obsah_obvod(10) obs = 314.16 obv = 62.832
whos Name Size Attributes obs 1x1 obv 1x1
Bytes
Class 8 8
double double
r – lokální proměnná (s ukončením funkce zaniká), výstupní proměnné zůstanou zachovány Funkce bez výstupu function obsah_obvod2(r) S=pi.*(r.^2); o=2.*pi.*r; disp('Obsah') disp(S) disp('Obvod') disp(o) Volání funkce obsah_obvod2(5) Obsah 78.5398 Obvod 31.4159 r, S, o – lokální proměnné (s ukončením funkce zanikají), výstupní proměnná není whos => po příkazu whos se nevypíše nic
Relační operátory == porovnání na rovnost (je rovno) ~= porovnání na nerovnost (není rovno) <, > je menší, je větší <=, >= je menší nebo rovno, je větší nebo rovno ~ negace (not) Příklady – viz http://edison.fel.zcu.cz – publikované sešity Logické operátory & – a zároveň (and) | – nebo (or) xor – „exkluzivní“ nebo ~ – negace (not) (3<5)&(4<6) ans = 1 (3<5)&(4>6) ans = 0 (3>5)&(4>6) ans = 0 Lze psát i takto: and((3>5),(4<6)) ans = 0 (3<5)|(4<6) ans = 1 (3<5)|(4>6) ans = 1 (3>5)|(4>6) ans = 0
– pravda a zároveň pravda – pravda a zároveň nepravda – nepravda a zároveň nepravda
– nepravda a zároveň pravda
– pravda nebo pravda – pravda nebo nepravda – nepravda nebo nepravda
Lze psát i takto: or((3>5),(4<6)) ans = 1
– nepravda nebo pravda
xor((3<5),(4<6)) ans = 0 xor((3>5),(4<6)) ans = 1 xor((3>5),(4>6)) ans = 0
– pravda nebo pravda exkluzivně
~0 ans =1 ~1 ans =0 ~5 ans =0 ~(3<5) ans =0 not(3<5) ans =0 Tedy pak: (~(3<5))&(4<6) ans = 0
– negace 0 (nepravdy) je 1 (pravda)
– pravda nebo pravda exkluzivně – pravda nebo pravda exkluzivně
– negace 1 (pravdy) je 0 (nepravda) – negace 5 (pravdy) je 0 (nepravda) – negace pravdy je nepravda – lze psát i takto, nepravda
– nepravda a zároveň pravda
Řízení běhu výpočtu - řídící příkazy: Podmíněný příkaz if logický_výraz příkaz; příkaz; příkaz; end; Středníky nejsou povinné. nebo plný tvar: if logický_výraz příkazy_je_li_podmínka_pravdivá; else příkazy_není_li_podmínka_pravdivá; end; Středníky nejsou povinné. nebo if logický_výraz příkazy1; elseif jiný_logický_výraz příkazy2; else příkazy3; end; Středníky nejsou povinné. Pozor else nebo elseif se vztahuje vždy k nejbližšímu if nad ním. (samozřejmě pokud se před ním nevyskytne end)
Příklad: Funkce pro porovnaní čísla s určitou mezí – vstupní parametry: číslo, mez testu,se kterou se bude číslo porovnávat function porovnani(cislo, mezTestu) if (cislo < mezTestu) disp('cislo je mensi nez zadana mez') elseif (cislo == mezTestu) disp('cislo je rovno zadane mezi') else disp('cislo je vetsi nez zadana mez') end Ve funkci – relační operátor == (rovná se ve smyslu porovnání) Volání funkce: porovnani(3, 2) cislo je vetsi nez zadana mez nebo porovnani(5, 5) cislo je rovno zadane mezi Příklad: Funkce, která vykreslí zadaný polynom p v zadaném rozsahu, najde a označí v grafu maximum a minimum v daném rozsahu a pokud jsou všechny kořeny reálné, vykreslí je také do grafu. Takže do příkladu, který byl uveden dříve, přidáme výpočet a zobrazení kořenů (do grafu jen reálné kořeny) function k=kresliPolynom(p,odkudX,kamX,prvkuX) x = linspace(odkudX,kamX,prvkuX); hodnotyPolynomu = polyval(p,x); plot(x,hodnotyPolynomu) title('Prubeh polynomu') xlabel('x')
ylabel('y') hold on [Max,indexMax]=max(hodnotyPolynomu); plot(x(indexMax),Max,'or') [Min,indexMin]=min(hodnotyPolynomu); plot(x(indexMin),Min,'ok') k = roots(p); if (isreal(k)==1) plot(k, polyval(p,k),'og') else disp('koreny polynomu nejsou realne') end hold off Volání funkce – např.: p(x) = x2 + 5x + 6 x = kresliPolynom([1,5,6],-4,-1,100) x = -3.0000 -2.0000
Příklady – viz http://home.zcu.cz/~lsroubov/PPEL nebo http://portal.zcu.cz
