műszaki dokumentáció
A módosított mátrix elemeinek (figyelembe véve a fenti, és annak szimmetria) explicit formulák kapunk:
a'rp = ponty - sarq (R = p, r = q !!)
a'rq = carq + Sarp (R! = p, r! = q)
a'pp = c 2 app + s 2 AQQ - 2csapq
a'qq = s 2 2 app + c AQQ + 2csapq
a'pq = (c 2 - s 2) APQ + cs (app - AQQ) Az ötlet a módszer Jacobi -, hogy megpróbálja nulla off-diagonális elemei egy sor forgatások. Annak érdekében, hogy visszaállítsa a'pq. forgási szög szerinti fenti képletek a következők: q = gyermekágy (2 f) = (c 2 - s 2) / (2cs) = (AQQ - app) / (2apq) .Polagaya t = s / c, meghatározására t megérkezik a következő egyenlet: t 2 + 2t q - 1 = 0.Menshy modulo gyökere ennek az egyenletnek felel meg a rotációs szögben kisebb p / 4; a választás, ez a szög minden szakaszában biztosítja a legstabilabb algoritmus. A képlet másodfokú egyenlet megoldása a nevezőben diszkrimináns kisebb gyökérzettel definiáljuk :. T = jel (q) / (| q | + sqrt (2 q + 1)) Ha az érték q olyan, hogy q 2 kifut, feltételezzük, t = 1 / (2 q). A t értéke határozza c és s: c = 1 / sqrt (t 2 +1), s = tc.Pri numerikus mátrix elemeinek módosítások kívánjuk csökkenteni kerekítési hibákat. Az ötlet az, hogy átírják őket, mint a korábbi érték, plusz egy kis korrekció adalékanyag. Ebben az esetben az együtthatók a mátrix transzformáció után fog kinézni: a'pp = app - tapq
a'rp = ARP - s (ARQ + GARP)
a'rq = ARQ + s (ARP - garq), ahol t (a tangensét fele forgási szög) úgy definiáljuk, mint t = s / (1 + c).
Konvergencia a Jacobi módszer becsülni lehet, figyelembe véve a négyzetösszeg az off-diagonális elemei S; transzformációs egyenletek biztosítani ezt az összeget a transzformáció után az alábbi összefüggés:
S = S r = s (| ars | 2!); S „= S - 2 | APQ | 2 .Poskolku ortogonális transzformáció, az összeget a diagonális elemeit, miközben növeli rendre | APQ | 2. összege a négyzetek az off-átlós elemek monoton csökken. Ahogy alulról korlátos nullával, és mint tudjuk, minden egyes alkalommal, hogy melyik akar APQ elem nulla, az összeg megy nullára.
Transzformáció után lánc fordul diagonális mátrix D. akár gépi nulla. A diagonális elemei sajátértékei az eredeti formájában a mátrix végre
D = V T AV. ahol V = P1P2P3. és Pi - Jacobi rotációs mátrix. A mátrix oszlopai V sajátvektorai. Ezek kiszámítási rekurzívan, minden egyes szakaszában a képletek: V „= VPI, ahol az eredeti V - azonosító mátrix. A konverziós elemeket a V általános képletű használt: v'rs = VRS (s = p, s = q !!)
v'rp = cvrp - svrq
v'rq = svrp + cvrq csökkentése kerekítési hibák fenti képletekben kell újraírni szempontjából g (mint a fenti a mátrix elemeinek A).
Az egyetlen megmaradt kérdés most az a stratégia, amely szeretné keresést elpusztult az elemeket. Az eredeti Jacobi algoritmus által 1846 minden szakaszában keresi a felső háromszög alakú a legnagyobb elem abszolút értékben és visszaáll. Ez egy nagyon racionális stratégia kézi számlálás, de elfogadhatatlan, hogy a számítógép, mint ahogy a tranzakciók száma minden egyes szakaszában a Jacobi forgás érdekében N 2 N helyett
A mi szempontunkból, a legjobb stratégia a ciklikus Jacobi módszer. ahol elemek eltávolítása rendezetten. Például ha csak elhalad sorok: P12, P13. p1n; P23. P2n ;. Belátható, hogy a konvergencia sebességét általában négyzet, mint az eredeti és a ciklikus Jacobi módszerrel. Azt mondjuk, hogy egy ilyen készlet az egymást követő, n (n-1) / 2 Jacobi átalakulás áthaladását.
A program alatt elhelyezett alapul algoritmusok [2, 3]; tartalmaz még két finomság:
- Az első három passzálással végzünk viszonyított relatív elforgatással az index (PQ) csak akkor, ha a | APQ |> e valamilyen küszöbértéket e = S0 / (5n 2), ahol S0 - összege abszolút off-diagonális elemei a felső háromszög mátrix.
- Miután négy halad abban az esetben, ha | APQ |<<|app | и |apq |<<|aqq |, полагаем |apq |=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10 D+2 раз, где D -- число значащих десятичных цифр.
Tipikus mátrixok igényelnek 6-10 menetben konvergencia eléréséhez, vagy a 2-től, 3n 5N 2 forgatások. Minden forgatás a végzés kötelezi 4n műveletek az összeadás és a szorzás, úgy, hogy a teljes költség nagyságrendileg 12n hogy 20N 3 3 műveleteket. Kiszámításakor a sajátvektorok a sajátértékek túlmenően, a műveletek száma növekszik a forgási 4n 6n, ami azt jelenti, hogy csak a többlet 50 százalék.
A program Jacobi
#include
#include „nrutil.h” / * Ez azonosít néhány dologra, mint a memória kiosztás * /
/ * Konvertálása az elemek forgása * /
#define ROTÁLJA (a, i, j, k, l) g = a [i] [j]; h = [K] [l]; a [i] [j] = gs * (H + g * tau ); a [K] [l] = h + s * (gh * tau)
/ * Maximális menetszám * /
#define MAXSWEEP 50
/ * Jacobi program kiszámítja az összes sajátvektor és sajátérték
valós szimmetrikus mátrix egy [1. n] [1. n]. A kilépő elpusztult
Overdiagonal minden elemét a. d [1. n] visszatér sajátértékei
mátrix, v [1. n] [1. n] - a oszlopok normalizált sajátvektorai
nrot - a fordulatok száma Jacobi, amelyek az ehhez szükséges programot.
* /
void Jacobi (float ** a, int n, float * d, float ** v, int * nrot) <
int j, IQ, ip, i;
float Tresh, théta, tau, t, sm, s, H, G, C, * b, * z;
/ * Memóriát kiosztani ideiglenes használatra vektorok (nyilatkozat nrutil.h) * /
b = vektor (1, n); Z = vektor (1, n);
/ * Inicializálása v, mint egy egység mátrix * /
a (IP = 1; ip<=n;ip++) <
az (IQ = 1; iq<=n;iq++) v[ip][iq]=0.;
v [ip] [ip] = 1.;
>
/ * Inicializálja a átlója b a, z - nulla * /
a (IP = 1; ip<=n;ip++)
/ * Az első fordulatok száma nullára * /
* Nrot = 0;
/ * Nem * több MAXSWEEP halad /
for (i = 1; i<=MAXSWEEP;i++) <
/ * Összegét számolja ki az abszolút off-diagonális elemei * /
A (sm = 0., ip = 1; ip<=n;ip++) for(iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]);
/ * Diagonális mátrix -> normál kimenet * /
ha (sm == 0). <
free_vector (Z, 1, N); free_vector (b, 1, N); return;
>
/ * A küszöbérték forgatható * /
Tresh = (i<4. 0.2*sm/(n*n). 0.);
/ * Pass végzi sorok minden sorának oszlop * /
a (IP = 1; ip<=n-1;ip++) for(iq=ip+1;iq<=n;iq++) <
/ * I nyomon követni egy kis eset tag után 4 átjárók * /
g = 100. * FAB (a [ip] [IQ]);
ha (i> 4 (Float) FAB (d [ip] + g) == (float) FAB (d [ip])
(Float) FAB (d [IQ] + g) == (float) FAB (d [IQ])) egy [ip] [IQ] = 0.;
/ * A esetében egy kis elem az első 3 menetben
(Process csak átlépte a küszöbértéket) * /
else if (FAB (a [ip] [IQ])> Tresh) <
h = D [ip] -d [IQ];
/ * Számítsuk ki a t = s / c a képlet a négyzetgyöke egyenlet * /
if ((float) (FAB (H) + g) == (float) FAB (h)) t = a [ip] [iq] / h;
más <
theta szögeknél = 0,5 * h / a [ip] [IQ];
t = 1 ./ (FAB (théta) + sqrt (1. + théta * theta));
ha (théta<0.) t = -t;
>
/ * Compute c, s, a tau, et al. Módosítása átlós. Reset (ip, iq) egy elem * /
c = 1 / sqrt (1 + t * t); s = t * C; tau = s / (1 + c); h = t * a [ip] [IQ];
z [ip] - = H; z [IQ] + = H; d [ip] - = H; d [IQ] + = H;
egy [ip] [IQ] = 0.;
/ * Forgás 1<=j
/ * Kapcsolja ip
/ * Kapcsolja át iq
/ * Kiegészítés a sajátvektor mátrix * /
a (j = 1; j<=n;j++)
/ * Növekvő számláló forgatások * /
++(* Nrot);
>
>
/ * Kitesz egy átlós és újrainicializálása z * /
a (IP = 1; ip<=n;ip++) <
b [ip] + = Z [ip]; d [ip] = b [ip]; z [ip] = 0.;
>
>
/ * Ha már itt tartunk, a fordulatok száma meghaladta a határértéket. nhiba funkció (output
hibafeltárást) ismertetett nyilatkozat nrutil.h. * /
n ( „Túl sok iterációinak rutin Jacobi”);
>
Figyeljük meg, hogy a program a fent bemutatott feltételezi, hogy a számítógép értelmezi az alább látható gép nulla nulla. Ha nem, akkor a program módosítani kell.
A kimenő e program, a sajátértékek nem rendezett nagyságrendű. Ha azt szeretnénk rendezni, akkor a következő program átrendezni hétvégén Jacobi program adatai, vagy más kutatási programok sajátvektorok és értékek (pl downstream). rendezési módszer - a közvetlen beillesztése - nagyságrendű az N 2- (lényegesen több, mint n kapcsolódás (N) a hatékony módszerek, de összehasonlítva a sorrendben N 3 eljárással Jacobi ez a gyorsulás elhanyagolható lenne).
programot eigsrt
/ * Eigsrt programot, fogadására bemeneten a sajátértékek és D [1. n] és a megfelelő
sajátvektorok v [1. n] [1. n], rendezni a sajátértékek a csökkenő
rendelni, így a megfelelő változtatásokat a tömbben vektorok. rendezési eljárás
közvetlen csatlakozással
* /
void eigsrt (float * d, float ** v, int n) <
int k, j, i;
float p;
/ * Pass összes értékét az övé. szám, kivéve az utolsót * /
for (i = 1; i
/ * Akár hogy egy permutáció és ahol a * /
a (j = i + 1, j<=n;j++) if(d[j]>= P) p = d [k = j];
/ * Ha igen, hogy ez * /
ha (k! = i) <
d [K] = d [i]; d [i] = p;
a (j = 1; j<=n;j++) <
p = v [j] [i]; v [j] [i] = v [j] [k]; v [j] [K] = p;
>
>
>