| |
12. Soustavy lineárních rovnic |
|
|
| |
- (rozšíření, lze zařadit kdykoli
poté, co se seznámíme s funkcemi)
- Motto: Rozdíl mezi matematikou a technikou je v tom,
že matematik trvá na existenci řešení problému a jejich
počtu, zatímco technik by chtěl za všech okolností vidět
jen jedno rozumné řešení.
První následek: Chce-li být matematik spokojen, používá
Matlabovské funkce det, rank a inv.
Druhý následek: Ostatní buď používají Matlabovskou funkci
pinv místo inv, nebo jsou nešťastni.
- Před započetím práce vytvořte několik regulárních,
přeurčených a nedourčených soustav lineárních rovnic a
nastudujte Matlabovské funkce rank, det, inv, pinv a
hilb.
|
|
|
| |
Příklad 12A.1: |
|
| |
Vytvořte funkci pro zjištění počtu řešení
soustavy lineárních rovnic.
Instrukce: Není důležité vyřešit, ale určit počet řešení pomocí
hodností matic. Krokujte pro různé případy a sledujte hodnosti
matic. |
|
| |
Řešení: |
|
| |
| function [r]=Kolik(A,b) |
| % Urceni poctu reseni soustavy linearnich rovnic Ax=b |
| % [r]=Kolik(A,b); |
| % r ... pocet reseni |
| % A ... matice soustavy (m,n) |
| % b ... sloupcovy vektor pravych stran (m,1) |
| r=0; |
% pocatecni pesimismus |
| ma=size(A,1); |
% urceni poctu radku matice A |
| na=size(A,2); |
% urceni poctu sloupcu matice A |
| mb=size(b,1); |
% urceni delky vektoru b |
| nb=size(b,2); |
% pocet sloupcu vektoru b |
| if ma~=mb | nb~=1 |
% nevhodne rozmery matice nebo vektoru? |
| return; |
% kolik je r? |
| end |
| ra=rank(A); |
% hodnost matice A |
| rab=rank([A b]); |
% hodnost rozsirene matice A|b |
| if ra~=rab |
% nestejne hodnosti? |
| return; |
% kolik je r? |
| end |
| if ra==na |
% hodnost rovna poctu neznamych? |
| r=1; |
% ma prave jedno reseni |
| else |
| r=inf; |
% nema reseni |
| end |
|
|
| |
| |
Příklad 12A.2: |
|
| |
Vytvořte funkci pro testování inverzní
matice.
Instrukce: Inverzní matici vynásobíme zleva původní maticí a
porovnáme s jednotkovou maticí. Vlivem zaokrouhlovacích chyb
není jejich rozdíl vždy nulová matice. Při jaké velikosti Hilbertovy
matice bude chyba větší než jedna? |
|
| |
Řešení: |
|
| |
| function [test]=Inverze(A) |
| % Numericky test inverze matice |
| % [test]=Inverze(A); |
| % test ... maximalni prvek matice A * inv(A) - I |
| % A ... matice (n,n) |
| % I ... jednotkova matice (n,n) |
| m=size(A,1); |
% urceni poctu radku matice A |
| n=size(A,2); |
% urceni poctu sloupcu matice A |
| if m~=n | det(A)==0 |
% nevhodna matice |
| test=inf; |
% koncime |
| return; |
| end |
| B=inv(A); |
% inverze matice |
| I=eye(n); |
% jednotkova matice |
| test=max(max(abs(A*B-I))); |
% vlastni test |
|
|
| |
| |
Příklad 12A.3: |
|
| |
Vytvořte funkci pro testování pseudoinverzní
matice.
Instrukce: Pseudoinverzní matici vynásobíme zleva a
zprava původní maticí a porovnáme s původní maticí. Vlivem
zaokrouhlovacích chyb není jejich rozdíl vždy nulová matice.
Při jaké velikosti Hilbertovy matice bude chyba větší než jedna? |
|
| |
Řešení: |
|
| |
| function [test]=PseudoInverze(A) |
| % Numericky test pseudoinverze P matice A |
| % [test]=PseudoInverze(A); |
| % test ... maximalni prvek matice A*P*A - A |
| % A ... matice (m,n) |
| % P ... pseudoinverzni matice (n,m) podle Moora a Penroseho |
| P=pinv(A); |
% pseudoinverze matice |
| test=max(max(abs(A*P*A-A))); |
% vlastni test |
|
|
| |
| |
Příklad 12A.4: |
|
| |
Vytvořte funkci pro řešení maticové
rovnice AX=B.
Instrukce: Použijeme pseudoinverzi matice soustavy.
Krokujte pro různé soustavy rovnic. |
|
| |
Řešení: |
|
| |
| function [X]=Vyres(A,B) |
| % Bezpecne reseni soustavy linearnich rovnic AX=B |
| % [X]=Vyres(A,B); |
| % X ... reseni soustavy (n,s) |
| % A ... matice soustavy (m,n) |
| % B ... matice pravych stran (m,s) |
| ma=size(A,1); |
% urceni poctu radku matice A |
| n=size(A,2); |
% urceni poctu sloupcu matice A |
| mb=size(B,1); |
% urceni poctu radku matice B |
| s=size(B,2); |
% urceni poctu sloupcu matice B |
| if ma~=mb |
% nevhodne rozmery matic |
| X=zeros(n,s); |
% rezignace |
| else |
| X=pinv(A)*B; |
% nadherne, ze |
| end |
|
|
| |
| |
Příklad 12A.5: |
|
| |
Vytvořte funkci pro proložení funkce
y = a + b*x + c*sin(x) + d/(|x|+1) metodou nejmenších čtverců.
Instrukce: Báze prostoru funkcí je tvořena funkcemi:
1, x, sin x, 1/x. Matici soustavy vytvoříme ze čtyř sloupcových
vektorů hodnot funkce báze. Pak použijeme pseudoinverzi.
Krokujte pro přeurčené i nedourčené problémy. |
|
| |
Řešení: |
|
| |
| function [par,LSQ]=ProlozLSQ(x,y) |
| % Prolozeni funkce y=a + b*x + c*sin(x) + d/(|x|+1) metodou nejmenších ctvercu. |
| % [par,LSQ]=ProlozLSQ(x,y); |
| % par ... vektor nalezenych parametru [a b c d] |
| % LSQ ... soucet ctvercu odchylek |
| % x ... vektor nezavisle promenne (1,m) |
| % y ... vektor zavisle promenne (1,m) |
| % test rozmeru |
| if size(x,1)~=1 | size(y,1)~=1 | size(x,2)~=size(y,2) |
| par=[]; |
% reseni neexistuje |
| return; |
| end |
| x=x'; y=y'; |
| % zacatek konkretni casti |
| phi1=ones(size(x)); |
% prvni sloupec |
| phi2=x; |
% druhy sloupec |
| phi3=sin(x); |
% treti sloupec s pouzitim vektorove funkce |
| phi4=phi1./(abs(x)+1); |
% ctvrty sloupec a teckova konvence |
| % konec konkretni casti |
| A=[phi1 phi2 phi3 phi4]; |
% sestaveni matice soustavy |
| par=pinv(A)*y; |
% bez pseudoinverze ani ranu |
| LSQ=norm(A*par-y).^2; |
| % staci urcit druhou mocninu normy chyboveho vektoru |
|
|
| |
| |
Příklady pro samostatné vypracování |
|
| |
Příklad 12B.1: |
|
| |
Modifikovaná Hilbertova matice má
stejné absolutní hodnoty prvků jako původní matice, ale znaménka
jejích prvků se šachovnicově střídají tak, že lichému součtu indexů
odpovídá záporné znaménko. Sestavte funkci ModiHilbert pro její
výpočet a použijte jí k testování inverze a pseudoinverze |
|
| |
| |
Příklad 12B.2: |
|
| |
Marquardtova metoda řešení soustavy lineárních
rovnic je dána vzorcem x=inv(A’*A+lambda*I)*A’*b pro lambda>0.
Součet čtverců odchylek je roven normě vektoru A*x-b.
Sestavte funkci, která vrací jak řešení, tak součet čtverců
odchylek. Pro jakou kladnou hodnotu parametru lambda je součet
čtverců odchylek minimální? |
|
| |
| |
Příklad 12B.3: |
|
| |
Bázi můžeme obecně chápat jako vektorovou
funkci vektorové proměnné. Vytvořte bázi z funkcí:
1, x1, x2, 1/x1, 1/x2, x1*x2, x1/x2, x2/x1, x1/(x1+x2)
jako Matlabovskou funkci se vstupním vekrorem délky 2 a
výstupním vektorem délky 9. Pak vytvořte obecnou funkci pro
proložení dat metodou nejmenších čtverců. Vstupní parametry
funkce by měly být: název báze, matice x a vektor y. Výstupními
parametry by měly být vektor konstant regrese a součet čtverců
odchylek. |
|
| |
| |
Příklad 12B.4: |
|
| |
Vytvořte bázi pro polynom n-tého
stupně jedné proměnné. |
|
| |
| |
Příklad 12B.5: |
|
| |
Vytvořte bázi pro lineární
funkci n-proměnných. |
|
| |
| |
Příklad 12B.6: |
|
| |
Vytvořte bázi pro kvadratickou
funkci n-proměnných. |
|
| |
| |
Příklad 12B.7: |
|
| |
Teplotní pole je dáno čtvercovou sítí o
rozměrech m, n. Na horním okraji pole je teplota 100 stupňů
Celsia. Ostatní body na okraji pole mají teplotu 20 stupňů Celsia.
Ve zbylých (m-2)*(n-2) bodech sítě je teplota rovná průměrné
teplotě v sousedních bodech. Vygenerujte příslušnou soustavu
lineárních rovnic a vyřešte jí. Sousedství v síti můžeme
chápat různým způsobem: 4 sousedi do kříže, 4 sousedi diagonálně,
8 sousedů se stejnými vahami nebo 8 sousedů s vahami 2 a 1.
Porovnejte všechny čtyři definice sousedství. Znázorněte
graficky. |
|
| |
| |
Seznam použitých příkazů |
|
| |
rank, det, inv, pinv a hilb
|
|
|