Télécharger la plaquette Cast3M


fr en

Télécharger xkm.eso

Retour à la liste

Numérotation des lignes : 

xkm

  1. C XKM       SOURCE    CB215821  16/04/21    21:18:42     8920
  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
  47.       CALL ZERO(A,2,2)
  48.       CALL ZDANUL(B,2)
  49. C
  50. C ...............................................................
  51. C DS. LA SUITE, LE 1ER PT. DE LA COURBE EST OMIS CAR TRAC(1)=SIGY
  52. C ET TRAC(2)=0 CE QUI POSE PB. POUR LE CALCUL DES LOGARITHMES
  53. C ...............................................................
  54.       DO 30 I=2,NCOURB
  55.            A(1,1)=A(1,1)+W(I)
  56.            A(1,2)=A(1,2)+W(I)*LOG(TRAC(2*I))
  57.            A(2,2)=A(2,2)+W(I)*LOG(TRAC(2*I))*LOG(TRAC(2*I))
  58.            B(1)=B(1)+W(I)*LOG(TRAC(2*I-1)-SIGY)
  59.            B(2)=B(2)+W(I)*LOG(TRAC(2*I-1)-SIGY)*LOG(TRAC(2*I))
  60.    30 CONTINUE
  61.       DET=A(1,1)*A(2,2)-A(1,2)*A(1,2)
  62.       C(1)=(B(1)*A(2,2)-B(2)*A(1,2))/DET
  63.       C(2)=(B(2)*A(1,1)-B(1)*A(1,2))/DET
  64.       XK=EXP(C(1))
  65.       XM=1.D0/C(2)
  66. C
  67. C ----------------------------------------------------------------
  68. C SI LA METHODE DE GAUSS-NEWTON NE CONVERGE PAS, ON RETIENT, COMME
  69. C SOLUTIONS DU PB., XK ET XM CALCULES PAR LA METHODE PRECEDENTE,
  70. C QU'ON STOCKE DS. XK0 ET XM0:
  71. C ----------------------------------------------------------------
  72.       XK0=XK
  73.       XM0=XM
  74. C
  75. C =======================
  76. C METHODE DE GAUSS-NEWTON
  77. C =======================
  78.       NITER=0
  79.    40 DO 50 I=2,NCOURB
  80.       F(I,1)=TRAC(2*I)**(1.D0/XM)
  81.       F(I,2)=-(XK/XM**2.D0)*LOG(TRAC(2*I))*(TRAC(2*I)**(1.D0/XM))
  82.    50 CONTINUE
  83. C
  84.       DO 60 I=1,2
  85.        DO 60 J=1,2
  86.         A(I,J)=0.D0
  87.         DO 60 K=2,NCOURB
  88.          A(I,J)=A(I,J)+F(K,I)*F(K,J)
  89.    60 CONTINUE
  90. C
  91.       DO 70 I=1,2
  92.         B(I)=0.D0
  93.         DO 70 K=2,NCOURB
  94.           B(I)=B(I)+F(K,I)*(TRAC(2*K-1)-SIGY-XK*TRAC(2*K)**(1.D0/XM))
  95.    70   CONTINUE
  96. C
  97.       DET=A(1,1)*A(2,2)-A(1,2)*A(2,1)
  98.       DXK=(B(1)*A(2,2)-B(2)*A(1,2))/DET
  99.       DXM=(B(2)*A(1,1)-B(1)*A(2,1))/DET
  100. C
  101. C -------------------------------------------------------------
  102. C CALCUL DE L'ERREUR AV. (Z1) ET APRES (Z2) LA MAJ. DE XK ET XM
  103. C -------------------------------------------------------------
  104.       Z1=0.D0
  105.       DO 80 K=2,NCOURB
  106.         ZZ1=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM))
  107.         ZZ1=ZZ1*ZZ1
  108.         Z1=Z1+ZZ1
  109.    80 CONTINUE
  110. C
  111.       XK=XK+DXK
  112.       XM=XM+DXM
  113. C
  114.       Z2=0.D0
  115.       DO 90 K=2,NCOURB
  116.         ZZ2=TRAC(2*K-1)-SIGY-XK*(TRAC(2*K)**(1.D0/XM))
  117.         ZZ2=ZZ2*ZZ2
  118.         Z2=Z2+ZZ2
  119.    90 CONTINUE
  120. C
  121. C ----------------------------
  122. C TESTS D'ARRET DES ITERATIONS
  123. C ----------------------------
  124.       NITER=NITER+1
  125. C      ARRET=ABS(Z2-Z1)/Z2
  126. C      IF (ARRET.LT.TOL) GOTO 100
  127.       ARRET=ABS(Z2-Z1)
  128.       IF (ARRET.LT.(TOL*Z2)) GOTO 100
  129.       IF (NITER.LT.MAXITE) GOTO 40
  130. C
  131. C ---------------------------------------------------------------
  132. C SI ON PASSE ICI, C'EST QUE G-NEWTON N'A PAS CONVERGE AU BOUT DE
  133. C MAXITE ITER.: LA SOL. EST CELLE CALCULEE PAR LA METHODE DES
  134. C MOINDRES CARRES PONDERES
  135. C ---------------------------------------------------------------
  136.       XK=XK0
  137.       XM=XM0
  138.   100 CONTINUE
  139.       RETURN
  140.       END
  141.  
  142.  
  143.  
  144.  
  145.  
  146.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales