Gauss

Size: px

Start display at page:

Download "Gauss"

きのこちゃわんや
5 years ago
Views:

1 15 1 LU LDL T 6 : 1g00p013-5

2 Gauss Gauss Gauss LU LU LU Cholesky Cholesky

3 4.3 Cholesky LDL T LDL T LDL T

4 A 49 3

6 7.1 LU LDL T

7 1 6

8 1.1 1 LDL T A A 1 1 LU LDL T 7

9 1.2 1 LU LDL T 1.3 2, Gauss 3,Gauss LU 4,LDL T Cholesky. 5,LDL T 6, 7, 8,. 8

10 2 Gauss 9

11 2.1 Gauss LU 2.2 Gauss x 1, x 2, x 3 a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 (2.1) a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 (2.2) a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 (2.3) 2x 1 + 5x 2 + 7x 3 = 23 (2.4) 4x x x 3 = 58 (2.5) 8x x x 3 = 132 (2.6) Gauss (2.4) (2.5) (2.6) 2x 1 + 5x 2 + 7x 3 = 23 (2.7) 13x x 3 = 12 (2.8) 29x x 3 = 40 (2.9) (2.8) (2.9). 2x 1 + 5x 2 + 7x 3 = 23 (2.10) 3x 2 + 6x 3 = 12 (2.11) 10 4x 3 = 4 (2.12)

12 Gauss (2.12) x 3 = 1 (2.11) x 2 = (12 6x 3 )/3 = (12 6 1)/3 = 2 (2.10) x 1 = (23 5x 2 7x 3 )/2 = ( )/2 = 3 a ii Gauss a 11 0 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 a 21 0 ( a 31 ) 11

13 0 a 11 0 x 1 a 12,,a 1n a 12,,a 1n a 11 a 12,,a 1n a Gauss a 11 x 1 + a 12 x 2 + a 13 x 3 = b 1 a 21 x 1 + a 22 x 2 + a 23 x 3 = b 2 a 31 x 1 + a 32 x 2 + a 33 x 3 = b 3 k = 1, 2,..., n 1 k i (i = k + 1,..., n) (A b) = (A (1) b (1) ) (A (1) b (1) ) (A (2) b (2) ) (A (k) b (k) ) (A (k+1) b (k+1) ) (A (n) b (n) ) A (n) = u 11 u 12 u 1n u 22 u 2n u nn, b(n) = y A (k) = (a (k) ij ), b (k) = (b (k) i ) Gauss 12

14 k = 1, 2,..., n 1 Gauss i = k + 1, k + 2,..., n m ik = a (k) ik /a(k) kk j = k + 1, k + 2,..., n b (k+1) i a (k+1) ij = b (k) i = a (k) ij m ik b (k) k m ik a (k) kj 1 a (1) 11 x 1 + a (1) 12 x a (1) 1,n 1x n 1 + a (1) 1n x n = b (1) 1 a (2) 22 x a (2) 2,n 1x n 1 + a (2) 2n x n = b (2) 2. a (n 1) n 1,n 1x n 1 + a (n 1) n 1,nx n a (n) nn x n = b (n 1) n 1 = b (n) n x n = b(n) n a (n) nn (2.13) 2 x n 1 x k = b(k) k a (k) k,k+1 x k+1 a (k) k,k+2 x k+2 a (k) kn x n a (k) kk (2.14) x 1 C 13

