Documentație tehnică

Pentru elementele matricei modificate (luând în considerare cele de mai sus și simetria sa) se obțin formule explicite:

a'rp = crap - sarq (r! = p, r! = q)
a'rq = carq + sarp (r! = p, r! = q)






a'pp = c 2 ap + s 2 aqq - 2csapq
a'qq = s 2 app + c 2 aqq + 2csapq
a'pq = (c 2 - s 2) apq + cs (app - aqq) Ideea metodei Jacobi este de a încerca să zeroze elementele off-diagonale cu o serie de rotații. Pentru a zero la a'pq. unghiul de rotație în conformitate cu formulele prezentate mai sus sunt următoarele: q = pat copii (2 f) = (c 2 - s 2) / (2CS) = (AQQ - app) / (2apq) .Polagaya t = s / c, pentru a determina t, ajungem la următoarea ecuație: t 2 + 2t q-1 = 0. Cea mai mică rădăcină a acestei ecuații corespunde rotației cu un unghi mai mic decât p / 4; Alegerea acestui unghi la fiecare etapă oferă cel mai stabil algoritm. Folosind formula de rezolvare a unei ecuații pătratice în discriminante numitor rădăcină mai mic este definit ca :. T = semnul (q) / (| q | + sqrt (2 q + 1)) Dacă valoarea q astfel încât q 2 va overflow, să presupunem t = 1 / (2 q). Din valoarea t sunt determinate c și s: c = 1 / sqrt (t 2 +1), s = tc.Pri elemente ale matricei numerice modificări ne propunem scăderea erorilor de rotunjire. Ideea este să le rescrieți ca o valoare anterioară plus un aditiv de corecție mic. În acest caz, coeficienții de matrice după conversie vor arăta astfel: a'pp = app - tapq
a'rp = arp - s (arq + garp)
a'rq = arq + s (arp - garq), unde t (tangenta jumătate a unghiului de rotație) este definită ca t = s / (1 + c).

Convergența metodei Jacobi poate fi estimată luând în considerare suma pătratelor elementelor off-diagonale ale lui S; ecuațiile de transformare dau pentru această sumă după transformare următoarea relație:

S = S r '= S (ars | 2); S '= S - 2 | apq | 2. Deoarece transformările sunt ortogonale, suma elementelor diagonale crește, respectiv, cu | apq | 2. Suma pătratelor elementelor off-diagonale scade monotonic. Deoarece este limitat de la zero cu zero și din moment ce putem alege întotdeauna ce fel de element apq vrem să zero, suma va tinde la zero.

După lanțul de transformări, matricea D este obținută în diagonală într-o mașină zero. Elementele diagonale formează valorile proprii ale matricei inițiale A. Deoarece

D = V T AV. unde V = P1P2P3. și Pi sunt matrici de rotație Jacobi. Coloanele matricei V sunt vectori proprii. Ei se calculează recurent, în fiecare etapă prin formulele: V '= VPi, unde inițial V este matricea de identitate. Pentru a converti elementele lui V, folosim formulele: v'rs = vrs (s! = P, s! = Q)
v'rp = cvrp - svrq
v'rq = svrp + cvrq Pentru a reduce eroarea de rotunjire, formulele de mai sus ar trebui să fie rescrise în termeni de g (ca mai sus pentru elementele matricei A).

Singura întrebare rămasă este strategia prin care trebuie să sortați obiectele distruse. Algoritmul original Jacobi din 1846, în fiecare etapă, a căutat cel mai mare element din regiunea triunghiulară superioară și la zero. Aceasta este o strategie foarte rațională pentru calculul manual, dar este inacceptabil pentru un calculator, deoarece face numărul de operații în fiecare etapă a rotației Jacobi a ordinului N 2 în loc de N.

