C XKM SOURCE CHAT 05/01/13 04:13:42 5004 SUBROUTINE XKM(NCOURB,TRAC,XK,XM,W,F) C C =============================================================== C CE SOUS-PROGRAMME EST APPELE PAR ENDOM1. C IL CALCULE LES COEFF. K ET M DE LA FONCTION DECRIVANT LA COURBE C DE TRACTION: Y=K*X**(1/M) OU X EST UNE DEFORM. ET Y UNE CONTR. C C ENTREES: C ------- C NCOURB = NBR. DE PTS. DE LA COURBE DE TRACTION C TRAC(2*NCOURB)= COURBE DE TRACTION C C SORTIES: C ------- C XK= COEFF. K C XM= COEFF. M C ================================================================ IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION TRAC(*),B(2),C(2),W(*) DIMENSION A(2,2),F(NCOURB,2) C DATA TOL,MAXITE/1.D-3,15/ C C C ========================================== C CAS OU LA COURBE DE TRACTION N'A QUE 2 PTS C ========================================== IF (NCOURB.EQ.2) THEN XK=(TRAC(3)-TRAC(1))/(TRAC(4)-TRAC(2)) XM=1.D0 GOTO 100 ENDIF C ======================================= C METHODE DES MOINDRES CARRES PONDERES C ON LINEARISE: LOG(Y)=LOG(K)+1/M*LOG(X) C C'EST DE LA FORME: YY=C1*F1(X)+C2*F2(X) C ON RESOUD: A(2,2)C(2)=B(2) C ======================================= SIGY=TRAC(1) DO 10 I=1,NCOURB W(I)=(TRAC(2*I-1)-SIGY)*(TRAC(2*I-1)-SIGY) 10 CONTINUE C CALL ZERO(A,2,2) CALL ZDANUL(B,2) C C ............................................................... C DS. LA SUITE, LE 1ER PT. DE LA COURBE EST OMIS CAR TRAC(1)=SIGY C ET TRAC(2)=0 CE QUI POSE PB. POUR LE CALCUL DES LOGARITHMES C ............................................................... DO 30 I=2,NCOURB A(1,1)=A(1,1)+W(I) A(1,2)=A(1,2)+W(I)*DLOG(TRAC(2*I)) A(2,2)=A(2,2)+W(I)*DLOG(TRAC(2*I))*DLOG(TRAC(2*I)) B(1)=B(1)+W(I)*DLOG(TRAC(2*I-1)-SIGY) B(2)=B(2)+W(I)*DLOG(TRAC(2*I-1)-SIGY)*DLOG(TRAC(2*I)) 30 CONTINUE DET=A(1,1)*A(2,2)-A(1,2)*A(1,2) C(1)=(B(1)*A(2,2)-B(2)*A(1,2))/DET C(2)=(B(2)*A(1,1)-B(1)*A(1,2))/DET XK=DEXP(C(1)) XM=1.D0/C(2) C C ---------------------------------------------------------------- C SI LA METHODE DE GAUSS-NEWTON NE CONVERGE PAS, ON RETIENT, COMME C SOLUTIONS DU PB., XK ET XM CALCULES PAR LA METHODE PRECEDENTE, C QU'ON STOCKE DS. XK0 ET XM0: C ---------------------------------------------------------------- XK0=XK XM0=XM C C ======================= C METHODE DE GAUSS-NEWTON C ======================= NITER=0 40 DO 50 I=2,NCOURB F(I,1)=TRAC(2*I)**(1.D0/XM) F(I,2)=-(XK/XM**2.D0)*DLOG(TRAC(2*I))*(TRAC(2*I)**(1.D0/XM)) 50 CONTINUE C DO 60 I=1,2 DO 60 J=1,2 A(I,J)=0.D0 DO 60 K=2,NCOURB A(I,J)=A(I,J)+F(K,I)*F(K,J) 60 CONTINUE C DO 70 I=1,2 B(I)=0.D0 DO 70 K=2,NCOURB B(I)=B(I)+F(K,I)*(TRAC(2*K-1)-SIGY-XK*TRAC(2*K)**(1.D0/XM)) 70 CONTINUE C DET=A(1,1)*A(2,2)-A(1,2)*A(2,1) DXK=(B(1)*A(2,2)-B(2)*A(1,2))/DET DXM=(B(2)*A(1,1)-B(1)*A(2,1))/DET C C ------------------------------------------------------------- C CALCUL DE L'ERREUR AV. (Z1) ET APRES (Z2) LA MAJ. DE XK ET XM C ------------------------------------------------------------- Z1=0.D0 DO 80 K=2,NCOURB ZZ1=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM)) ZZ1=ZZ1*ZZ1 Z1=Z1+ZZ1 80 CONTINUE C XK=XK+DXK XM=XM+DXM C Z2=0.D0 DO 90 K=2,NCOURB ZZ2=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM)) ZZ2=ZZ2*ZZ2 Z2=Z2+ZZ2 90 CONTINUE C C ---------------------------- C TESTS D'ARRET DES ITERATIONS C ---------------------------- NITER=NITER+1 C ARRET=DABS(Z2-Z1)/Z2 C IF (ARRET.LT.TOL) GOTO 100 ARRET=DABS(Z2-Z1) IF (ARRET.LT.(TOL*Z2)) GOTO 100 IF (NITER.LT.MAXITE) GOTO 40 C C --------------------------------------------------------------- C SI ON PASSE ICI, C'EST QUE G-NEWTON N'A PAS CONVERGE AU BOUT DE C MAXITE ITER.: LA SOL. EST CELLE CALCULEE PAR LA METHODE DES C MOINDRES CARRES PONDERES C --------------------------------------------------------------- XK=XK0 XM=XM0 100 CONTINUE RETURN END