15 for (i=1;i<=n;i++){ pivot = a[i][i]; Gauss C for (j=i+1;j<=n;j++){ p=a[j][i]/pivot; a[j][i]=0.0; for (k=i+1;k<=n;k++) a[j][k] -= p*a[i][k]; b[j] -= p*b[i]; x[n] = b[n]/a[n][n]; for (i=n-1;i>=1;i--){ x[i]=b[i]; for (j=i+1;j<=n;j++) x[i] -= a[i][j]*x[j]; x[i] /= a[i][i]; / (n i) + (n i)(n i + 1) + (n i) + 1 = (n i)(n i + 3) + 1 / (n i)(n i + 1) (n i 1) = (n i)(n i + 2) / n [(n i)(n i + 3) + 1] = 1 + (n 2 2ni + i 2 + 3n 3i + 1) i=1 n 1 i=1 14

16 n 1 = 1 + (n 2 n 1 + 3n + 1) 1 (2n + 3) i=1 = 1 + (n 2 + 3n 1)(n 1) (2n + 3) + (n 1)n(2n 1) 6 / n 1 i=1 (n i)(n i + 2) = n 1 i=1 = (n 2 + 2n) (n 2 2ni + i 2 + 2n 2i) n 1 i=1 n 1 1 (2n + 2) = (n 2 + 2n)(n 1) (2n + 2) = 2n3 + 3n 2 5n 6 i=1 i=1 = 2n3 + 6n 2 3n 6 i + n 1 i 2 i=1 (n 1)n 2 + i + n 1 i 2 i=1 (n 1)n 2 (n 1)n(2n 1) 6 n 3 /3 n 3 /3 1 Gauss O(n 3 ) 2.5 Gauss LU 15

17 3 LU 16

18 3.1 Gauss LU LU 3.2 LU Ax = b x = A 1 b Gauss b Gauss M 1 = 1 m m m n A (1) = A M 1 A (1) = A (2) = a (1) 11 a (1) 12 a (1) 1n 0 a (2) 22 a (2) 2n a (2) n2 a (2) nn 17

19 M 1 A (1) k M K =, m jk = a (k) jk.. m k+1,k 1 /a(k) kk (3.1).. m k+2,k m n,k 0 1 k = 1, 2,..., n 1 Gauss i = k + 1, k + 2,..., n m ik = a (k) ik /a(k) kk j = k + 1, k + 2,..., n b (k+1) i a (k+1) ij = b (k) i = a (k) ij m ik b (k) k m ik a (k) kj M k A (k) = A (k+1) (3.2) M n 1 M n 2...M 2 M 1 A = A (n) 18

20 = U = a (1) 11 a (1) 12 a (1) 13 a (1) 1n a (2) 22 a (2) 23 a (2) 2n 0 a (1) 33 a (3) 3n.... a (n) nn (3.3) U 0 M n 1 M n 2...M 1 b = b (n) = b (1) 1 b (2) 2. b (n) n (3.4) Ax = b 3.3 U Ux = b (n) (3.5) M k Mk =.. m k+1,k 1.. m k+2,k m n,k 0 1, m jk = a (k) jk /a(k) kk (3.6) M k I L Mk 1 L = M 1 1 M M 1 n 1 (3.7) 19

21 1 m m 31 m 32 1 L = m n1 m n2 m n3 m n,n 1 1, m jk = a (k) jk /a(k) kk (3.8) L Gauss A A = LU (3.9) L U A LU A LU L,U Ax = b Ax = b LUx = b y Ly = b (3.10) Ux = y (3.11) x Ly = b m kk = 1 y 1 = b 1 k = 2, 3,..., n y k = b k k 1 j=1 m kj y j Ux = y 20

22 x n = y n /a (n) nn k = n 1, n 2,..., 1 n x k = y k a (k) kj x j j=k+1 /a (k) kk 3.3 LU 3 3 A A = Gauss A = b A L U 21

23 A L U = (3.12) A L U LU LU for( i = 1; i <= n; i++ ){ for( j = 1; j <= n; j++ ){ LU C if( i == j ){ l[i][j] = 1; else{ l[i][j] = 0; u[i][j] = a[i][j]; for( i = 1; i <= n; i++ ){ for( j = i+1; j <= n; j++ ){ l[j][i] = u[j][i] / u[i][i]; for( k = i; k <= n; k++ ){ u[j][k] = u[j][k] - u[i][k] * l[j][i]; Gauss a[1][1] 0 a[1][1] a[2][1] a kk a[k][k] 0 0 a kk 22

24 3.4 Gauss LU Cholesky 23

25 4 Cholesky 24

26 4.1 LU A Cholesky 4.2 Cholesky A A = S T S 1 A LDL T Cholesky A > 0 a (k) kk D 1 2 = a (1) 11 0 a (2) a (n) nn (4.1) S = D 1 2 L T A A = S T S (4.2) A S s 11 s 12 s 1n s 22 s 2n S = s nn S T S ij a ij (4.3) s 1i s ij + s 2i s 2j + + s ii s ij = a ij, (i j) (4.4) 25

27 s ii = a ii s ij = (a ij i 1 k=1 i 1 k=1 s 2 ki, (s 11 = a 11 ) s ki s kj )/s ii, (s 1j = a 1j /s 11 ) (4.5) s 11 = a 11 s 12, s 13,..., s 1n, s 22, s 23,... S S S T y = b, Sx = y Ax = b A 4.3 Cholesky Cholesky s ii = a ii s ij = (a ij i 1 k=1 i 1 k=1 s 2 ki, (s 11 = a 11 ) s ki s kj )/s ii, (s 1j = a 1j /s 11 ) (4.6) C 26

28 for (i = 1; i <= n; i++) { s = a[i][i]; for (k = 1; k < i; k++) s -= sqr(l[i][k]); if (s < 0) { Cholesky C exit(0); L[i][i] = sqrt(s); for (j =i + 1; j <= n; j++) { s = a[j][i]; for (k = 1; k < i; k++) s -= L[j][k] * L[i][k]; L[j][i] = s / L[i][i]; 4.4 A Cholesky LDL T 27

29 5 LDL T 28

30 5.1 A Cholesky LDL T 5.2 LDL T A LU U U U = DV (5.1) D = a (1) 11 0 a (2) 22 0 a (3) a (n) nn (5.2) V = 1 v 12 v 13 v 1n 1 v 23 v 2n 1 v 3n , v kj = a (k) kj /a(k) kk, k < j (5.3) A A t = A (n k) (n k) Gauss 1 29

31 A (k) s = a (k) kk a (k) k+1,k. a (k) nk a (k) k,k+1 a (k) kn a (k) k+1,k+1 a (k) k+1,n..... a (k) n,k+1 a (k) nn (5.4) A (k) s = A (k)t s (5.5) a (k) ij = a (k) ji, k i, j n (5.6) k 1 a ij = a (k 1) ij a(k 1) i,k 1 a (k 1) k 1,k 1 a (k 1) = a (k 1) a (k) = a (k) αβ βα ij ji a (k 1) k 1,j (5.7) A = A T a (1) αβ = a(1) βα A v kj = a(k) kj a (k) kk = a(k) jk a (k) kk = m jk (5.8) V L V = L T (5.9) A A A = LDL T (5.10) A LDL T 30

32 5.3 A LDL T D ii d i = a (i) ii L ij l ij A = LDL T LDL T A k=1,2,...,n i l kj d j l ij = a ki i = 1, 2,..., k 1 j=1 k l ki d i l ki = a kk i=1 l ii = 1 d 1 = a 11 k = 2, 3,..., n LDL T i = 1, 2,..., k 1 i 1 l ki = a ki l kj l ij d j /d i d k = a kk j=1 k 1 lkid 2 i i=1 0 d 1 = a 11 l 21, d 2, l 31, l 32, d 3,... D L A LDL T Ax = LDL T x = Ly, y = DL T x Ly = b, DL T x = y Ax = b 31

33 5.4 LDL T LDL T 1 1 2x 1 + x 2 + x 3 = 8 (5.11) x 1 + 3x 2 + 2x 3 = 11 (5.12) x 1 + 2x 2 + 4x 3 = 16 (5.13) x = x x (5.14) Ax = b (5.15) A LDL T A = LDL T (5.16) 1 d 11 1 l 21 l 31 = l 21 1 d 22 1 l 32 (5.17) l 31 l 32 1 d 33 1 d 11 d 11 l 21 d 11 l 31 = d 11 d 11 l d 22 d 11 l 21 l 31 + d 22 l 32 (5.18) d 11 d 11 l 21 l 31 + d 22 l 32 d 11 l d 22 l d = (5.19)

34 LU AX = b LDL T x = b LDL T x = b (5.20) DL T x = y (5.21) Ly = b (5.22) = y 1 y 2 y (5.23) y 1 y 2 y 3 = (5.24) DL T x = y x 1 x 2 x 3 = (5.25) x 1 x 2 x 3 = (5.26) 33

35 x 1 x 2 x 3 = (5.27) L D d 11 = a 11, l j1 = a j1 /d 1 (j = 2,..., n) d 22 = a 22 d 11 l 2 21, l j2 = (a 2j d 11 l 21 l j1 /d 22 )(j = 3,..., n) d 33 = a 33 d 11 l 2 31 d 22 l 2 32, l j3 = (a 3j d 11 l 31 l j1 d 2 l 32 l j2 )/d 33 (j = 4,..., n) d[1]=a[1][1]; for(k=2;k<=n;k++){ s=0;t=0; for(i=1;i<=k-1;i++){ LDL T C for(j=1;j<=i-1;j++){ s+=l[k][j]*l[i][j]*d[j]; l[k][i]=(a[k][i] - s)/d[i]; for(i=1;i<=k-1;i++){ t+=l[k][i]*l[k][i]*d[i]; d[k]=a[k][k]-t ; 34

36 5.5 A LDL T LU LDL T 35

37 6 36

38 6.1 C LU LDL T. 6.2,,. CPU memory OS Celeron 2.0GHz 512MB Windows XP cygwin+gcc 6.1: 6.3 LU 1. LDL T A 2. A a 11 2,3,...,n 37

39 A = n n n n n n n n n A 3. A x x i = i, (i = 1,..., n) 4. Ax = b b 5. A LU 6. Ly = b y 7. Ux = y x LDL T 1. A, x, b LU 2. A LDL T 3. Ly = b y 4. DL T x = y x

40 6.4 C LU LDL T 39

41 7 40

42 7.1 C LU LDL T "LU.txt" "LDLT.txt" time dimension 7.1: LU LDL T 100,LDL T LU LU n 3 /3 + O(n 2 ) LDL T n 3 /6 + O(n 2 ) 41

43 7.3 C LU LDL T 42

44 8 43

45 , Gauss 3,Gauss LU 4,LDL T Cholesky. 5,LDL T 6, 7, 8.3 A LDL T LU 44

46 8.4 LU LDL T LU n 3 /3 + O(n 2 ) LDL T n 3 /6 + O(n 2 ) LDL T LU 45

47 46

48 ,,..,,. 47

49 [1] :,, 2003, [2] :, 1999 [3] :,, 2002, [4] G.H.Golub and C.F.Van Loan: Matrix Computations,Johns Hopkins,

50 Appendix A 49

51 ///////////////LDL^T_Factorization/////////////// #include <stdio.h> #include <stdlib.h> #include <time.h> #define N 1000 /* */ int main(void); int main(void){ int i,j,k; double s; double t; double sum; static double a[n+1][n+1]; static double b[n+1]; static double w[n+1][n+1]; static double x[n+1]; static double y[n+1]; static double l[n+1][n+1]; static double d[n+1]; static double ld[n+1][n+1]; unsigned long t1,t2 ; /*a_{11 a[1][1] */ /*A */ for(i=1;i<=n;i++){ for(j=i;j<=n;j++){ a[i][j]=j; for(j=i-1;j>=1;j--){ a[i][j]=a[j][i]; /* X */ 50

52 for(i=1;i<=n;i++){ x[i]=i; /*Ax=b B */ for(i=1;i<=n;i++){ for(k=1,sum=0;k<=n;k++){ sum+=a[i][k]*x[k]; b[i] = sum; /* */ for(i=1;i<=n;i++){ x[i]=0; t1=clock(); /* LDL^T */ /*L */ for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ if(i==j){ else{ l[i][j]=1; l[i][j]=0; d[1]=a[1][1]; for(k=2;k<=n;k++){ for(i=1;i<=k-1;i++){ for(j=1,s=0;j<=i-1;j++){ 51

53 s+=l[k][j]*l[i][j]*d[j]; l[k][i]=(a[k][i] - s)/d[i]; for(i=1,t=0;i<=k-1;i++){ t+=l[k][i]*l[k][i]*d[i]; d[k]=a[k][k]-t ; /* LDL^T */ for(i=1;i<=n;i++){ y[i]=b[i]; for(j=1;j<i;j++){ y[i]-=l[i][j]*y[j]; /* DL^T*x=y */ /*DL^T */ for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ for(i=n;i>=1;i--){ w[i][j]=d[i]*l[j][i]; x[i]=y[i]; for(j=n;j>i;j--){ x[i]-=w[i][j]*x[j]; x[i]/=w[i][i]; t2=clock(); /*l[j][i] ^ */ 52

54 printf(" ; %f[ms]\n", (t2-t1)/ ); return(0); ///////////////LU_Factorization/////////////// #include <stdio.h> #include <stdlib.h> #include <time.h> #define N 1000 int main(void); int main(void){ int i,j,k; double s,t,sum; static double a[n+1][n+1]; static double b[n+1]; static double u[n+1][n+1]; static double w[n+1][n+1]; static double x[n+1]; static double y[n+1]; static double l[n+1][n+1]; static double d[n+1]; static double ld[n+1][n+1]; unsigned long t1,t2 ; /*a_{11 a[1][1] */ /*A */ for(i=1;i<=n;i++){ for(j=i;j<=n;j++){ a[i][j]=j; for(j=i-1;j>=1;j--){ a[i][j]=a[j][i]; 53

55 /* X */ for(i=1;i<=n;i++){ x[i]=i; /*Ax=b B */ for(i=1;i<=n;i++){ for(k=1,sum=0;k<=n;k++){ sum+=a[i][k]*x[k]; b[i] = sum; /* */ for(i=1;i<=n;i++){ x[i]=0; t1=clock(); /* LU */ for( i = 1; i <= N; i++ ){ for( j = 1; j <= N; j++ ){ if( i == j ){ l[i][j] = 1; else{ l[i][j] = 0; u[i][j] = a[i][j]; for( i = 1; i <= N; i++ ){ for( j = i+1; j <= N; j++ ){ l[j][i] = u[j][i] / u[i][i]; for( k = i; k <= N; k++ ){ u[j][k] = u[j][k] - u[i][k] * l[j][i]; 54

56 for(i=1;i<=n;i++){ y[i]=b[i]; for(j=1;j<i;j++){ y[i]-=l[i][j]*y[j]; for(i=n;i>=1;i--){ x[i]=y[i]; for(j=n;j>i;j--){ x[i]-=u[i][j]*x[j]; x[i]/=u[i][i]; t2=clock(); printf(" ; %f[ms]\n", (t2-t1)/ ); return(0); 55

Ｃ言語による数値計算プログラミング演習

Ｃ言語による数値計算プログラミング演習 5. 行列の固有値問題 n n 正方行列 A に対する n 個の固有値 λ i (i=1,,,n) と対応する固有ベクトル u i は次式を満たす Au = λ u i i i a11 a1 L a1 n u1i a1 a a n u i A =, ui = M O M M an 1 an L ann uni これらはまとめて, つぎのように書ける 5.1 ヤコビ法 = Λ, = [ u1 u u