小弟最近在學有限元素法
在下面程式中，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