4. laborgyakorlat
Numerikus számítások fizikusoknak
Contents
clear all, close all
Gauss-módszer főelemkiválasztással
A lebegőpontos számításokból eredő hibákat csökkenthetjük részleges- vagy teljes főelemkiválasztással. A részleges főelemkiválasztásnál sorcserét hajtunk végre úgy, hogy a főelem legyen a legnagyobb abszolút értékű elem az adott oszlop hátra lévő részében (A(k:n,k)-ban). Teljes főelemkiválasztásnál sor és oszlopot is cserélünk és az A(k:n,k:n+1) almátrixban lévő abszolút értékben legynagyobb elemet cseréljük a főelem helyére.
Általános LU-felbontás
A részleges főelemkiválasztással elvégzett Gauss-módszerrel egy
alakú ún. általános LU-felbontáshoz jutunk, ahol P megfelelő permutációs mátrix. L elemei most 1-nél nem nagyobb abszolút értékűek.
H=hilb(5) [L,U,P]=lu(H) % Hilbert-mátrix általános LU-felbontása P\(L*U) %=H
H = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 L = 1.0000 0 0 0 0 0.3333 1.0000 0 0 0 0.5000 1.0000 1.0000 0 0 0.2000 0.8000 -0.9143 1.0000 0 0.2500 0.9000 -0.6000 0.5000 1.0000 U = 1.0000 0.5000 0.3333 0.2500 0.2000 0 0.0833 0.0889 0.0833 0.0762 0 0 -0.0056 -0.0083 -0.0095 0 0 0 0.0007 0.0015 0 0 0 0 -0.0000 P = 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 1 0 ans = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111
Cholesky-felbontás
Cholesky-felbontás szimmetrikus, pozitív definit mátrixokra. Műveletszáma kb. fele a Gauss-módszerének
G=chol(H)
G'*G %=H
G = 1.0000 0.5000 0.3333 0.2500 0.2000 0 0.2887 0.2887 0.2598 0.2309 0 0 0.0745 0.1118 0.1278 0 0 0 0.0189 0.0378 0 0 0 0 0.0048 ans = 1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111
Teli- és ritka mátrixok
Egy mátrix ritka, ha a nemnulla elemek száma O(n) nagyságrendű. A sparse paranccsal lehet egy mátrixot ritkává (ekkor a Matlab tudja, hogy ritka mátrixról van szó) alakítani. A full parancs ritkából csinál teli mátrixot. spy ábrázolja a mátrix nemnulla elemeinek szerkezetét.
A=eye(20); whos A A=sparse(A); whos A spy(A)
Name Size Bytes Class Attributes A 20x20 3200 double Name Size Bytes Class Attributes A 20x20 324 double sparse
![](nsz2017_4_01.png)
A=toeplitz([2,1,zeros(1,8)]) spy(A)
A = 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2 1 0 0 0 0 0 0 0 0 1 2
![](nsz2017_4_02.png)
Futási idők összehasonlítása
A=rand(1000); % véletlen mátrix [L,U,P]=lu(A); B=A'*A; % szimmetrikus, pozitív definit mátrix C=toeplitz([2,-1,zeros(1,998)]); % tridiagonális mátrix D=sparse(C); b=rand(1000,1); tic, A\b; altalanos_ido=toc tic, U\(L\(P*b)); altalanos_LUval_ido=toc tic, B\b; szimm_poz_def_ido=toc tic, C\b; teli_tridiag_matrix_ido=toc tic, D\b; ritka_tridiag_matrix_ido=toc
altalanos_ido = 0.5663 altalanos_LUval_ido = 0.0116 szimm_poz_def_ido = 0.3485 teli_tridiag_matrix_ido = 0.3140 ritka_tridiag_matrix_ido = 2.4637e-004
Ciklusok
Azokat a programrészeket, melyeket a programmal ismételtetni szeretnénk ciklusoknak hívjuk. Kétféle ciklust használunk: for ill. while ciklust.
Példa for ciklusra
osszeg=0; for i=1:2:999 osszeg=osszeg+i^2; end osszeg
osszeg = 166666500
Sokszor érdemes kiváltani vektor- ill. mátrixműveletekkel
sum([1:2:999].^2) % gyorsabb megvalósítás
ans = 166666500
Példa while ciklusra
osszeg=0; i=1; while i<=999 osszeg=osszeg+i^2; i=i+2; end osszeg
osszeg = 166666500
Egyszerű logikai műveletek
1 az igaz, 0 a hamis logikai értéket jelenti
x=3; y=5; x==y % egyenlő x<y % kisebb x>y % nagyobb x<=y % kisebb vagy egyenlő x>=y % nagyobb vagy egyenlő x~=y % nem egyenlő x>2 & x<4 % és x>2 | y>6 % vagy ~(y==6) % tagadás
ans = 0 ans = 1 ans = 0 ans = 1 ans = 0 ans = 1 ans = 1 ans = 1 ans = 1
Iterációs egyenletrendszer-megoldás
Egy vektoriteráció határértékeként állítjuk elő a megoldást
Jacobi-iteráció
A=[7,2,3;-1,6,4;1,4,8]; b=[2;5;7]; x_pontos=A\b % Ez a Matlab által adott megoldás d=diag(A); % Az A diagonális elemeinek vektora. B=A-diag(d); % Ugyanaz, mint A, csak a főátló helyén nullák állnak. x=[1;2;3]; % Egy tetszőleges kezdővektor, ahonnét az iterációt indítjuk. for i=1:25 x=(b-B*x)./d end
x_pontos = -0.1193 0.3303 0.7248 x = -1.5714 -1.0000 -0.2500 x = 0.6786 0.7381 1.5714 x = -0.5986 -0.1012 0.4211 x = 0.1341 0.4528 1.0004 x = -0.2724 0.1887 0.6318 x = -0.0390 0.3667 0.8147 x = -0.1682 0.2837 0.6965 x = -0.0939 0.3410 0.7542 x = -0.1349 0.3149 0.7163 x = -0.1112 0.3333 0.7344 x = -0.1243 0.3252 0.7222 x = -0.1167 0.3311 0.7279 x = -0.1209 0.3286 0.7240 x = -0.1185 0.3305 0.7258 x = -0.1198 0.3297 0.7246 x = -0.1190 0.3303 0.7251 x = -0.1194 0.3301 0.7247 x = -0.1192 0.3303 0.7249 x = -0.1193 0.3302 0.7248 x = -0.1192 0.3303 0.7248 x = -0.1193 0.3303 0.7248 x = -0.1193 0.3303 0.7248 x = -0.1193 0.3303 0.7248 x = -0.1193 0.3303 0.7248 x = -0.1193 0.3303 0.7248