24小时热门版块排行榜    

Znn3bq.jpeg
查看: 543  |  回复: 1

kimileegdut

捐助贵宾 (小有名气)

[求助] 求助:迭代计算有问题

编的程序迭代部分代码如下,在运行的时候出现问题,无论开始给一维数组A3赋非0的任何值,程序都只是迭代计算2次,而且结果也不太对,想请教一下各位大神,这段代码中哪里出现了问题?事关毕业,麻烦各位了!

********************************************************
!1)赋初值给A3
    A3=1.D0
!
!2)计算U3XT,U3YT:含有待求系数列阵A3和边界位移列阵U3B
  itera_of_A3O
    U3X=MATMUL(EX,U3B)+MATMUL(DX,A3)
    U3Y=MATMUL(EY,U3B)+MATMUL(DY,A3)
    DO i=1,M
      DO t=1,M
        IF(i==t)THEN
          U3XT(i,t)=U3X(i)
          U3YT(i,t)=U3Y(i)
        ELSE
          U3XT(i,t)=0
          U3YT(i,t)=0
        END IF
      END DO
    END DO
!
!3)计算系数矩阵C1和C2:含有待求的系数列阵A3
    temp1C1=MATMUL(U3XT,temp1B1)
    C1A=MATMUL(temp1C1,U3B)
    temp3C1=MATMUL(U3XT,A11)
    C1B=MATMUL(temp3C1,A3)
    temp2C1=MATMUL(EXY,U3B)+MATMUL(DXY,A3)
    C1C=0.5*(1.+MIU)*MATMUL(U3YT,temp2C1)
    C1=-C1A-C1B-C1C
!
    temp1C2=MATMUL(U3YT,temp2B2)
    C2A=MATMUL(temp1C2,U3B)
    temp2C2=MATMUL(U3YT,A22)
    C2B=MATMUL(temp2C2,A3)
    C2C=0.5*(1.D0+MIU)*MATMUL(U3XT,temp2C1)
    C2=-C2A-C2B-C2C
!
!4)A1,A2:两个待求的系数矩阵
!A1
    A1A=A11-MATMUL(A12,A22N)
    temp1A1=MATMUL(A12,A22N)
    temp2A1=C2-B2
    temp3A1=MATMUL(temp1A1,temp2A1)
    A1B=C1-B1-temp3A1
    A1=MATMUL(A1A,A1B)
!
!A2
    A2B=C1-B1-MATMUL(A11,A1)
    A2=MATMUL(A12N,A2B)
!
!5)U1XT,U1YT,U2XT,U2YT,U3,U3XX,U3YY,U3XY
    U1X=MATMUL(EX,U1B)+MATMUL(DX,A1)
    U1Y=MATMUL(EY,U1B)+MATMUL(DY,A1)
    U2X=MATMUL(EX,U2B)+MATMUL(DX,A2)
    U2Y=MATMUL(EY,U2B)+MATMUL(DY,A2)
    U3=MATMUL(E,U3B)+MATMUL(D,A3)
    U3XX=MATMUL(EXX,U3B)+MATMUL(DXX,A3)
    U3YY=MATMUL(EYY,U3B)+MATMUL(DYY,A3)
    U3XY=MATMUL(EXY,U3B)+MATMUL(DXY,A3)
    DO i=1,M
      DO t=1,M
        IF(i==t)THEN
          U1XT(i,t)=U1X(i)
          U1YT(i,t)=U1Y(i)
          U2XT(i,t)=U2X(i)
          U2YT(i,t)=U2Y(i)
        ELSE
          U1XT(i,t)=0
          U1YT(i,t)=0
          U2XT(i,t)=0
          U2YT(i,t)=0
        END IF
      END DO
    END DO
!
!6)利用简单迭代法x=f(x)计算
    DD1=U1XT+0.5*MATMUL(U3XT,U3XT)+MIU*U2YT+0.5*MIU*MATMUL(U3YT,U3YT)
    DD2=(1.D0-MIU)*(U1YT+U2XT+MATMUL(U3XT,U3YT))
    DD3=U2YT+0.5*MATMUL(U3YT,U3YT)+MIU*U1XT+0.5*MIU*MATMUL(U3XT,U3XT)
    DD4=MATMUL(DD1,EXX)+MATMUL(DD2,EXY)+MATMUL(DD3,EYY)
    K1=MATMUL(DD4,U3B)
    DD5=MATMUL(DD1,DXX)+MATMUL(DD2,DXY)+MATMUL(DD3,DYY)
    K2=(1.D0/D0)*(q-2.D0*ks*U3)
    K3=K1+K2
