亿万人中,你我独一无二分享 http://blog.sciencenet.cn/u/stone1913

博文

FORTRAN programs of gradients method

已有 3020 次阅读 2016-1-15 17:07 |个人分类:经验交流|系统分类:科研笔记| method, gradients

Thefollowing are the FORTRAN programs (FORTRAN 90) thatwere used to compute the real root of a nonlinear equations system, byusing the gradients method.

 

CALL GRADIENTS (X01,X02, YT1, YT2, XP, YP, ZP, X, Y, Z, CT)

 

SUBROUTINE GRADIENTS (X01,X02, YT1, YT2, XP, YP, ZP, X, Y, Z, CT)

C   X01, X02 are the given initial values

C   YT1, YT2 are the final solutions

C   XP, YP, ZP are the coordinates of thecalculation points

C   X, Y, Z are the coordinates of the eightnodal points of the a curved quadrilateral element

C   CT is the anisotropic tensor

C   EPS is a stopping criterion

IMPLICIT NONE

DOUBLE PRECISION  X(8), Y(8), Z(8), X01, X02, YT1,YT2, XP, YP, ZP, EPS, Y1, Y2

DOUBLE PRECISION F, F1, F2, CT(3, 3)

INTEGER L

 

L=5000

   EPS=1.0E-35

   YT1=X01

   YT2=X02        

 

   CALL DSNSE(EPS, YT1,YT2, Y1, Y2, L, XP, YP, ZP, X, Y, Z, CT)

   IF(L.GT.0) THEN

   WRITE(11,*) 'iteration number = ',5000-L  

   ENDIF

   RETURN

   END

 

 

SUBROUTINE DSNSE(EPS,YT1, YT2, Y1, Y2, L, XP, YP, ZP, X, Y, Z, CT)

   IMPLICIT NONE

      DOUBLE PRECISION EPS, YT1, YT2, Y1, Y2, X(8), Y(8), Z(8), XP, YP, ZP, F, D, S, CT(3, 3)

   INTEGER L

 

5   CALL FS(YT1, YT2, F, Y1, Y2, XP, YP, ZP, X, Y, Z, CT)

   IF(DSQRT(F).GE.EPS) THEN

   L=L-1

   IF(L.EQ.0) RETURN

    D=Y1**2+Y2**2

     IF(D.EQ.0)THEN

    L=-1

     RETURN

     ENDIF

    S=F/D

    YT1=YT1-S*Y1

    YT2=YT2-S*Y2

     GOTO 5

     ENDIF

     RETURN

     END

 

SUBROUTINE FS(YT1,YT2, F, Y1, Y2, XP, YP, ZP, X, Y, Z, CT)

   IMPLICIT NONE

      DOUBLE PRECISION YT1, YT2, Y1, Y2, X(8), Y(8), Z(8), XP, YP, ZP, F, XK, YK, ZK

   DOUBLE PRECISION F1, F2, DF1, DF2, DX1, DY1, DZ1, DX2, DY2, DZ2, CT(3, 3)

   DOUBLE PRECISION Q1, Q2, Q3, Q4, Q5, Q6, Q7, Q8, DF11, DF12, DF21, DF22, DX11, DY11, DZ11

   DOUBLE PRECISION DX12, DY12, DZ12, RR1, RR2, RR3

   DOUBLE PRECISION DX21, DY21, DZ21

   DOUBLE PRECISION DX22, DY22, DZ22

 

   DX1 =(1-YT2)*(2.D0*YT1+YT2)/4.D0*X(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*X(2)

    +           + (1+YT2)*(2.D0*YT1+YT2)/4.D0*X(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*X(4)

    +           + YT1*(YT2-1.D0)*X(5)

    +           + (1.D0-YT2**2)/2.D0*X(6)

    +           -  YT1*(1.D0+YT2)*X(7)

    +           +(YT2**2-1.D0)/2.D0*X(8)  

 

      DY1 = (1-YT2)*(2.D0*YT1+YT2)/4.D0*Y(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*Y(2)

    +           +(1+YT2)*(2.D0*YT1+YT2)/4.D0*Y(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*Y(4)

    +           + YT1*(YT2-1.D0)*Y(5)

    +           + (1.D0-YT2**2)/2.D0*Y(6)

    +           -  YT1*(1.D0+YT2)*Y(7)

    +           +(YT2**2-1.D0)/2.D0*Y(8)  

 

       DZ1 = (1-YT2)*(2.D0*YT1+YT2)/4.D0*Z(1)

    +           +(1-YT2)*(2.D0*YT1-YT2)/4.D0*Z(2)

    +           +(1+YT2)*(2.D0*YT1+YT2)/4.D0*Z(3)

    +           +(1+YT2)*(2.D0*YT1-YT2)/4.D0*Z(4)

    +           + YT1*(YT2-1.D0)*Z(5)

    +           + (1.D0-YT2**2)/2.D0*Z(6)

    +           -  YT1*(1.D0+YT2)*Z(7)

    +           + (YT2**2-1.D0)/2.D0*Z(8)

 

       DX2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*X(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*X(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*X(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*X(4)

    +           + (YT1**2-1.D0)/2.D0*X(5)

    +           - YT2*(1.D0+YT1)*X(6)

    +           + (1.D0-YT1**2)/2.D0*X(7)

    +           +YT2*(YT1-1.D0)*X(8)  

 

       DY2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*Y(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*Y(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*Y(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*Y(4)

    +           + (YT1**2-1.D0)/2.D0*Y(5)

    +           - YT2*(1.D0+YT1)*Y(6)

    +           + (1.D0-YT1**2)/2.D0*Y(7)

    +           + YT2*(YT1-1.D0)*Y(8)  

 

       DZ2 = (1-YT1)*(2.D0*YT2+YT1)/4.D0*Z(1)

    +           -(1+YT1)*(-2.D0*YT2+YT1)/4.D0*Z(2)

    +           +(1+YT1)*(2.D0*YT2+YT1)/4.D0*Z(3)

    +           +(-1+YT1)*(-2.D0*YT2+YT1)/4.D0*Z(4)

    +           + (YT1**2-1.D0)/2.D0*Z(5)

    +           - YT2*(1.D0+YT1)*Z(6)

    +           + (1.D0-YT1**2)/2.D0*Z(7)

    +           + YT2*(YT1-1.D0)*Z(8)

       

       Q1 = (1-YT1)*(1-YT2)*(-YT1-YT2-1.D0)/4.D0

       Q2 = (1+YT1)*(1-YT2)*(YT1-YT2-1.D0)/4.D0

       Q3 = (1+YT1)*(1+YT2)*(YT1+YT2-1.D0)/4.D0

       Q4 = (1-YT1)*(1+YT2)*(-YT1+YT2-1.D0)/4.D0

       Q5 = (1-YT1**2)*(1-YT2)/2.D0

       Q6 = (1+YT1)*(1-YT2**2)/2.D0

       Q7 = (1+YT2)*(1-YT1**2)/2.D0

       Q8 = (1-YT1)*(1-YT2**2)/2.D0

 

       XK=Q1*X(1)+Q2*X(2)+Q3*X(3)+Q4*X(4)+Q5*X(5)+Q6*X(6)+Q7*X(7)+Q8*X(8)

       YK=Q1*Y(1)+Q2*Y(2)+Q3*Y(3)+Q4*Y(4)+Q5*Y(5)+Q6*Y(6)+Q7*Y(7)+Q8*Y(8)

       ZK=Q1*Z(1)+Q2*Z(2)+Q3*Z(3)+Q4*Z(4)+Q5*Z(5)+Q6*Z(6)+Q7*Z(7)+Q8*Z(8)

 

        RR1=XK-XP

        RR2=YK-YP

        RR3=ZK-ZP

        F1=CT(1,1)*DX1*RR1+CT(1,2)*DX1*RR2+CT(1,3)*DX1*RR3

    +      +CT(2,1)*DY1*RR1+CT(2,2)*DY1*RR2+CT(2,3)*DY1*RR3

    +      +CT(3,1)*DZ1*RR1+CT(3,2)*DZ1*RR2+CT(3,3)*DZ1*RR3

 

       F2=CT(1,1)*DX2*RR1+CT(1,2)*DX2*RR2+CT(1,3)*DX2*RR3

    +     +CT(2,1)*DY2*RR1+CT(2,2)*DY2*RR2+CT(2,3)*DY2*RR3

    +     +CT(3,1)*DZ2*RR1+CT(3,2)*DZ2*RR2+CT(3,3)*DZ2*RR3

       

       F=F1*F1+F2*F2

 

       DX11=(1-YT2)/2.D0*X(1)+(1-YT2)/2.D0*X(2)+(1+YT2)/2.D0*X(3)+(1+YT2)/2.D0*X(4)+(YT2-1)*X(5)+(-YT2-1)*X(7)

       DY11=(1-YT2)/2.D0*Y(1)+(1-YT2)/2.D0*Y(2)+(1+YT2)/2.D0*Y(3)+(1+YT2)/2.D0*Y(4)+(YT2-1)*Y(5)+(-YT2-1)*Y(7)

       DZ11=(1-YT2)/2.D0*Z(1)+(1-YT2)/2.D0*Z(2)+(1+YT2)/2.D0*Z(3)+(1+YT2)/2.D0*Z(4)+(YT2-1)*Z(5)+(-YT2-1)*Z(7)

 

       DF11=CT(1,1)*DX11*RR1+CT(1,1)*DX1*DX1+CT(1,2)*DX11*RR2+CT(1,2)*DX1*DY1+CT(1,3)*DX11*RR3+CT(1,3)*DX1*DZ1

    +           +CT(2,1)*DY11*RR1+CT(2,1)*DY1*DX1+CT(2,2)*DY11*RR2+CT(2,2)*DY1*DY1+CT(2,3)*DY11*RR3+CT(2,3)*DY1*DZ1

    +          +CT(3,1)*DZ11*RR1+CT(3,1)*DZ1*DX1+CT(3,2)*DZ11*RR2+CT(3,2)*DZ1*DY1+CT(3,3)*DZ11*RR3+CT(3,3)*DZ1*DZ1        

 

       DX12=(1-2.D0*YT1-2.D0*YT2)/4.D0*X(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*X(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*X(3)

    +          +(-1+2.D0*YT1-2.D0*YT2)/4.D0*X(4)

    +           +YT1*X(5)-YT2*X(6)-YT1*X(7)+YT2*X(8)

       DY12=(1-2.D0*YT1-2.D0*YT2)/4.D0*Y(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*Y(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*Y(3)

    +           +(-1+2.D0*YT1-2.D0*YT2)/4.D0*Y(4)

    +           +YT1*Y(5)-YT2*Y(6)-YT1*Y(7)+YT2*Y(8)

       DZ12=(1-2.D0*YT1-2.D0*YT2)/4.D0*Z(1)+(-1-2.D0*YT1+2.D0*YT2)/4.D0*Z(2)+(1+2.D0*YT1+2.D0*YT2)/4.D0*Z(3)

    +          +(-1+2.D0*YT1-2.D0*YT2)/4.D0*Z(4)

    +           +YT1*Z(5)-YT2*Z(6)-YT1*Z(7)+YT2*Z(8)

 

       DF12=CT(1,1)*DX12*RR1+CT(1,1)*DX1*DX2+CT(1,2)*DX12*RR2+CT(1,2)*DX1*DY2+CT(1,3)*DX12*RR3+CT(1,3)*DX1*DZ2    !!!

    +          +CT(2,1)*DY12*RR1+CT(2,1)*DY1*DX2+CT(2,2)*DY12*RR2+CT(2,2)*DY1*DY2+CT(2,3)*DY12*RR3+CT(2,3)*DY1*DZ2

    +          +CT(3,1)*DZ12*RR1+CT(3,1)*DZ1*DX2+CT(3,2)*DZ12*RR2+CT(3,2)*DZ1*DY2+CT(3,3)*DZ12*RR3+CT(3,3)*DZ1*DZ2    

       

       DX21=DX12

       DY21=DY12

       DZ21=DZ12        

       

       DF21=CT(1,1)*DX21*RR1+CT(1,1)*DX2*DX1+CT(1,2)*DX21*RR2+CT(1,2)*DX2*DY1+CT(1,3)*DX21*RR3+CT(1,3)*DX2*DZ1    !!!

    +           +CT(2,1)*DY21*RR1+CT(2,1)*DY2*DX1+CT(2,2)*DY21*RR2+CT(2,2)*DY2*DY1+CT(2,3)*DY21*RR3+CT(2,3)*DY2*DZ1

    +          +CT(3,1)*DZ21*RR1+CT(3,1)*DZ2*DX1+CT(3,2)*DZ21*RR2+CT(3,2)*DZ2*DY1+CT(3,3)*DZ21*RR3+CT(3,3)*DZ2*DZ1        

 

       DX22=(1-YT1)/2.D0*X(1)+(1+YT1)/2.D0*X(2)+(1+YT1)/2.D0*X(3)+(1-YT1)/2.D0*X(4)+(-YT1-1)*X(6)+(YT1-1)*X(8)

       DY22=(1-YT1)/2.D0*Y(1)+(1+YT1)/2.D0*Y(2)+(1+YT1)/2.D0*Y(3)+(1-YT1)/2.D0*Y(4)+(-YT1-1)*Y(6)+(YT1-1)*Y(8)

       DZ22=(1-YT1)/2.D0*Z(1)+(1+YT1)/2.D0*Z(2)+(1+YT1)/2.D0*Z(3)+(1-YT1)/2.D0*Z(4)+(-YT1-1)*Z(6)+(YT1-1)*Z(8)

 

       DF22=CT(1,1)*DX22*RR1+CT(1,1)*DX2*DX2+CT(1,2)*DX22*RR2+CT(1,2)*DX2*DY2+CT(1,3)*DX22*RR3+CT(1,3)*DX2*DZ2    !!!

    +          +CT(2,1)*DY22*RR1+CT(2,1)*DY2*DX2+CT(2,2)*DY22*RR2+CT(2,2)*DY2*DY2+CT(2,3)*DY22*RR3+CT(2,3)*DY2*DZ2

    +          +CT(3,1)*DZ22*RR1+CT(3,1)*DZ2*DX2+CT(3,2)*DZ22*RR2+CT(3,2)*DZ2*DY2+CT(3,3)*DZ22*RR3+CT(3,3)*DZ2*DZ2

 

       Y1=2.D0*(F1*DF11+F2*DF21)

       Y2=2.D0*(F1*DF12+F2*DF22)

     RETURN

     END

 

 




https://blog.sciencenet.cn/blog-425262-950615.html

上一篇:64位Intel Fortran 调用64位Matlab
收藏 IP: 123.235.86.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 00:24

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部