1 : C XKM SOURCE CHAT 05/01/13 04:13:42 5004 2 : SUBROUTINE XKM(NCOURB,TRAC,XK,XM,W,F) 3 : C 4 : C =============================================================== 5 : C CE SOUS-PROGRAMME EST APPELE PAR ENDOM1. 6 : C IL CALCULE LES COEFF. K ET M DE LA FONCTION DECRIVANT LA COURBE 7 : C DE TRACTION: Y=K*X**(1/M) OU X EST UNE DEFORM. ET Y UNE CONTR. 8 : C 9 : C ENTREES: 10 : C ------- 11 : C NCOURB = NBR. DE PTS. DE LA COURBE DE TRACTION 12 : C TRAC(2*NCOURB)= COURBE DE TRACTION 13 : C 14 : C SORTIES: 15 : C ------- 16 : C XK= COEFF. K 17 : C XM= COEFF. M 18 : C ================================================================ 19 : IMPLICIT INTEGER(I-N) 20 : IMPLICIT REAL*8 (A-H,O-Z) 21 : DIMENSION TRAC(*),B(2),C(2),W(*) 22 : DIMENSION A(2,2),F(NCOURB,2) 23 : C 24 : DATA TOL,MAXITE/1.D-3,15/ 25 : C 26 : C 27 : C ========================================== 28 : C CAS OU LA COURBE DE TRACTION N'A QUE 2 PTS 29 : C ========================================== 30 : IF (NCOURB.EQ.2) THEN 31 : XK=(TRAC(3)-TRAC(1))/(TRAC(4)-TRAC(2)) 32 : XM=1.D0 33 : GOTO 100 34 : ENDIF 35 : 36 : C ======================================= 37 : C METHODE DES MOINDRES CARRES PONDERES 38 : C ON LINEARISE: LOG(Y)=LOG(K)+1/M*LOG(X) 39 : C C'EST DE LA FORME: YY=C1*F1(X)+C2*F2(X) 40 : C ON RESOUD: A(2,2)C(2)=B(2) 41 : C ======================================= 42 : SIGY=TRAC(1) 43 : DO 10 I=1,NCOURB 44 : W(I)=(TRAC(2*I-1)-SIGY)*(TRAC(2*I-1)-SIGY) 45 : 10 CONTINUE 46 : C 46 : CALL ZERO(A,2,2) 46 : CALL ZDANUL(B,2) 47 : C 48 : C ............................................................... 49 : C DS. LA SUITE, LE 1ER PT. DE LA COURBE EST OMIS CAR TRAC(1)=SIGY 50 : C ET TRAC(2)=0 CE QUI POSE PB. POUR LE CALCUL DES LOGARITHMES 51 : C ............................................................... 52 : DO 30 I=2,NCOURB 53 : A(1,1)=A(1,1)+W(I) 54 : A(1,2)=A(1,2)+W(I)*DLOG(TRAC(2*I)) 55 : A(2,2)=A(2,2)+W(I)*DLOG(TRAC(2*I))*DLOG(TRAC(2*I)) 56 : B(1)=B(1)+W(I)*DLOG(TRAC(2*I-1)-SIGY) 57 : B(2)=B(2)+W(I)*DLOG(TRAC(2*I-1)-SIGY)*DLOG(TRAC(2*I)) 58 : 30 CONTINUE 59 : DET=A(1,1)*A(2,2)-A(1,2)*A(1,2) 60 : C(1)=(B(1)*A(2,2)-B(2)*A(1,2))/DET 61 : C(2)=(B(2)*A(1,1)-B(1)*A(1,2))/DET 62 : XK=DEXP(C(1)) 63 : XM=1.D0/C(2) 64 : C 65 : C ---------------------------------------------------------------- 66 : C SI LA METHODE DE GAUSS-NEWTON NE CONVERGE PAS, ON RETIENT, COMME 67 : C SOLUTIONS DU PB., XK ET XM CALCULES PAR LA METHODE PRECEDENTE, 68 : C QU'ON STOCKE DS. XK0 ET XM0: 69 : C ---------------------------------------------------------------- 70 : XK0=XK 71 : XM0=XM 72 : C 73 : C ======================= 74 : C METHODE DE GAUSS-NEWTON 75 : C ======================= 76 : NITER=0 77 : 40 DO 50 I=2,NCOURB 78 : F(I,1)=TRAC(2*I)**(1.D0/XM) 79 : F(I,2)=-(XK/XM**2.D0)*DLOG(TRAC(2*I))*(TRAC(2*I)**(1.D0/XM)) 80 : 50 CONTINUE 81 : C 82 : DO 60 I=1,2 83 : DO 60 J=1,2 84 : A(I,J)=0.D0 85 : DO 60 K=2,NCOURB 86 : A(I,J)=A(I,J)+F(K,I)*F(K,J) 87 : 60 CONTINUE 88 : C 89 : DO 70 I=1,2 90 : B(I)=0.D0 91 : DO 70 K=2,NCOURB 92 : B(I)=B(I)+F(K,I)*(TRAC(2*K-1)-SIGY-XK*TRAC(2*K)**(1.D0/XM)) 93 : 70 CONTINUE 94 : C 95 : DET=A(1,1)*A(2,2)-A(1,2)*A(2,1) 96 : DXK=(B(1)*A(2,2)-B(2)*A(1,2))/DET 97 : DXM=(B(2)*A(1,1)-B(1)*A(2,1))/DET 98 : C 99 : C ------------------------------------------------------------- 100 : C CALCUL DE L'ERREUR AV. (Z1) ET APRES (Z2) LA MAJ. DE XK ET XM 101 : C ------------------------------------------------------------- 102 : Z1=0.D0 103 : DO 80 K=2,NCOURB 104 : ZZ1=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM)) 105 : ZZ1=ZZ1*ZZ1 106 : Z1=Z1+ZZ1 107 : 80 CONTINUE 108 : C 109 : XK=XK+DXK 110 : XM=XM+DXM 111 : C 112 : Z2=0.D0 113 : DO 90 K=2,NCOURB 114 : ZZ2=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM)) 115 : ZZ2=ZZ2*ZZ2 116 : Z2=Z2+ZZ2 117 : 90 CONTINUE 118 : C 119 : C ---------------------------- 120 : C TESTS D'ARRET DES ITERATIONS 121 : C ---------------------------- 122 : NITER=NITER+1 123 : C ARRET=DABS(Z2-Z1)/Z2 124 : C IF (ARRET.LT.TOL) GOTO 100 125 : ARRET=DABS(Z2-Z1) 126 : IF (ARRET.LT.(TOL*Z2)) GOTO 100 127 : IF (NITER.LT.MAXITE) GOTO 40 128 : C 129 : C --------------------------------------------------------------- 130 : C SI ON PASSE ICI, C'EST QUE G-NEWTON N'A PAS CONVERGE AU BOUT DE 131 : C MAXITE ITER.: LA SOL. EST CELLE CALCULEE PAR LA METHODE DES 132 : C MOINDRES CARRES PONDERES 133 : C --------------------------------------------------------------- 134 : XK=XK0 135 : XM=XM0 136 : 100 CONTINUE 137 : RETURN 138 : END 139 : 140 : 141 :
© Cast3M 2003 - All rights reserved.
Disclaimer