Pentru scopurile noastre, cea mai bună strategie este metoda ciclică Jacobi. unde elementele sunt șterse într-o ordine strictă. De exemplu, puteți trece pur și simplu prin linii: P12, P13. p1n; P23. P2n ;. Se poate arăta că rata de convergență va fi patratică în cazul general, atât pentru metoda inițială, cât și pentru metoda ciclică Jacobi. Noi numim un astfel de set de transformări succesive n (n-1) / 2 Jacobi printr-un pasaj.

Programul de mai jos se bazează pe algoritmi de la [2, 3]; conține două subtilități:

• In primele trei pasaje vom efectua o rotație în raport cu indexul (PQ) numai atunci când | APQ |> e pentru un anumit prag e = S0 / (5n 2), unde S0 - suma absolute a elementelor off-diagonală a matricei superior triunghiulară.
• După patru treceri în cazul în care | apq |<<|app | и |apq |<<|aqq |, полагаем |apq |=0 и пропускаем ротацию. Критерием сравнения является превышение диагональных элементов в 10 D+2 раз, где D -- число значащих десятичных цифр.

Matricile tipice necesită 6 până la 10 pasaje pentru realizarea convergenței sau de la 3n 2 la 5n 2 rotații. Fiecare rotație necesită operațiuni de înmulțire și adăugare de 4n, deci costurile totale sunt de ordinul a 12n 3 până la 20n 3 operațiuni. Atunci când se calculează vectorii proprii împreună cu valorile proprii, numărul de operații pe rotație crește de la 4n la 6n, ceea ce înseamnă un exces de numai 50%.

Programul Jacobi

#include
#include "nrutil.h" / * Unele utilitare precum alocarea memoriei sunt definite aici * /

/ * Conversia elementelor cu rotație * /
#define ROTATE (a, i, j, k, l) g = a [i] [j]; h = a [k] [l]; o [i] [j] = gs * (h + g * tau ), o [k] [l] = h + s * (gh * tau)
/ * numărul maxim de treceri * /
#define MAXSWEEP 50

/ * Programul jacobi calculează toate valorile proprii și vectorii proprii
a matricei reale simetrice a [1. n] [1. n]. La ieșire sunt distruse
toate elementele off-diagonale a. d [1. n] returnează valorile proprii
matrice, v [1. n] [1. n] - în coloane vectorii proprii normalizați,
nrot este numărul de rotații Jacobi necesare pentru acest program.
* /
void jacobi (float ** a, int n, float * d, float ** v, int * nrot) <
int j, iq, ip, i;
float tresh, theta, tau, t, sm, s, h, g, c, * b, * z;
/ * alocarea memoriei pentru vectorii folosiți temporar (declarații în nrutil.h) * /
b = vector (1, n); z = vector (1, n);
/ * initialize v ca o unitate matrice * /
pentru (ip = 1, ip<=n;ip++) <
pentru (iq = 1; iq<=n;iq++) v[ip][iq]=0.;
v [ip] [ip] = 1;
>
/ * initializeaza b cu diagonala a, z - zero * /
pentru (ip = 1, ip<=n;ip++)
/ * la început numărul de rotații este zero * /
* nrot = 0;
/ * nu mai trece MAXSWEEP passes * /
pentru (i = 1; i<=MAXSWEEP;i++) <
/ * calcula suma modulelor elementelor off-diagonale * /
pentru (sm = 0, ip = 1; ip<=n;ip++) for(iq=ip+1;iq<=n;iq++) sm += fabs(a[ip][iq]);
/ * matrice diagonală -> ieșire normală * /
dacă (sm == 0.) <
free_vector (z, 1, n); vectorul liber (vector b, 1, n); return;
>
/ * valoarea pragului elementului pentru rotire * /
tresh = (i<4. 0.2*sm/(n*n). 0.);
/ * trece prin linii, în fiecare rând pe coloane * /
pentru (ip = 1, ip<=n-1;ip++) for(iq=ip+1;iq<=n;iq++) <
/ * urmări cazul unui element mic după 4 treceri * /
g = 100. * fabs (a [ip] [iq]);
dacă (i> 4 (float) fabs (d [ip] + g) == (float) fabs (d [ip])
(float) fabs (d [iq] + g) == (float) fabs (d [iq])) a [ip] [iq] = 0;
/ * și cazul unui element mic pe primele 3 treceri
(procesați numai cele care au depășit pragul) * /
altfel dacă (fabs (a [ip] [iq])> tresă) <
h = d [ip] -d [iq];
/ * calculează valoarea lui t = s / c prin formula rădăcinii ecuației patrate * /
dacă (float) (fabs (h) + g) == (float) fabs (h)) t = a [ip] [iq] / h;
altfel <
theta = 0,5 * h / a [ip] [iq];
t = 1 ./ (fabs (theta) + sqrt (1. + theta * theta));
dacă (theta<0.) t = -t;
>
/ * calculează c, s, tau, etc. Schimbați diagonala. Zero (ip, iq) - element * /
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;
a [ip] [iq] = 0;
/ * rândul său, la 1<=j pentru (j = 1;<=ip-1;j++)
/ * rândul său, la ip pentru (j = ip + 1;<=iq-1;j++)
/ * rotație cu iq pentru (j = iq + 1;<=n;j++)
/ * aditiv pentru matricea vectorilor proprii * /
pentru (j = 1;<=n;j++)
/ * creșterea numărătorului de rotație * /
++(* nrot);
>
>
/ * adăugați până la diagonală și reinitializați z * /
pentru (ip = 1, ip<=n;ip++) <
b [ip] + = z [ip]; d [ip] = b [ip]; z [ip] = 0;
>
>
/ * dacă suntem aici, atunci numărul de rotații a depășit limita. Funcția de gros
diagnoza erorilor) este descrisă în declarația din nrutil.h. * /
nerror ("Prea multe iterații în rutina jacobi");
>

