real q,g,sum,sun,m,n,h
dimension X(1000),Y(1000),l(1000),z(1000),A(1000),b(1000),f(1000)
OPEN(UNIT=6,FILE='simpson.txt')
write(*,*) '반지름을 입력하세요'
read *, r
write(*,*) 'simpson 승수 s1를 입력하세요'
read *,n
write(*,*) 'simpson 승수 s2를 입력하세요'
read *,m
d=2./3.*r**3*3.1416
x(1)=0
h=r*2/n
sum=0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
do 888 i=1,n+1
y(i)=ver(r,x(i))
sun=0
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ! 반원 넓이 구하기!
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
j=1
do while (j<m+1)
z(1)=0
g=y(i)*2/m
a(j)=ver(y(i),z(j))
z(j+1)=z(j)+g
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ! Simpson 승수 !
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (MOD(j,2)==0) then
f(j)=4
ELSE
f(j)=2
end if
f(1)=1
f(m+1)=1
b(j)=a(j)*f(j)*g/3
sun=sun+b(j)
j=j+1
end do
continue
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ! Simpson 승수 !
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (MOD(i,2)==0) then
l(i)=4
ELSE
l(i)=2
end if
l(1)=1
l(n+1)=1
x(i+1)=x(i)+h
sum=sum+sun*l(i)*h/3
v=d-sum
888 continue
write(*,*) '실제 Volume=',d,' simpson공식 계산 =',sum
write(*,*) '오차는=', v,'입니다'
stop
end
function ver(t,u) result(dy)
intent(in) :: t,u
dy=sqrt(t**2-(u-t)**2)
end function ver
r=10으로 주고 m=100 으로 주면 에러가 뜨네요..
m=97~100까지 에러가 뜨고 다른 수를 넣으면 되고.. 무엇이 문제일까요 ㅠㅠ
r을 m등분하는 부분에서 에러가 뜨는것 같은데용 아무리 생각을 해봐도 sqrt안의 계산값은 0보다 작을 수가 없는데요...
r**2-(x-r)**2이니까요.. x는 1부터 r까지 증가를 하는데 아우
첫댓글 sqrt안이 0보다 작은 결과가 나오면 그렇습니다. r=10으로 주고 m=100 으로 주면 에러가 뜨는것은 분명 0보다 작은 값이 나오기 때문입니다. 기회가 되시면 손으로 직접 풀어보시면 에러를찾을수도 있습니다. 저도 이런 에러가 나왔을때는 손으로 직접풀어서, 계산과정중 어떤 결과가 나오는지를 직접확인하곤 했었습니다.
곰곰히 생각해보니.. 1/3= 0.333333333333인데 만약 반올림 되어서 이런 숫자가 0.334로 된다면 결국 r보다 커지게 되어 -가 나오는것 같습니다. 그런데그것을 해결할 방법을 모르겠네요... 임시 방편으로 절대값을 취해주니 에러는 안뜨는데 홀수에서 값이 이상하게 나오네요
그럼 F25.20 이런식으로 자릿수를 정하거나, IF문을 이용해보세요. r보다 커질경우 어떻게 처리를 해야 할지를 알려주면 됩니다. 또한, 그런경우가 생길경우 write문을 이용해서 화면에 값을 표시히주는 방법도 있습니다.