!
!求DD5的逆阵DD5N
    DWZ_DD5=DD5
!
!生成单位阵
    EDD5=0
    DO i=1,M
      EDD5(i,i)=1
    END DO
    DD5N=EDD5
!调用子程序
    CALL gaussr_inverse(DWZ_DD5,M,M,DD5N,M,M)
!
!判断逆矩阵是否有效
    DWZ_DD5=MATMUL(DD5,DD5N)
    dm_DD5=0
    DO i=1,M
        DO t=1,M
            d_DD5=DWZ_DD5(i,t)-EDD5(i,t)
            IF(abs(d_DD5) > abs(dm_DD5) ) THEN
               dm_DD5=d_DD5
            END IF
        END DO
     END DO
     IF(ABS(dm_DD5)<error_inverse) GO TO 74
     WRITE(*,*)'Error of DD5 inverse'
74 CONTINUE
!
!迭代式
    A3I=-MATMUL(DD5N,K3)
    MAX_A3L=0
    DO i=1,M
      A3L=A3(i)-A3I(i)
      IF(ABS(A3L)>ABS(MAX_A3L))THEN
         MAX_A3L=A3L
      END IF
    END DO
    n_itera=n_itera+1
    IF(MAX_A3L<=itera_A3_error)EXIT itera_of_A3
    A3=A3I
  END DO itera_of_A3
回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐,对您同样有帮助:

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

kimileegdut

捐助贵宾 (小有名气)

itera_of_A3:DO
2楼2015-05-17 15:34:48
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 kimileegdut 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 博士申请 +5 星…… 2026-05-18 6/300 2026-05-18 23:49 by 糊糊涂涂好
[基金申请] 今年审到国自然15份,谈谈感受 +16 国自然国社科中 2026-05-17 16/800 2026-05-18 14:58 by gy116024
[基金申请] 青C资助名额大幅增加! +12 西葫芦炒鸡蛋 2026-05-13 16/800 2026-05-18 10:02 by Equinoxhua
[硕博家园] 我在等一个没有答案的答案 +3 Love_MH 2026-05-17 3/150 2026-05-18 02:22 by 竹林孤影
[找工作] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +4 ky2p12rrjj 2026-05-15 4/200 2026-05-17 19:47 by Equinoxhua
[硕博家园] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +6 l7k6xnh0yc 2026-05-14 7/350 2026-05-17 19:42 by Equinoxhua
[考研] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +4 xx7gd5zq4e 2026-05-15 6/300 2026-05-17 19:36 by Equinoxhua
[考研] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +6 cjf4bx70cj 2026-05-14 7/350 2026-05-17 18:49 by Equinoxhua
[考博] 2026博士还有哪些学校有名额 +6 小王求读研 2026-05-15 7/350 2026-05-17 16:54 by 知音湖畔
[硕博家园] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +4 v9tggjlwd0 2026-05-15 4/200 2026-05-17 08:06 by 11n4dfd8yn
[硕博家园] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +5 cjf4bx70cj 2026-05-14 7/350 2026-05-17 06:55 by 11n4dfd8yn
[找工作] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +3 k37jurhrau 2026-05-16 3/150 2026-05-17 01:37 by ue3ir18jc3
[考研] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +4 x0mp7owy2b 2026-05-15 4/200 2026-05-17 00:25 by ue3ir18jc3
[基金申请] 请问大佬b0816评完了吗 +4 市民华南虎 2026-05-12 8/400 2026-05-16 19:54 by Equinoxhua
[有机交流] 求助2,4-二氯-5-嘧啶甲醛的合成方法 20+3 光吃不拉 2026-05-14 6/300 2026-05-16 19:46 by Equinoxhua
[有机交流] 如何实现卤原子转化 +3 BT20230424 2026-05-15 5/250 2026-05-16 16:20 by czyzsu
[考研] 售SCI一区T0P文章,我:8.O.5.5.1.O.5.4,科目齐全,可+急 +4 l7k6xnh0yc 2026-05-14 6/300 2026-05-16 11:29 by h3oerqvkv9
[文学芳草园] 风把牡丹吹跑了 +5 myrtle 2026-05-12 9/450 2026-05-15 15:27 by myrtle
[考博] 材料类只有一篇综述能申博么 +4 乐逍遥谷 2026-05-13 4/200 2026-05-14 12:05 by zhyzzh
[论文投稿] 求助大佬sci投稿哪个好中 +3 江沅188 2026-05-12 4/200 2026-05-13 14:35 by 江沅188
信息提示
请填处理意见