Rețineți că programul prezentat mai sus presupune că calculatorul tratează un număr sub zero mașinii ca zero. Dacă nu, programul trebuie modificat.

La ieșirea din acest program, valorile proprii nu sunt în ordine de mărime. Dacă este necesară sortarea, următorul program va modifica ieșirea programului jacobi sau a altor programe pentru căutarea vectorilor proprii și valori (de exemplu, localizate mai jos). Metoda de sortare - introducerea directă - este de ordinul N 2 (mult mai mare decât Nlog (N) pentru metode eficiente, dar această accelerare ar fi nesemnificativă comparativ cu ordinea N 3 pentru procedura jacobi).

Programul eigsrt

/ * Program eigsrt, obținând la valorile proprii d [1. n] și corespunzătoare
eigenvectors v [1. n] [1. n] sortează valorile proprii în ordinea descrescătoare
ordine, făcând transformările adecvate în matricea de vectori. Metoda de sortare
inserarea directă
* /
void eigsrt (float * d, float ** v, int n) <
int k, j, i;
float p;
/ * trece prin toate valorile proprii. numere, cu excepția ultimului * /
pentru (i = 1; i p = d [k = i];
/ * dacă să facă o permutare și unde * /
pentru (j = i + 1;<=n;j++) if(d[j]>= p) p = d [k = j];
/ * dacă este cazul, produceți-o * /
dacă (k! = i) <
d [k] = d [i]; d [i] = p;
pentru (j = 1;<=n;j++) <
p = v [j] [i]; v [j] [i] = v [j] [k]; v [j] [k] = p;
>
>
>

Articole similare

• Calculul numărului de posturi - reconstrucția stației de întreținere a autovehiculelor - altele

• Maeștri de delphi, salvarea și încărcarea datelor în obiecte pe exemplul colecțiilor (documentație,

• Genius atunci când aparatul foto devine o jucărie, o documentație de la A la Z

Trimiteți-le prietenilor: