| 查看: 1406 | 回复: 1 | ||
| 【悬赏金币】回答本帖问题，作者高湘伦将赠送您 30 个金币 | ||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | ||
[求助]
高分求黄平《润滑数值计算方法》中多重网格积分法原理程序解读的注释
|
||
|
DIMENSION X(1100),P(1100),H(1100),W(2200) COMMON /COMAK/AK(0:1100) DATA NW,PAI/2200,3.14159265/ OPEN(8,FILE='OUT.DAT',STATUS='UNKNOWN') N=129 CALL SUBAK(N) DX=3./(N-1) DO I=1,N X(I)=-1.5+(I-1)*DX P(I)=0.0 IF(ABS(X(I)).LE.1.0)P(I)=SQRT(1.-X(I)*X(I)) ENDDO K=3 CALL DISP(N,NW,K,DX,P,W) WX=0.5*PAI*DX*ALOG(DX) DO I=1,N W(I)=WX DO J=1,N IJ=IABS(I-J) W(I)=W(I)+AK(IJ)*P(J)*DX ENDDO ENDDO DO 30 I=1,N H(I)=1.24+0.5*X(I)*X(I)-W(I)/PAI 30 CONTINUE DO I=1,N WRITE(8,40)X(I),P(I),H(I) ENDDO 40 FORMAT(1X,6(E12.6,1X)) STOP END SUBROUTINE DISP(N,NW,KMAX,DX,P1,W) DIMENSION P1(N),W(NW),P(2200),AK1(0:50),AK2(0:50) COMMON /COMAK/AK(0:1100) DATA NMAX,KMIN/2200,1/ N2=N M=3+2*ALOG(FLOAT(N)) K1=N+KMAX DO 10 I=1,N 10 P(K1+I)=P1(I) DO 20 KK=KMIN,KMAX-1 K=KMAX+KMIN-KK N1=(N2+1)/2 CALL DOWNP(NMAX,N1,N2,K,P) 20 N2=N1 DX1=DX*2**(KMAX-KMIN) CALL WI(NMAX,N1,KMIN,KMAX,DX,DX1,P,W) DO 30 K=KMIN+1,KMAX N2=2*N1-1 DX1=DX1/2. CALL AKCO(M+5,KMAX,K,AK1) CALL AKIN(M+6,AK1,AK2) CALL WCOS(NMAX,N1,N2,K,W) CALL CORR(NMAX,N2,K,M,1,DX1,P,W,AK1) CALL WINT(NMAX,N2,K,W) CALL CORR(NMAX,N2,K,M,2,DX1,P,W,AK2) 30 N1=N2 DO 40 I=1,N 40 W(I)=W(K1+I) RETURN END SUBROUTINE DOWNP(NMAX,N1,N2,K,P) DIMENSION P(NMAX) K1=N1+K-1 K2=N2+K-1 DO 10 I=3,N1-2 I2=2*I+K2 10 P(K1+I)=(16.*P(I2)+9.*(P(I2-1)+P(I2+1))-(P(I2-3)+P(I2+3)))/32. P(K1+2)=0.25*(P(K2+3)+P(K2+5))+0.5*P(K2+4) P(K1+N1-1)=0.25*(P(K2+N2-2)+P(K2+N2))+0.5*P(K2+N2-1) RETURN END SUBROUTINE WCOS(NMAX,N1,N2,K,W) DIMENSION W(NMAX) K1=N1+K-1 K2=N2+K DO 10 I=1,N1 II=2*I-1 10 W(K2+II)=W(K1+I) RETURN END SUBROUTINE WINT(NMAX,N,K,W) DIMENSION W(NMAX) K2=N+K DO 10 I=4,N-3,2 II=K2+I 10 W(II)=(9.*(W(II-1)+W(II+1))-(W(II-3)+W(II+3)))/16. I1=K2+2 I2=K2+N-1 W(I1)=0.5*(W(I1-1)+W(I1+1)) W(I2)=0.5*(W(I2-1)+W(I2+1)) RETURN END SUBROUTINE CORR(NMAX,N,K,M,I1,DX,P,W,AK) DIMENSION P(NMAX),W(NMAX),AK(0:M) K1=N+K IF(I1.EQ.2)GOTO 20 DO 10 I=1,N,2 II=K1+I J1=MAX0(1,I-M) J2=MIN0(N,I+M) DO 10 J=J1,J2 IJ=IABS(I-J) 10 W(II)=W(II)+AK(IJ)*DX*P(K1+J) RETURN 20 DO 30 I=2,N,2 II=K1+I J1=MAX0(1,I-M) J2=MIN0(N,I+M) DO 30 J=J1,J2 IJ=IABS(I-J) 30 W(II)=W(II)+AK(IJ)*DX*P(K1+J) RETURN END SUBROUTINE WI(NMAX,N,KMIN,KMAX,DX,DX1,P,W) DIMENSION P(NMAX),W(NMAX) COMMON /COMAK/AK(0:1100) K1=N+1 K=2**(KMAX-KMIN) C=ALOG(DX) DO 10 I=1,N II=K1+I W(II)=0.0 DO 10 J=1,N IJ=K*IABS(I-J) 10 W(II)=W(II)+(AK(IJ)+C)*DX1*P(K1+J) RETURN END SUBROUTINE AKCO(KA,KMAX,K,AK1) DIMENSION AK1(0:KA) COMMON /COMAK/AK(0:1100) J=2**(KMAX-K) DO 10 I=0,KA II=J*I 10 AK1(I)=AK(II) RETURN END SUBROUTINE AKIN(KA,AK1,AK2) DIMENSION AK1(KA),AK2(KA) DO 10 I=4,KA-3 10 AK2(I)=(9.*(AK1(I-1)+AK1(I+1))-(AK1(I-3)+AK1(I+3)))/16. AK2(1)=(9.*AK1(2)-AK1(4))/8. AK2(2)=(9.*(AK1(1)+AK1(3))-(AK1(3)+AK1(5)))/16. AK2(3)=(9.*(AK1(2)+AK1(4))-(AK1(2)+AK1(6)))/16. DO 20 I=1,KA 20 AK2(I)=AK1(I)-AK2(I) DO 30 I=1,KA-1,2 I1=I+1 AK1(I)=0.0 30 AK1(I1)=AK2(I1) RETURN END SUBROUTINE SUBAK(MM) COMMON /COMAK/AK(0:1100) DO 10 I=0,MM 10 AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.) RETURN END |
回复此楼» 猜你喜欢
网传南京林业大学加拿大博士首聘期不合格，自杀
已经有8人回复
-
准研三生态毒理学，爱化学生物新能源碳中和，电池等，25博士岗位求捞啦！
已经有14人回复
-
24/25申博，科研助理找导师
已经有7人回复
-
F06送审了吗？
已经有7人回复
-
基金委治打招呼顽疾越治越严重
已经有46人回复
-
南京林业大学老师首聘期考核不过，服药自杀
已经有3人回复
-
国社科系统提交后，下载申请书，里面的思路图不显示，是什么问题
已经有3人回复
-
求助，请问我是否该转专业
已经有9人回复
-
投稿竟然倒在了Scientific Reports上。
已经有9人回复
-
再谈一作和二作
已经有29人回复
100