[問題] 輸出值顯示問號

看板Fortran作者dregsgod (五胚)時間5年前 (2019/03/09 04:25)推噓1(1推 0噓 0→)

留言1則, 1人參與討論串1/1

小弟最近在學有限元素法在下面程式中，n<24內都可以正常顯示輸出的值但一旦n超過24時，輸出的值除了原本給定的邊界值計算出來的值均會變成??????????? 想請問版上的大神們有遇過這種情形嗎? PROGRAM FEM_HW1 IMPLICIT NONE real*8,parameter :: min = 0.d0 !Boundary real*8,parameter :: Max = 1.d0 !Boundary integer,parameter :: n = 24 integer :: i real*8 :: h,a,b,dx,x real*8 :: aa11,aa12,aa22,bb1,bb2 real*8 :: a11,a12,a22,b1,b2 real*8 :: ans real*8 :: gl(0:n) real*8 :: BB(0:n),DD(0:n),AA(0:n),CC(0:n) = 0 real*8 :: y(0:n) aa11(x,b,dx) = x*((b-x)/dx)**2 !Q = x aa22(x,a,dx) = x*((x-a)/dx)**2 !Q = x aa12(x,a,b,dx) = x*(a-x)/dx*(x-b)/dx !Q = x bb1(x,b,dx) = -x*(b-x)/dx !G = -x bb2(x,a,dx) = -x*(x-a)/dx !G = -x open(unit = 9,FILE='FEM_HW1.txt') dx = (Max - min)/dble(n) do i = 0,n gl(i) = dble(i)*dx+min end do do i = 0,n-1 h = (gl(i+1)-gl(i))/3.d0 call simpson(gl(i),gl(i+1),aa11(gl(i),gl(i+1),dx) C,aa11(gl(i)+h,gl(i+1),dx),aa11(gl(i)+h*2,gl(i+1),dx) C,aa11(gl(i)+h*3,gl(i+1),dx),ans) a11 = 1.d0/dx - ans call simpson(gl(i),gl(i+1),aa22(gl(i),gl(i),dx) C,aa22(gl(i)+h,gl(i),dx),aa22(gl(i)+h*2,gl(i),dx) C,aa22(gl(i)+h*3,gl(i),dx),ans) a22 = 1.d0/dx - ans call simpson(gl(i),gl(i+1),aa12(gl(i),gl(i),gl(i+1),dx) C,aa12(gl(i)+h,gl(i),gl(i+1),dx),aa12(gl(i)+h*2,gl(i),gl(i+1),dx) C,aa12(gl(i)+h*3,gl(i),gl(i+1),dx),ans) a12 = -1.d0/dx - ans call simpson(gl(i),gl(i+1),bb1(gl(i),gl(i+1),dx) C,bb1(gl(i)+h,gl(i+1),dx),bb1(gl(i)+h*2,gl(i+1),dx) C,bb1(gl(i)+h*3,gl(i+1),dx),ans) b1 = -ans call simpson(gl(i),gl(i+1),bb2(gl(i),gl(i),dx) C,bb2(gl(i)+h,gl(i),dx),bb2(gl(i)+h*2,gl(i),dx) C,bb2(gl(i)+h*3,gl(i),dx),ans) b2 = -ans BB(i+1) = BB(i+1) + a12 DD( i ) = DD( i ) + a11 AA( i ) = AA( i ) + a12 DD(i+1) = DD(i+1) + a22 CC( i ) = CC( i ) + b1 CC(i+1) = CC(i+1) + b2 end do BB(0) = 0.d0 DD(0) = 1.d0 AA(0) = 0.d0 CC(0) = 0.d0 !B.C. BB(n) = 0.d0 DD(n) = 1.d0 AA(n) = 0.d0 CC(n) = 0.d0 !B.C. CC( 1 ) = CC(1) - BB(1)*CC(0) BB( 1 ) = 0.d0 CC(n-1) = CC(n-1) - AA(n-1)*CC(n) AA(n-1) = 0.d0 call SY(1,n-1,BB(1:n-1),DD(1:n-1),AA(1:n-1),CC(1:n-1)) y = CC do i = 0,n write(*,"((F12.8))")y(i) write(9,"((F12.8))")y(i) end do close(9) STOP END SUBROUTINE simpson(min,Max,a,b,c,d,I) IMPLICIT NONE real*8 :: min,Max,a,b,c,d real*8 :: I I = (Max-min)*(a+b*3+c*3+d)/8 RETURN END SUBROUTINE SY(IL,IU,BB,DD,AA,CC) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION AA(1), BB(1), CC(1), DD(1) LP = IL + 1 DO 10 I = LP, IU R = BB(I)/DD(I-1) DD(I) = DD(I) - R*AA(I-1) 10 CC(I) = CC(I) - R*CC(I-1) CC(IU) = CC(IU)/DD(IU) DO 20 I = LP,IU J = IU - I + IL 20 CC(J) = (CC(J)-AA(J)*CC(J+1))/DD(J) RETURN END -- ※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.116.20.17 ※ 文章網址: https://www.ptt.cc/bbs/Fortran/M.1552076717.A.FB2.html