aritmitiki-analisi2 - ΑΡΙΘΜΗΤΙΚΗ

Info iconThis preview shows page 1. Sign up to view the full content.

View Full Document Right Arrow Icon
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 ΕΡΩΤΗΜΑ 1: Το πρόγραμμα που ακολουθεί εκτελεί μόνο του όλα τα ζητούμενα και τα αποτελέσματα αποθηκεύονται σε ένα αρχείο. Ένα πρώτο σχόλιο είναι ότι χρησιμοποιούνται τέσσερις επιπλέον πίνακες εκτός από τους πίνακες που χρησιμοποιεί το πρόγραμμα DGAUSS της δισκέτας. Ένας πίνακας είναι ο δισδιάστατος πίνακας G, στον οποίο αποθηκεύεται στήλη στήλη το κάθε αποτέλεσμα της μεθόδου για κάθε διαφορετικό πίνακα Β που εισάγεται σε αυτήν με μια αναδρομή. Ο δεύτερος πίνακας είναι και αυτός δισδιάστατος, ο C, ο οποίος προκύπτει από τον πολλαπλασιασμό του Α με τον αντίστροφό του. Μετά αφαιρούμε από κάθε στοιχείο της διαγωνίου τη μονάδα, για να προκύψει ο ζητούμενος πίνακας υπόλοιπο. Οι άλλοι δυο πίνακες Υ και LPIV είναι πίνακες που χρησιμοποιούνται για την εύρεση του δείκτη κατάστασης (σύμφωνα με το πρόγραμμα CONDE της δισκέτας ). Ένα δεύτερο σχόλιο είναι ότι χρησιμοποιείται το IOPT ως μετρητής για να αποφύγω την τριγωνοποίηση του πίνακα Α για περισσότερες φορές από μία. Για το λόγο αυτό υπάρχει διαφορετική συνθήκη στην υπορουτίνα DGAUSS για το ΙΟΡΤ. Ένα τρίτο σχόλιο είναι ότι το πρόγραμμα με ένα επαναληπτικό βρόγχο για n=3,..,8 εκτελεί τα ζητούμενα, εκτιμώντας επιπλέον το δείκτη κατάστασης (CONDEST). Ακολουθεί ο κώδικας, και στη συνέχεια το αρχείο στο οποίο αποθηκεύτηκαν τα αποτελέσματα του προγράμματος. (Επειδή κατά το COPY-PASTE από το Digital V.Fortran στο Word τα κενά παρουσιάζονται κάπως διαφοροποιημένα, δηλαδή άλλες γραμμές είναι πιο έξω από κάποιες άλλες κλπ. παρέχεται δισκέτα με τον κώδικα των προγραμμάτων για διευκόλυνση στην ανάγνωση). C PROGRAM ASK1 C C PARAMETER(ND=10) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION A,B,C,G,D DIMENSION A(ND,ND),B(ND),L(ND),C(ND,ND),G(ND,ND),Y(ND),LPIV(ND) OPEN(UNIT=6,FILE='ASK1.TXT') DO 120 N=3,8 WRITE(6,*)'N=',N DO 6 I=1,N DO 6 J=1,N A(I,J)=1.D0/(I+J-1.D0) 6 CONTINUE IOPT=0 DO 7 M=1,N DO 8 K=1,N B(K)=0.D0 8 CONTINUE B(M)=1 CALL DGAUSS(ND,N,A,B,IOPT,D,IREG,L) IOPT=IOPT+1 DO 9 I=1,N 1 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 G(I,M)=B(I) 9 CONTINUE 7 CONTINUE WRITE(6,*)'OUTPUT:' WRITE(6,*)'A^-1:' DO 11 I=1,N WRITE(6,*) (G(I,J),J=1,N) 11 WRITE(6,*)' ' DO 12 I=1,N DO 12 J=1,N A(I,J)=1.D0/(I+J-1.D0) 12 CONTINUE DO 13 I=1,N DO 13 J=1,N C(I,J)=0 DO 13 K=1,N C(I,J)=C(I,J)+A(I,K)*G(K,J) 13 CONTINUE DO 14 I=1,N C(I,I)=C(I,I)-1.D0 14 CONTINUE WRITE(6,*)'R=AX-I:' DO 15 I=1,N WRITE(6,*) (C(I,J),J=1,N) 15 WRITE(6,*)' ' C COMPUTE NORM1 OF A C AN=0.D0 DO 20 J=1,N S=0.D0 DO 30 I=1,N 30 S=S+DABS(A(I,J)) IF(S.GT.AN) AN=S 20 CONTINUE C C TRIANGULAR DECOMPOSITION OF A=L*U, WHERE L AND U ARE 'STORED' C IN THE LOWER AND UPPER-TRIANGULAR PART OF A, BY GAUSS' METHOD C CALL DGAUSSCONDE(ND,N,A,Y,1,D,IREG,LPIV) C IF(IREG.EQ.0) STOP C C ESTIMATE THE CONDITION NUMBER OF A : C - SOLVE THE SYSTEM AT*Y=E BY SOLVING SUCCESSIVELY THE C LOWER AND UPPER-TRIANGULAR SYSTEMS UT*Y'=E AND LT*Y=Y'. 2 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C C C C C C C C E IS A VECTOR WITH COMPONENTS +1 OR -1 CHOSEN SO AS TO MAXIMIZE THE COMPONENTS OF Y' IN ABSOLUTE VALUE. - SOLVE THE SYSTEM U*Z=Y. - COMPUTE THE ESTIMATE CONDEST=(NORM1 OF A)*(NORM1 OF Z)/(NORM1 OF Y). SOLVE THE SYSTEM UT*Y'=E BY FORWORD-SUBSTITUTION DO 50 K=1,N S=0.D0 IF(K.EQ.1) GO TO 45 KM=K-1 DO 40 J=1,KM S=S+A(J,K)*Y(J) 40 CONTINUE 45 E=-1.D0 IF(S.LT.0.D0) E=1.D0 Y(K)=(E-S)/A(K,K) 50 CONTINUE C C SOLVE THE SYSTEM LT*Y=Y' BY BACK-SUBSTITUTION (THE DIAGONAL C ELEMENTS OF L ARE EQUAL TO 1) AND RE-ORDER THE COMPONENTS OF Y. C DO 60 I=1,N-1 K=N-I S=0.D0 KP=K+1 DO 55 J=KP,N S=S+A(J,K)*Y(J) 55 CONTINUE Y(K)=Y(K)-S C M=LPIV(K) IF(M.EQ.K) GO TO 60 S=Y(M) Y(M)=Y(K) Y(K)=S 60 CONTINUE C C COMPUTE NORM1 OF Y C YN=0.D0 DO 65 I=1,N YN=YN+DABS(Y(I)) 65 CONTINUE C 3 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C C SOLVE THE SYSTEM U*Z=Y BY BACH-SUBSTITUTION CALL DGAUSSCONDE(ND,N,A,Y,2,D,IREG,LPIV) C C C COMPUTE NORM1 OF Z ZN=0.D0 DO 70 I=1,N ZN=ZN+DABS(Y(I)) 70 CONTINUE C C COMPUTE THE ESTIMATE OF THE CONDITION NUMBER C CE=AN*ZN/YN IF(CE.LT.1.D0) CE=1.D0 WRITE(6,*) 'CONDEST=',CE 120 CONTINUE STOP END C------------------------------------------------------------------SUBROUTINE DGAUSS(ND,N,A,B,IOPT,D,IREG,L) C C DECOMPOSITION OF THE MATRIX A AND/OR SOLUTION OF THE C LINEAR SYSTEM AX=B BY GAUSS' ELIMINATION METHOD C C INPUT: C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS C N: DIMENSION OF ARRAYS C A: SYSTEM MATRIX (NXN) C B: RIGHT-HAND SIDE VECTOR (N) C IOPT = 0, DECOMPOSITION OF A AND SOLUTION OF THE SYSTEM C = 1, DECOMPOSITION OF A ONLY C = 2, SOLUTION OF THE SYSTEM (SKIPS DECOMPOSITION) C C OUTPUT: C A: DECOMPOSED MATRIX (LOWER\UPPER PART), IF IOPT=0 OR 1 C B:=X: SOLUTION VECTOR STORED IN B (N) C D: DETERMINANT OF A C IREG = 0, SINGULAR MATRIX A C = 1, REGULAR MATRIX A C L: VECTOR OF PIVOTING INDICES (N) C DOUBLE PRECISION A,B,D,S,EE DIMENSION A(ND,N),B(ND),L(ND) EE=1.D-12 IREG=1 NM=N-1 4 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 IF(IOPT.NE.0) GO TO 60 C C TRIANGULAR DECOMPOSITION OF A WITH MAXIMUM COLUMN PIVOTING C (ELIMINATION WITH MULTIPLIERS STORED IN THE LOWER PART OF A) C D=1.D0 DO 10 K=1,NM KP=K+1 M=K DO 20 I=KP,N IF(DABS(A(I,K)).GT.DABS(A(M,K))) M=I 20 CONTINUE L(K)=M D=D*A(M,K) IF(DABS(A(M,K)).LE.EE) GO TO 100 IF(M.EQ.K) GO TO 30 D=-D DO 40 J=K,N S=A(M,J) A(M,J)=A(K,J) 40 A(K,J)=S 30 DO 50 I=KP,N S=A(I,K)/A(K,K) A(I,K)=S DO 50 J=KP,N 50 A(I,J)=A(I,J)-S*A(K,J) 10 CONTINUE D=D*A(N,N) IF(DABS(A(N,N)).LE.EE) GO TO 100 IF(IOPT.EQ.1) RETURN C C ELIMINATION OPERATIONS ON B C 60 DO 70 K=1,NM KP=K+1 M=L(K) S=B(M) B(M)=B(K) B(K)=S DO 70 I=KP,N 70 B(I)=B(I)-A(I,K)*S C C BACK-SUBSTITUTION C B(N)=B(N)/A(N,N) DO 80 I=1,NM 5 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 K=N-I KP=K+1 S=0.D0 DO 90 J=KP,N 90 S=S+A(K,J)*B(J) 80 B(K)=(B(K)-S)/A(K,K) RETURN 100 IREG=0 RETURN END C--------------------------------------------------------------------------SUBROUTINE DGAUSSCONDE(ND,N,A,B,IOPT,D,IREG,LPIV) C C DECOMPOSITION OF THE MATRIX A AND/OR SOLUTION OF THE C LINEAR SYSTEM A*X=B BY GAUSS' ELIMINATION METHOD C C INPUT: C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS C N: DIMENSION OF ARRAYS C A: SYSTEM MATRIX (NxN) C B: RIGHT-HAND SIDE VECTOR (N) C IOPT = 0, DECOMPOSITION OF A AND SOLUTION OF THE SYSTEM C = 1, DECOMPOSITION OF A ONLY C = 2, SOLUTION OF THE SYSTEM (SKIPS DECOMPOSITION) C C OUTPUT: C A: DECOMPOSED MATRIX (LOWER\UPPER PART), IF IOPT=0 OR 1 C B:=X: SOLUTION VECTOR STORED IN B (N) C D: DETERMINANT OF A C IREG = 0, SINGULAR MATRIX A C = 1, REGULAR MATRIX A C LPIV: VECTOR OF PIVOTING INDICES (N) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(ND,N),B(ND),LPIV(ND) EE=1.D-14 IREG=1 NM=N-1 IF(IOPT.EQ.2) GO TO 60 C C TRIANGULAR DECOMPOSITION OF A WITH MAXIMUM COLUMN PIVOTING C (ELIMINATION WITH MULTIPLIERS STORED IN THE LOWER PART OF A) C D=1.D0 DO 10 K=1,NM KP=K+1 6 ΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 M=K DO 20 I=KP,N IF(DABS(A(I,K)).GT.DABS(A(M,K))) M=I 20 CONTINUE LPIV(K)=M D=D*A(M,K) IF(DABS(A(M,K)).LE.EE) GO TO 100 IF(M.EQ.K) GO TO 30 D=-D DO 40 J=K,N S=A(M,J) A(M,J)=A(K,J) 40 A(K,J)=S 30 DO 50 I=KP,N S=A(I,K)/A(K,K) A(I,K)=S DO 50 J=KP,N 50 A(I,J)=A(I,J)-S*A(K,J) 10 CONTINUE D=D*A(N,N) IF(DABS(A(N,N)).LE.EE) GO TO 100 IF(IOPT.EQ.1) RETURN C C C ELIMINATION OPERATIONS ON B 60 DO 70 K=1,NM KP=K+1 M=LPIV(K) S=B(M) B(M)=B(K) B(K)=S DO 70 I=KP,N 70 B(I)=B(I)-A(I,K)*S C C C BACK-SUBSTITUTION B(N)=B(N)/A(N,N) DO 80 I=1,NM K=N-I KP=K+1 S=0.D0 DO 90 J=KP,N 90 S=S+A(K,J)*B(J) 80 B(K)=(B(K)-S)/A(K,K) RETURN 100 IREG=0 RETURN 7 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 END Ακολουθεί το αρχείο με το πρόγραμμα εκτελεσμένο. N= 3 OUTPUT: A^-1: 9.00000000000003 -36.0000000000001 30.0000000000001 -36.0000000000001 192.000000000001 -180.000000000001 30.0000000000001 -180.000000000001 180.000000000001 R=AX-I: 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 8.881784197001252E-016 -7.105427357601002E-015 0.000000000000000E+000 0.000000000000000E+000 -7.105427357601002E-015 0.000000000000000E+000 CONDEST= 680.809278350518 N= 4 OUTPUT: A^-1: 15.9999999999993 -119.999999999992 -139.999999999986 239.999999999979 -119.999999999992 1679.99999999984 1199.99999999990 -2699.99999999976 239.999999999980 -4199.99999999961 -2699.99999999976 6479.99999999940 -139.999999999987 2799.99999999974 1679.99999999984 -4199.99999999961 R=AX-I: 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 3.552713678800501E-015 -5.684341886080801E-014 -1.136868377216160E-013 0.000000000000000E+000 0.000000000000000E+000 -5.684341886080801E-014 0.000000000000000E+000 -1.136868377216160E-013 8 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 0.000000000000000E+000 2.842170943040401E-014 1.136868377216160E-013 5.684341886080801E-014 CONDEST= 21523.4924623096 N= 5 OUTPUT: A^-1: 24.9999999999919 -299.999999999882 -1399.99999999946 629.999999999774 1049.99999999958 -299.999999999879 26879.9999999943 4799.99999999840 -12599.9999999979 -18899.9999999950 1049.99999999955 -117599.999999987 -18899.9999999947 56699.9999999967 79379.9999999856 -1399.99999999941 179199.999999992 26879.9999999937 -88200.0000000006 -117599.999999986 629.999999999739 -88199.9999999999 -12599.9999999975 44100.0000000023 56699.9999999955 R=AX-I: 0.000000000000000E+000 -4.547473508864641E-013 -1.818989403545856E-012 0.000000000000000E+000 0.000000000000000E+000 -1.421085471520200E-014 0.000000000000000E+000 0.000000000000000E+000 -1.818989403545856E-012 -2.728484105318785E-012 0.000000000000000E+000 -4.547473508864641E-013 -3.637978807091713E-012 0.000000000000000E+000 -9.094947017729282E-013 0.000000000000000E+000 2.273736754432321E-013 0.000000000000000E+000 0.000000000000000E+000 -9.094947017729282E-013 1.421085471520200E-014 0.000000000000000E+000 0.000000000000000E+000 -1.818989403545856E-012 0.000000000000000E+000 CONDEST= 694086.551600326 N= 6 OUTPUT: A^-1: 36.0000000012132 -630.000000036420 -7560.00000067314 7560.00000075274 3360.00000025373 -2772.00000029968 -630.000000035890 -88200.0000074889 14700.0000010760 9 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 211680.000019846 -220500.000022196 83160.0000088327 3360.00000024775 -1411200.00013675 -88200.0000074209 1512000.00015291 564480.000051622 -582120.000060838 -7560.00000065281 3628800.00035992 211680.000019543 -3969000.00040240 -1411200.00013590 1552320.00016008 7560.00000072711 -3969000.00040057 -220500.000021758 4410000.00044780 1512000.00015127 -1746360.00017813 -2772.00000028842 1552320.00015880 83160.0000086282 -1746360.00017751 -582120.000059976 698544.000070610 R=AX-I: 1.136868377216160E-013 0.000000000000000E+000 0.000000000000000E+000 2.910383045673370E-011 0.000000000000000E+000 0.000000000000000E+000 1.136868377216160E-013 1.818989403545856E-012 0.000000000000000E+000 1.455191522836685E-010 5.820766091346741E-011 0.000000000000000E+000 -5.684341886080801E-014 0.000000000000000E+000 2.910383045673370E-011 8.731149137020111E-011 -2.910383045673370E-011 0.000000000000000E+000 -1.136868377216160E-013 -1.818989403545856E-012 7.275957614183426E-012 8.731149137020111E-011 2.910383045673370E-011 1.455191522836685E-011 1.136868377216160E-013 -1.818989403545856E-012 2.910383045673370E-011 0.000000000000000E+000 2.910383045673370E-011 -1.455191522836685E-011 0.000000000000000E+000 0.000000000000000E+000 7.275957614183426E-012 2.910383045673370E-011 0.000000000000000E+000 -1.455191522836685E-011 CONDEST= 22618622.5955126 N= 7 OUTPUT: A^-1: 49.0000000598698 -1176.00000240191 -29400.0000902577 48510.0001655510 12012.0000469736 8820.00002319296 -38808.0001430726 -1176.00000238850 1128960.00359705 -504504.001871147 37632.0000957828 -1940400.00659609 -317520.000924539 1596672.00569963 8820.00002298274 -10584000.0345871 -317520.000921351 18711000.0634155 2857680.00889129 -15717240.0547908 10 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 5045040.01798604 -29400.0000892284 40320000.1342113 -20180160.0697760 1128960.00357622 -72765000.2460521 -10584000.0345058 62092800.2125720 48510.0001633697 -72765000.2456298 37837800.1276783 -1940400.00654657 133402500.450284 18711000.0631577 -115259760.388990 -38808.0001410007 62092800.2119284 -33297264.1101437 1596672.00564938 -115259760.388479 -15717240.0544964 100590336.335582 12012.0000462473 -20180160.0694924 11099088.0361122 -504504.001852735 37837800.1273779 5045040.01787080 -33297264.1100290 R=AX-I: -4.547473508864641E-013 0.000000000000000E+000 2.328306436538696E-010 4.656612873077393E-010 -9.313225746154785E-010 9.313225746154785E-010 -4.656612873077393E-010 0.000000000000000E+000 2.182787284255028E-011 2.328306436538696E-010 -2.328306436538696E-009 2.793967723846436E-009 -4.656612873077393E-010 4.656612873077393E-010 0.000000000000000E+000 -3.637978807091713E-011 0.000000000000000E+000 -9.313225746154785E-010 1.862645149230957E-009 1.862645149230957E-009 -2.328306436538696E-010 2.273736754432321E-013 -1.455191522836685E-011 5.820766091346741E-010 -6.984919309616089E-010 -4.656612873077393E-010 -4.656612873077393E-010 2.328306436538696E-010 2.273736754432321E-013 -7.275957614183426E-012 2.910383045673370E-010 -1.396983861923218E-009 9.313225746154785E-010 -1.396983861923218E-009 2.328306436538696E-010 0.000000000000000E+000 -2.182787284255028E-011 0.000000000000000E+000 0.000000000000000E+000 -9.313225746154785E-010 0.000000000000000E+000 -3.492459654808044E-010 -9.094947017729282E-013 -7.275957614183426E-012 -5.820766091346741E-011 -6.984919309616089E-010 -4.656612873077393E-010 0.000000000000000E+000 -3.492459654808044E-010 11 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 CONDEST= 741927130.071823 N= 8 OUTPUT: A^-1: 64.0000001567259 -2016.00000531168 -92400.0000847057 221759.999954030 192191.999539107 -51479.9998248257 20160.0000391360 -288287.999617279 -2016.00000501411 4656959.99159789 -10594583.9400957 84672.0001023525 -11642399.9626484 2882879.98067319 -952559.999710740 15567551.9300742 20160.0000326871 -58211999.7651100 141261118.877496 -952559.999572029 149687999.174872 -38918879.6515659 11430719.9731289 -204324118.609099 -92400.0000395089 304919998.204538 -776936152.500370 4656959.99029175 -800414994.121587 216215997.712116 -58211999.7578406 1109908790.47180 221759.999810831 -800414994.038503 2118916776.74468 -11642399.9579623 2134439981.13480 -594593992.975400 149687999.139968 -2996753730.04424 -288287.999390060 1109908790.25877 -3030050987.54586 15567551.9221033 -2996753729.78701 856215349.058711 -204324118.542110 4249941648.63007 192191.999362710 -776936152.292235 2175421219.95031 -10594583.9336307 2118916776.40922 -618377751.618841 141261118.819409 -3030050987.33009 -51479.9997715968 216215997.639875 -618377751.582082 2882879.97865969 -594593992.844291 176679357.493075 -38918879.6326882 856215348.946199 R=AX-I: -1.818989403545856E-012 -1.164153218269348E-010 -1.862645149230957E-009 7.450580596923828E-009 1.490116119384766E-008 -5.960464477539063E-008 2.980232238769531E-008 3.725290298461914E-009 -9.094947017729282E-013 1.164153218269348E-010 0.000000000000000E+000 1.490116119384766E-008 -4.470348358154297E-008 -4.470348358154297E-008 2.980232238769531E-008 0.000000000000000E+000 -3.637978807091713E-012 1.746229827404022E-010 0.000000000000000E+000 7.450580596923828E-009 -7.450580596923828E-009 -4.470348358154297E-008 12 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 -1.490116119384766E-008 3.725290298461914E-009 1.818989403545856E-012 -8.731149137020111E-011 1.396983861923218E-009 -7.450580596923828E-009 3.725290298461914E-008 -5.960464477539063E-008 1.490116119384766E-008 -1.862645149230957E-008 1.818989403545856E-012 2.910383045673370E-011 1.862645149230957E-009 0.000000000000000E+000 4.470348358154297E-008 0.000000000000000E+000 -7.450580596923828E-009 -1.117587089538574E-008 0.000000000000000E+000 1.746229827404022E-010 -4.190951585769653E-009 1.117587089538574E-008 -7.450580596923828E-009 -5.960464477539063E-008 1.490116119384766E-008 1.862645149230957E-009 -3.637978807091713E-012 1.455191522836685E-010 -1.862645149230957E-009 0.000000000000000E+000 -7.450580596923828E-009 2.980232238769531E-008 -2.980232238769531E-008 5.587935447692871E-009 1.818989403545856E-012 1.164153218269348E-010 9.313225746154785E-010 1.862645149230957E-009 2.980232238769531E-008 -1.490116119384766E-008 -7.450580596923828E-009 -1.117587089538574E-008 CONDEST= 24446819273.4855 ΕΡΩΤΗΜΑ 2: Ένα πρώτο σχόλιο για το πρόγραμμα αυτού του δεύτερου ερωτήματος είναι ότι χρησιμοποιούνται τρεις επιπλέον πίνακες εκτός από αυτούς του προγράμματος DGAUSS. Ένας πίνακας είναι ο μονοδιάστατος Β1 που τον χρησιμοποιούμε για να κατασκευάσουμε τον πίνακα Β (ο Β προκύπτει με πολλαπλασιασμό του Α με το διάνυσμα Β1(1 1 1....). Ο δεύτερος πίνακας είναι ο C, ένας μονοδιάστατος πίνακας στον οποίο αποθηκεύονται οι νόρμες των γραμμών του αρχικού πίνακα. Ο τρίτος μονοδιάστατος πίνακας είναι ο πίνακας G, ο οποίος χρησιμοποιείται για την εύρεση της νόρμας του σφάλματος. Ένα δεύτερο σχόλιο που πρέπει να γίνει αφορά το οδηγό στοιχείο που χρησιμοποιούμε το οποίο τώρα δίνεται από τον τύπο amk=max(i/c). Για το λόγο αυτό 1≤in χρησιμοποιείται διαφορετικό κριτήριο εύρεσης οδηγού στοιχείου στην υπορουτίνα από το πρόγραμμα DGAUSS. Ακολουθεί ο κώδικας του ζητούμενου προγράμματος. C C C PROGRAM ASK2 PARAMETER(ND=10) DOUBLE PRECISION A,B,D,B1,G,NORM,C DIMENSION A(ND,ND),B(ND),L(ND),B1(ND),G(ND),C(ND) READ(5,*) N,IOPT 13 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 DO 6 I=1,N DO 6 J=1,N A(I,J)=1.D0/(I+J-1.D0) 6 CONTINUE DO 8 I=1,N B1(I)=1 8 CONTINUE DO 9 I=1,N B(I)=0 DO 9 K=1,N B(I)=B(I)+A(I,K)*B1(K) 9 CONTINUE DO 4 K=1,N C(K)=A(K,1) DO 5 J=1,N-1 IF(DABS(A(K,J+1)).GT.DABS(C(K))) C(K)=A(K,J+1) 5 CONTINUE 4 CONTINUE WRITE(6,*) ‘INPUT:’ WRITE(6,*) ‘N=’,N,’ IOPT=’,IOPT WRITE(6,*) ‘A=’ DO 10 I=1,N WRITE(6,*) (A(I,J),J=1,N) 10 WRITE(6,*) ‘ ‘ WRITE(6,*) ‘B=’,(B(I),I=1,N) WRITE(6,*) ‘ ‘ WRITE(6,*) ‘OUTPUT:’ CALL DGAUSS1(ND,N,A,B,C,IOPT,D,IREG,L) WRITE(6,*) ‘D=’,D,’ IREG=’,IREG IF(IREG.EQ.0) STOP WRITE(6,*) ‘X=’,(B(I),I=1,N) DO 21 K=1,N G(K)=B(K)-1 21 CONTINUE NORM=DABS(G(1)) DO 22 K=1,N-1 IF(DABS(G(K+1)).GT.NORM) NORM=DABS(G(K+1)) 22 CONTINUE WRITE(6,*) ‘NORM X-Z=’,NORM STOP END C------------------------------------------------------------------SUBROUTINE DGAUSS1(ND,N,A,B,C,IOPT,D,IREG,L) C C DECOMPOSITION OF THE MATRIX A AND/OR SOLUTION OF THE C LINEAR SYSTEM AX=B BY GAUSS’ ELIMINATION METHOD C 14 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C C C C C C C C C C C C C C C C C INPUT: ND: MAXIMUM COLUMN DIMENSION OF ARRAYS N: DIMENSION OF ARRAYS A: SYSTEM MATRIX (NXN) B: RIGHT-HAND SIDE VECTOR (N) IOPT = 0, DECOMPOSITION OF A AND SOLUTION OF THE SYSTEM = 1, DECOMPOSITION OF A ONLY = 2, SOLUTION OF THE SYSTEM (SKIPS DECOMPOSITION) OUTPUT: A: DECOMPOSED MATRIX (LOWER\UPPER PART), IF IOPT=0 OR 1 B:=X: SOLUTION VECTOR STORED IN B (N) D: DETERMINANT OF A IREG = 0, SINGULAR MATRIX A = 1, REGULAR MATRIX A L: VECTOR OF PIVOTING INDICES (N) DOUBLE PRECISION A,B,D,S,EE,C DIMENSION A(ND,N),B(ND),L(ND),C(ND) EE=1.D-12 IREG=1 NM=N-1 IF(IOPT.EQ.2) GO TO 60 C C TRIANGULAR DECOMPOSITION OF A WITH MAXIMUM COLUMN PIVOTING C (ELIMINATION WITH MULTIPLIERS STORED IN THE LOWER PART OF A) C D=1.D0 DO 10 K=1,NM KP=K+1 M=K DO 20 I=KP,N IF(DABS(A(I,K)/C(I)).GT.DABS(A(M,K)/C(M))) M=I 100 CONTINUE L(K)=M D=D*A(M,K) IF(DABS(A(M,K)).LE.EE) GO TO 100 IF(M.EQ.K) GO TO 30 D=-D DO 40 J=K,N S=A(M,J) A(M,J)=A(K,J) 40 A(K,J)=S 30 DO 50 I=KP,N S=A(I,K)/A(K,K) A(I,K)=S 15 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 DO 50 J=KP,N 50 A(I,J)=A(I,J)-S*A(K,J) 10 CONTINUE D=D*A(N,N) IF(DABS(A(N,N)).LE.EE) GO TO 100 IF(IOPT.EQ.1) RETURN C C C ELIMINATION OPERATIONS ON B 60 DO 70 K=1,NM KP=K+1 M=L(K) S=B(M) B(M)=B(K) B(K)=S DO 70 I=KP,N 70 B(I)=B(I)-A(I,K)*S C C C BACK-SUBSTITUTION B(N)=B(N)/A(N,N) DO 80 I=1,NM K=N-I KP=K+1 S=0.D0 DO 90 J=KP,N 90 S=S+A(K,J)*B(J) 80 B(K)=(B(K)-S)/A(K,K) RETURN 100 IREG=0 RETURN END Ακολουθεί το εκτελεσμένο πρόγραμμα το οποίο υπολογίζει τη νόρμα του σφάλματος x − z ∞ του προγράμματος αυτού. 7 0 INPUT: N= 7 IOPT= A= 1.00000000000000 0.250000000000000 0.142857142857143 0.500000000000000 0.200000000000000 0 0.500000000000000 0.200000000000000 0.333333333333333 0.166666666666667 0.333333333333333 0.166666666666667 0.250000000000000 0.142857142857143 16 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 0.125000000000000 0.333333333333333 0.166666666666667 0.111111111111111 0.250000000000000 0.142857142857143 0.200000000000000 0.125000000000000 0.250000000000000 0.142857142857143 0.100000000000000 0.200000000000000 0.125000000000000 0.166666666666667 0.111111111111111 0.200000000000000 0.166666666666667 0.125000000000000 0.111111111111111 9.090909090909091E-002 0.142857142857143 0.100000000000000 0.166666666666667 0.142857142857143 0.111111111111111 0.100000000000000 8.333333333333333E-002 0.125000000000000 9.090909090909091E-002 0.142857142857143 0.125000000000000 0.111111111111111 0.100000000000000 9.090909090909091E-002 8.333333333333333E-002 7.692307692307693E-002 B= 2.59285714285714 1.71785714285714 1.32896825396825 1.09563492063492 0.936544011544012 0.819877344877345 0.730133755133755 OUTPUT: D= 4.835802605930008E-025 IREG= 1 X= 0.999999999994110 1.00000000023215 0.999999997781628 1.00000000857461 0.999999984348117 1.00000001347937 0.999999995586050 NORM X-Z= 1.565188267882434E-008 Ακολουθεί ο κώδικας του προγράμματος DGAUSS της δισκέτας με τη διαφορά ότι στο κυρίως πρόγραμμα υπάρχει ένας επιπλέον πίνακας G που χρησιμοποιείται για την εύρεση της νόρμας x − z ∞ της μεθόδου GAUSS. C C C C PROGRAM DGAUS GAUSS' ELIMINATION METHOD PARAMETER(ND=10) DOUBLE PRECISION A,B,B1,D,NORM,G DIMENSION A(ND,ND),B(ND),L(ND),G(ND),B1(ND) READ(5,*) N,IOPT DO 6 I=1,N DO 6 J=1,N 17 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 A(I,J)=1.D0/(I+J-1.D0) 6 CONTINUE DO 8 I=1,N B1(I)=1 8 CONTINUE DO 9 I=1,N B(I)=0 DO 9 K=1,N B(I)=B(I)+A(I,K)*B1(K) 9 CONTINUE WRITE(6,*) 'INPUT:' WRITE(6,*) 'N=',N,' IOPT=',IOPT WRITE(6,*) 'A=' DO 10 I=1,N WRITE(6,*) (A(I,J),J=1,N) 10 WRITE(6,*) ' ' WRITE(6,*) 'B=',(B(I),I=1,N) WRITE(6,*) ' ' WRITE(6,*) 'OUTPUT:' CALL DGAUSS(ND,N,A,B,IOPT,D,IREG,L) WRITE(6,*) 'D=',D,' IREG=',IREG IF(IREG.EQ.0) STOP WRITE(6,*) 'X=',(B(I),I=1,N) DO 21 K=1,N G(K)=B(K)-1 21 CONTINUE NORM=DABS(G(1)) DO 22 K=1,N-1 IF(DABS(G(K+1)).GT.NORM) NORM=DABS(G(K+1)) 22 CONTINUE WRITE(6,*) 'NORM X-Z=',NORM STOP END C------------------------------------------------------------------SUBROUTINE DGAUSS(ND,N,A,B,IOPT,D,IREG,L) C C DECOMPOSITION OF THE MATRIX A AND/OR SOLUTION OF THE C LINEAR SYSTEM AX=B BY GAUSS' ELIMINATION METHOD C C INPUT: C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS C N: DIMENSION OF ARRAYS C A: SYSTEM MATRIX (NXN) C B: RIGHT-HAND SIDE VECTOR (N) C IOPT = 0, DECOMPOSITION OF A AND SOLUTION OF THE SYSTEM C = 1, DECOMPOSITION OF A ONLY C = 2, SOLUTION OF THE SYSTEM (SKIPS DECOMPOSITION) 18 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C C C C C C C C C OUTPUT: A: DECOMPOSED MATRIX (LOWER\UPPER PART), IF IOPT=0 OR 1 B:=X: SOLUTION VECTOR STORED IN B (N) D: DETERMINANT OF A IREG = 0, SINGULAR MATRIX A = 1, REGULAR MATRIX A L: VECTOR OF PIVOTING INDICES (N) DOUBLE PRECISION A,B,D,S,EE DIMENSION A(ND,N),B(ND),L(ND) EE=1.D-14 IREG=1 NM=N-1 IF(IOPT.EQ.2) GO TO 60 C C TRIANGULAR DECOMPOSITION OF A WITH MAXIMUM COLUMN PIVOTING C (ELIMINATION WITH MULTIPLIERS STORED IN THE LOWER PART OF A) C D=1.D0 DO 10 K=1,NM KP=K+1 M=K DO 20 I=KP,N IF(DABS(A(I,K)).GT.DABS(A(M,K))) M=I 20 CONTINUE L(K)=M D=D*A(M,K) IF(DABS(A(M,K)).LE.EE) GO TO 100 IF(M.EQ.K) GO TO 30 D=-D DO 40 J=K,N S=A(M,J) A(M,J)=A(K,J) 40 A(K,J)=S 30 DO 50 I=KP,N S=A(I,K)/A(K,K) A(I,K)=S DO 50 J=KP,N 50 A(I,J)=A(I,J)-S*A(K,J) 10 CONTINUE D=D*A(N,N) IF(DABS(A(N,N)).LE.EE) GO TO 100 IF(IOPT.EQ.1) RETURN C C ELIMINATION OPERATIONS ON B 19 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C 60 DO 70 K=1,NM KP=K+1 M=L(K) S=B(M) B(M)=B(K) B(K)=S DO 70 I=KP,N 70 B(I)=B(I)-A(I,K)*S C C C BACK-SUBSTITUTION B(N)=B(N)/A(N,N) DO 80 I=1,NM K=N-I KP=K+1 S=0.D0 DO 90 J=KP,N 90 S=S+A(K,J)*B(J) 80 B(K)=(B(K)-S)/A(K,K) RETURN 100 IREG=0 RETURN END Aκολουθεί στην επόμενη σελίδα το εκτελεσμένο εκτελεσμένο πρόγραμμα το οποίο υπολογίζει τη νόρμα του σφάλματος x − z ∞ της κλασικής μεθόδου Gauss. 7 0 INPUT: N= 7 IOPT= A= 1.00000000000000 0.250000000000000 0.142857142857143 0 0.500000000000000 0.200000000000000 0.333333333333333 0.166666666666667 0.500000000000000 0.200000000000000 0.125000000000000 0.333333333333333 0.166666666666667 0.250000000000000 0.142857142857143 0.333333333333333 0.166666666666667 0.111111111111111 0.250000000000000 0.142857142857143 0.200000000000000 0.125000000000000 20 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 0.250000000000000 0.142857142857143 0.100000000000000 0.200000000000000 0.125000000000000 0.166666666666667 0.111111111111111 0.200000000000000 0.166666666666667 0.125000000000000 0.111111111111111 9.090909090909091E-002 0.142857142857143 0.100000000000000 0.166666666666667 0.142857142857143 0.111111111111111 0.100000000000000 8.333333333333333E-002 0.125000000000000 9.090909090909091E-002 0.142857142857143 0.125000000000000 0.111111111111111 0.100000000000000 9.090909090909091E-002 8.333333333333333E-002 7.692307692307693E-002 B= 2.59285714285714 1.71785714285714 1.32896825396825 1.09563492063492 0.936544011544012 0.819877344877345 0.730133755133755 OUTPUT: D= 4.835802607684246E-025 IREG= 1 X= 0.999999999990736 1.00000000037119 0.999999996413681 1.00000001397250 0.999999974341236 1.00000002220175 0.999999992702105 NORM X-Z= 2.565876389137856E-008 Παρατηρούμε λοιπόν ότι η νόρμα του σφάλματος x − z ∞ της κλασικής μεόδου Gauss είναι σχεδόν η διπλάσια της νόρμας του σφάλματος της μεθόδου που χρησιμοποιήσαμε. Αυτό ήταν αναμενόμενο αφού η μέθοδος που χρησιμοποιήσαμε δεν χρησιμοποιεί ως οδηγό στοιχείο το μεγαλύτερο κατά απόλυτη τιμή στοιχείο κάθε στήλης του πίνακα, αλλά το στοιχείο με τη μεγαλύτερη δύναμη, το στοιχείο που δίνεται δηλαδή από τον τύπο amk=max(i/c). Συμπεραίνουμε δηλαδή ότι η μέθοδος που χρησιμοποιήσαμε αποτελεί βελτίωση της 1≤in κλασικής μεθόδου του Gauss όσον αφορά την οδήγηση. ΕΡΩΤΗΜΑ 3: Το ερώτημα αυτό ζητά να συγκρίνουμε τις μεθόδους Jacobi και Gauss-Seidel ως προς την ταχύτητα σύγκλισης. Τα προγράμματα που χρησιμοποιούμε είναι αυτά της δισκέτας, με τη διαφορά ότι ο Β κατασκευάζεται μέσα στο πρόγραμμα και ότι τα δεδομένα διαβάζονται από αρχεία. Ένα αρχικό σχόλιο είναι ότι χρησιμοποιήσαμε ως ΕΡS=10 −15 , γιατί για μικρότερο η μέθοδος Gauss-Seidel σταματούσε σε επανάληψη μικρότερη της 21 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 δέκατης πέμπτης και έτσι δεν θα μπορούσαμε να συκρίνουμε τις νόρμες x15 ∞ , όπως εζητείτο. Ακολουθεί ο κώδικας της μεθόδου Jacobi. C C C C PROGRAM JACOB JACOBI'S METHOD PARAMETER(ND=10) DOUBLE PRECISION EPS,A,B,X,Y,E DIMENSION A(ND,ND),B(ND),X(ND),Y(ND) OPEN(UNIT=5,FILE='JACOB.TXT') READ(5,*) N,NIT,EPS READ(5,*) ((A(I,J),J=1,N),I=1,N) READ(5,*) (X(I),I=1,N) DO 5 I=1,N B(I)=0 5 CONTINUE WRITE(6,*) 'INPUT:' WRITE(6,*) 'N=',N,' NIT=',NIT,' EPS=',EPS WRITE(6,*) 'A=' DO 10 I=1,N WRITE(6,*) (A(I,J),J=1,N) 10 WRITE(6,*) ' ' WRITE(6,*) 'B=',(B(I),I=1,N) WRITE(6,*) 'X=',(X(I),I=1,N) WRITE(6,*) ' ' WRITE(6,*) 'OUTPUT:' CALL JACOBI(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y) WRITE(6,*) 'ITEST=',ITEST,' KIT=',KIT WRITE(6,*) 'X=',(X(I),I=1,N) WRITE(6,*) 'E=',E STOP END C----------------------------------------------------------------SUBROUTINE JACOBI(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y) C C SOLUTION OF THE LINEAR SYSTEM AX=B BY JACOBI'S METHOD C C INPUT: C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS C N: DIMENSION OF ARRAYS C A: SYSTEM MATRIX (NXN) C B: RIGHT-HAND SIDE VECTOR (N) C NIT: MAXIMUM NUMBER OF ITERATIONS C EPS: CONVERGENCE TEST NUMBER C 22 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C C C C C C C OUTPUT: X: APPROXIMATE SOLUTION VECTOR (N) KIT: NUMBER OF EXECUTED ITERATIONS E: ERROR BOUND ITEST = 1, CONVERGENCE TEST SATISFIED = 2, CONVERGENCE TEST NOT SATISFIED DOUBLE PRECISION EPS,A,B,X,Y,CN,S,E DIMENSION A(ND,N),B(ND),X(ND),Y(ND) ITEST=1 C C C ITERATION MATRIX NORM CN=0.D0 DO 10 I=1,N S=0.D0 DO 20 J=1,N IF(J.EQ.I) GO TO 20 S=S+DABS(A(I,J)) 20 CONTINUE S=S/DABS(A(I,I)) IF(S.GT.CN) CN=S 10 CONTINUE C C C ITERATIONS DO 30 K=1,NIT KIT=K DO 60 I=1,N 60 Y(I)=X(I) E=0.D0 DO 40 I=1,N X(I)=B(I) DO 50 J=1,N IF(J.EQ.I) GO TO 50 X(I)=X(I)-A(I,J)*Y(J) 50 CONTINUE X(I)=X(I)/A(I,I) S=DABS(X(I)-Y(I)) IF(S.GT.E) E=S 40 CONTINUE WRITE(6,*) ' X=',(X(I),I=1,N) C C C ERROR TEST E=E*CN/(1.D0-CN) IF(E.LE.EPS) RETURN 23 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 30 CONTINUE ITEST=2 RETURN END Ακολουθεί το πρόγραμμα εκτελεσμένο. INPUT: N= 4 NIT= A= 8.00000000000000 1.00000000000000 40 EPS= 1.000000000000000E-014 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 1.00000000000000 B= 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 X= 10.0000000000000 6.00000000000000 4.00000000000000 2.00000000000000 OUTPUT: X= -1.50000000000000 -2.00000000000000 -2.25000000000000 -2.50000000000000 X= 0.843750000000000 0.781250000000000 0.750000000000000 0.718750000000000 X= -0.281250000000000 -0.289062500000000 -0.292968750000000 -0.296875000000000 X= 0.109863281250000 0.108886718750000 0.108398437500000 0.107910156250000 X= -4.064941406250000E-002 -4.077148437500000E-002 -4.083251953125000E-002 -4.089355468750000E-002 X= 1.531219482421875E-002 1.529693603515625E-002 1.528930664062500E-002 1.528167724609375E-002 X= -5.733489990234375E-003 -5.735397338867188E-003 -5.736351013183594E-003 -5.737304687500000E-003 X= 2.151131629943848E-003 2.150893211364746E-003 2.150774002075195E-003 2.150654792785645E-003 X= -8.065402507781982E-004 -8.065700531005859E-004 -8.065849542617798E-004 -8.065998554229736E-004 X= 3.024693578481674E-004 3.024656325578690E-004 3.024637699127197E-004 24 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 3.024619072675705E-004 X= -1.134239137172699E-004 -1.134243793785572E-004 -1.134246122092009E-004 -1.134248450398445E-004 X= 4.253422957845032E-005 4.253417137078941E-005 4.253414226695895E-005 4.253411316312850E-005 X= -1.595030335010961E-005 -1.595031062606722E-005 -1.595031426404603E-005 -1.595031790202484E-005 X= 5.981367849017261E-006 5.981366939522559E-006 5.981366484775208E-006 5.981366030027857E-006 X= -2.243012431790703E-006 -2.243012545477541E-006 -2.243012602320960E-006 -2.243012659164378E-006 X= 8.411297258703598E-007 8.411297116595051E-007 8.411297045540778E-007 8.411296974486504E-007 X= -3.154236392077792E-007 -3.154236409841360E-007 -3.154236418723144E-007 -3.154236427604928E-007 X= 1.182838657021179E-007 1.182838654800733E-007 1.182838653690510E-007 1.182838652580287E-007 X= -4.435644951339413E-008 -4.435644954114970E-008 -4.435644955502749E-008 -4.435644956890528E-008 X= 1.663366858313531E-008 1.663366857966586E-008 1.663366857793114E-008 1.663366857619641E-008 X= -6.237625716724177E-009 -6.237625717157858E-009 -6.237625717374698E-009 -6.237625717591538E-009 X= 2.339109644015512E-009 2.339109643961302E-009 2.339109643934197E-009 2.339109643907092E-009 X= -8.771661164753237E-010 -8.771661164821000E-010 -8.771661164854881E-010 -8.771661164888762E-010 X= 3.289372936820580E-010 3.289372936812110E-010 3.289372936807875E-010 3.289372936803640E-010 X= -1.233514851302953E-010 -1.233514851304012E-010 -1.233514851304541E-010 -1.233514851305071E-010 X= 4.625680692392030E-011 4.625680692390706E-011 4.625680692390045E-011 4.625680692389383E-011 X= -1.734630259646267E-011 -1.734630259646432E-011 -1.734630259646515E-011 -1.734630259646598E-011 X= 6.504863473674431E-012 6.504863473674224E-012 6.504863473674121E-012 6.504863473674017E-012 X= -2.439323802627795E-012 -2.439323802627821E-012 -2.439323802627834E-012 -2.439323802627847E-012 X= 9.147464259854377E-013 9.147464259854345E-013 9.147464259854329E-013 9.147464259854313E-013 X= -3.430299097445373E-013 -3.430299097445377E-013 -3.430299097445379E-013 -3.430299097445381E-013 X= 1.286362161542017E-013 1.286362161542017E-013 1.286362161542017E-013 1.286362161542016E-013 X= -4.823858105782561E-014 -4.823858105782563E-014 -4.823858105782563E-014 -4.823858105782563E-014 25 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 X= 1.808946789668461E-014 1.808946789668461E-014 1.808946789668461E-014 1.808946789668461E-014 X= -6.783550461256728E-015 -6.783550461256728E-015 -6.783550461256728E-015 -6.783550461256728E-015 X= 2.543831422971273E-015 2.543831422971273E-015 2.543831422971273E-015 2.543831422971273E-015 ITEST= 1 KIT= 36 X= 2.543831422971273E-015 2.543831422971273E-015 2.543831422971273E-015 2.543831422971273E-015 E= 5.596429130536801E-015 Ακολουθεί ο κώδικας της μεθόδου Gauss-Seidel. C C C C PROGRAM GAUSE GAUSS-SEIDEL METHOD PARAMETER(ND=10) DOUBLE PRECISION EPS,A,B,X,Y,E DIMENSION A(ND,ND),B(ND),X(ND),Y(ND) OPEN(UNIT=5,FILE='GAUSE.TXT') READ(5,*) N,NIT,EPS READ(5,*) ((A(I,J),J=1,N),I=1,N) READ(5,*) (X(I),I=1,N) DO 5 I=1,N B(I)=0 5 CONTINUE WRITE(6,*) 'INPUT:' WRITE(6,*) 'N=',N,' NIT=',NIT,' EPS=',EPS WRITE(6,*) 'A=' DO 10 I=1,N WRITE(6,*) (A(I,J),J=1,N) 10 WRITE(6,*) ' ' WRITE(6,*) 'B=',(B(I),I=1,N) WRITE(6,*) 'X=',(X(I),I=1,N) WRITE(6,*) ' ' WRITE(6,*) 'OUTPUT:' CALL GAUSEID(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y) WRITE(6,*) 'ITEST=',ITEST,' KIT=',KIT WRITE(6,*) 'X=',(X(I),I=1,N) WRITE(6,*) 'E=',E STOP END C----------------------------------------------------------------SUBROUTINE GAUSEID(ND,N,A,B,NIT,EPS,X,KIT,E,ITEST,Y) C 26 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 C SOLUTION OF THE LINEAR SYSTEM AX=B BY THE GAUSS-SEIDEL METHOD C C INPUT: C ND: MAXIMUM COLUMN DIMENSION OF ARRAYS C N: DIMENSION OF ARRAYS C A: SYSTEM MATRIX (NXN) C B: RIGHT-HAND SIDE VECTOR (N) C NIT: MAXIMUM NUMBER OF ITERATIONS C EPS: CONVERGENCE TEST NUMBER C C OUTPUT: C X: APPROXIMATE SOLUTION VECTOR (N) C KIT: NUMBER OF EXECUTED ITERATIONS C E: ERROR BOUND C ITEST = 1, CONVERGENCE TEST SATISFIED C = 2, CONVERGENCE TEST NOT SATISFIED C DOUBLE PRECISION EPS,A,B,X,Y,CN,S,E DIMENSION A(ND,N),B(ND),X(ND),Y(ND) ITEST=1 C C NORM OF THE JACOBI ITERATION MATRIX C CN=0.D0 DO 10 I=1,N S=0.D0 DO 20 J=1,N IF(J.EQ.I) GO TO 20 S=S+DABS(A(I,J)) 20 CONTINUE S=S/DABS(A(I,I)) IF(S.GT.CN) CN=S 10 CONTINUE C C ITERATIONS C DO 30 K=1,NIT KIT=K DO 25 I=1,N 25 Y(I)=X(I) E=0.D0 DO 40 I=1,N X(I)=B(I) DO 50 J=1,N IF(J.EQ.I) GO TO 50 X(I)=X(I)-A(I,J)*X(J) 27 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 50 CONTINUE X(I)=X(I)/A(I,I) S=DABS(X(I)-Y(I)) IF(S.GT.E) E=S 40 CONTINUE WRITE(6,*) ' X=',(X(I),I=1,N) C C C ERROR TEST E=E*CN/(1.D0-CN) IF(E.LE.EPS) RETURN 30 CONTINUE ITEST=2 RETURN END Ακολουθεί το πρόγραμμα εκτελεσμένο. INPUT: N= 4 NIT= A= 8.00000000000000 1.00000000000000 40 EPS= 1.000000000000000E-014 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 8.00000000000000 1.00000000000000 1.00000000000000 B= 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 X= 10.0000000000000 6.00000000000000 4.00000000000000 2.00000000000000 OUTPUT: X= -1.50000000000000 -0.562500000000000 7.812500000000000E-003 0.256835937500000 X= 3.723144531250000E-002 -3.773498535156250E-002 -3.204154968261719E-002 4.068136215209961E-003 X= 8.213549852371216E-003 2.469982951879501E-003 -1.843958627432585E-003 -1.104946772102267E-003 X= 5.986530595691875E-005 3.611300116972416E-004 8.549393180601328E-005 -6.331115618252170E-005 28 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 X= -4.791409841509164E-005 3.216415348950008E-006 1.350110490608292E-005 3.899572270007340E-006 X= -2.577136565630033E-006 -1.852942576307528E-006 6.631335899127765E-008 5.454707228682854E-007 X= 1.551448118059956E-007 -9.586611170819484E-008 -7.559367787076078E-008 2.039372221620000E-009 X= 2.117755216966695E-008 6.547094184934229E-009 -3.720502322027647E-009 -3.000518004071691E-009 X= 2.174076764563867E-011 8.374099448067125E-010 2.676709114524175E-010 -1.408527029880961E-010 X= -1.205285191588792E-010 -7.862111631802668E-013 3.277092916376945E-011 1.106797514478626E-011 X= -5.381586643171930E-012 -4.807164708172972E-012 -1.099029741801694E-013 1.287331790690634E-012 X= 4.537169864578135E-013 -2.038932253710348E-013 -1.921444439721766E-013 -7.209914639325260E-015 X= 5.040594799781708E-014 1.861855132671060E-014 -7.726823085650303E-015 -7.662209529859672E-015 X= -4.036898389000781E-016 1.974090306801257E-015 7.614761327448116E-016 -2.914845750807488E-016 X= -3.055102330581650E-016 -2.056016557573724E-017 7.719437171433137E-017 3.110950336494635E-017 ITEST= 1 KIT= 15 X= -3.055102330581650E-016 -2.056016557573724E-017 7.719437171433137E-017 3.110950336494635E-017 E= 1.196790283426196E-015 Παρατηρούμε ότι η x15 ∞ της μεθόδου Jacobi ισούται με 2.243012659164378E-006 ενώ η x15 ∞ της μεθόδου Gauss-Seidel ισούται με 3.055102330581650E-016. Η νόρμα της μεθόδου Gauss-Seidel είναι πάρα πολύ μικρότερη από αυτή της μεθόδου Jacobi( της τάξης των 10 −10 ). Παρατηρούμε ότι η νόρμα της Jacobi φτάνει την τάξη της νόρμας της μεθόδου Gauss-Seidel στην τριακοστή έκτη επανάληψη ( x36 ∞ =2.543831422971273E015 ).Δηλαδή η μέθοδος Gauss-Seidel είναι φανερό ότι συγκλίνει πολύ πιο γρήγορα από ότι η Jacobi (μάλιστα με τη διπλάσια σχεδόν ταχύτητα). ΕΡΩΤΗΜΑ 4: Το πρόγραμμα που ακολουθεί βρίσκει εφαρμογή κατά την αριθμητική επίλυση μερικών διαφορικών εξισώσεων με πεπερασμένα στοιχεία όπου προκύπτουν πίνακες με μορφή ταινίας, δηλαδή πίνακες τέτοιοι ώστε aij=0∀,ijμε−>M. Το πρόγραμμά μας εκμεταλλεύεται τη συγκεκριμένη δομή του πίνακα, δηλαδή οι μετρητές πηγαίνουν ανάλογα με το πλάτος της ταινίας και δεν φτάνουν μέχρι τα μηδενικά. Ένα δεύτερο στοιχείο, χαρακτηριστικό των πινάκων αυτών είναι ότι το μεγαλύτερο στοιχείο καταπόλυτη τιμή είναι το διαγώνιο στοιχείο, οπότε δεν χρειάζεται οδήγηση. Το πρώτο σχόλιο επί του προγράμματος είναι 29 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 ότι οι μετρητές προχωρούν για τριγωνοποίηση μέχρι Μ στοιχεία κάτω από το κάθε στοιχείο της διαγωνίου, αφού όλα τα υπόλοιπα είναι μηδέν. Ένα δεύτερο σχόλιο είναι ότι οι πράξεις στο διάνυσμα Β γίνονται συγχρόνως με την τριγωνοποίηση του πίνακα Α. Ένα τρίτο σχόλιο είναι ότι το Μ εισάγεται στο πρόγραμμα από τον χρήστη. Τέλος, πρέπει να πούμε ότι τα δεδομένα το πρόγραμμα τα διαβάζει από ένα αρχείο και τα αποτελέσματα αποθηκεύονται σε ένα άλλο αρχείο. Ακολουθεί ο κώδικας του προγράμματος. C C C C PROGRAM MATRBAND MATRBAND ELIMINATION METHOD PARAMETER(ND=10) DOUBLE PRECISION A,B,D DIMENSION A(ND,ND),B(ND) OPEN(UNIT=5,FILE='MATR.TXT') OPEN(UNIT=6,FILE='MATR1.TXT') READ(5,*) N,M READ(5,*) ((A(I,J),J=1,N),I=1,N) READ(5,*) (B(I),I=1,N) WRITE(6,*) 'INPUT:' WRITE(6,*) 'N=',N,' M=',M WRITE(6,*) 'A=' DO 10 I=1,N WRITE(6,*) (A(I,J),J=1,N) 10 WRITE(6,*) ' ' WRITE(6,*) 'B=',(B(I),I=1,N) WRITE(6,*) ' ' WRITE(6,*) 'OUTPUT:' CALL MATRBAND(ND,N,N,A,B,D,IREG) WRITE(6,*) 'D=',D,' IREG=',IREG IF(IREG.EQ.0) STOP WRITE(6,*) 'X=',(B(I),I=1,N) STOP END C------------------------------------------------------------------SUBROUTINE MATRBAND(ND,N,M,A,B,D,IREG) DOUBLE PRECISION A,B,D,S,EE,P DIMENSION A(ND,N),B(ND) EE=1.D-12 IREG=1 C C D=1.D0 DO 10 K=1,N-1 D=D*A(K,K) IF (DABS(A(K,K)).LE.EE) GO TO 100 DO 50 I=K+1,K+M 30 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 IF(I.LE.N)THEN P=A(I,K)/A(K,K) DO 40 J=K+1,K+M IF(J.LE.N) THEN A(I,J)=A(I,J)-P*A(K,J) END IF 40 CONTINUE B(I)=B(I)-P*B(K) END IF 50 CONTINUE 10 CONTINUE D=D*A(N,N) IF(DABS(A(N,N)).LE.EE) GO TO 100 C C B(N)=B(N)/A(N,N) DO 80 I=1,N-1 K=N-I S=0.D0 DO 90 J=K+1,K+M IF(J.LE.N) THEN S=S+A(K,J)*B(J) END IF 90 CONTINUE B(K)=(B(K)-S)/A(K,K) 80 CONTINUE RETURN 100 IREG=0 RETURN END Ακολουθεί το αρχείο στο οποίο αποθηκεύτηκαν τα αποτελέσματα του προγράμματος. INPUT: N= 9 M= 3 A= 4.00000000000000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 31 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 -1.00000000000000 4.00000000000000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 4.00000000000000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 -1.00000000000000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 0.000000000000000E+000 -1.00000000000000 0.000000000000000E+000 -1.00000000000000 4.00000000000000 B= 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 1.000000000000000E-002 OUTPUT: D= 100352.000000000 IREG= 1 X= 6.875000000000001E-003 8.750000000000003E-003 6.875000000000001E-003 8.750000000000001E-003 1.125000000000000E-002 8.750000000000001E-003 6.875000000000000E-003 8.750000000000001E-003 6.875000000000001E-003. Για επαλήθευση των αποτελεσμάτων αυτών λύσαμε το σύστημα με το Mathematica και πήραμε τα ακόλουθα αποτελέσματα: @ 8 NSolve 4 x1 - x2 - x4 0.01, - x1 + 4 x2 - x3 - x5 0.01, - x2 + 4 x3 - x6 0.01, - x1 + 4 x4 - x5 - x7 0.01, - x2 - x4 + 4 x5 - x6 - x8 0.01, - x3 - x5 + 4 x6 - x9 0.01, - x4 + 4 x7 - x8 0.01, - x5 - x7 + 4 x8 - x9 0.01, - x6 - x8 + 4 x9 0.01 , x1, x2, x3, x4, x5, x6, x7, x8, x9 < 8 < D 32 8 8 ΑΡΙΘΜΗΤΙΚΗ ΑΝΑΛΥΣΗ-ΕΡΓΑΣΤΗΡΙΑΚΗ ΑΣΚΗΣΗ 2 < < x1 ® 0.006875, x2 ® 0.00875, x3 ® 0.006875, x4 ® 0.00875, x5 ® 0.01125, x6 ® 0.00875, x7 ® 0.006875, x8 ® 0.00875, x9 ® 0.006875 Και όπως παρατηρούμε είμαστε σωστοί. . 33 ...
View Full Document

Ask a homework question - tutors are online