사다리꼴 또는 심슨 규칙을 사용한 적분 계산

Question

기본적으로 할당은 "f (x) = 0 ....에 대한 해결책 찾기"라고 말합니다. x = 0에서 x - 0.45까지 f (x) = integrate (1/sqrt (2 * pi)) exp (-t^2/2) dt. 죄송합니다. 여기에 정수를 넣는 방법을 찾을 수 없었습니다.

내 교수. 우리에게 시작하는 법을 알려주었습니다. 여기

start x 
do i = 1,2,3... 
    fp = 1/sqrt(2*pi)exp(-x^2/2) 
    f = use trap,or simpson's rule to find the integration than subtract 0.45 
    x = x - (f/fp) 
end do 

내가 내 실수가 될 수있는 곳을 찾기 위해 인쇄했다

program main 
implicit none 

integer :: n, k, i 
double precision :: h, a, fp, f, x1, x2, pi, blub, integ, e, dx, j, m 

a = 0 
n = 25 
x1 = 0.5 
pi = 4.0d0 * atan(1.0d0) 

do i = 1, n 
    e = exp((-x1)**2.0d0/2.0d0) 
    print*, e 

    fp = (1.0d0/sqrt(2.0d0 * pi))* e 

    !start of the trapezoid 
    dx = (x1-a)/n !find dx by subtracting starting point from endpoint and divide by number of steps taken 
    m = (blub(a) + blub(x1))/2 !find the mean value of the integral 


    j = 0 
     do k=1, n-1 
      h = i 
      j = j + blub(h) !calculate the each step of the integral from one to n and add all together. 
     end do 

    integ = dx*(m+j) !result of the trapezoid method 

    print*, "integ: ", integ 
    f = integ - 0.45 
    a = x1 
    x1 = a - (f/fp) 
    print*, "x1: ",x1 
    print*, "a: ",a 
    print*, "f: ",f 
    print*, "fp: ", fp 
    print*, (x1-a)/n 
end do 

stop 
end 

double precision function blub (x) 
double precision :: x 

blub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0) 

return 
end 

을 한 것입니다. 나는 내 dx가 n보다 작은 x1을 기반으로 작아지고 있음을 깨달았다. 그 때문에 내 integ는 0.45에 거의 영향을 미치지 않는 e^-310과 같은 숫자가되고, 따라서 f는 동일하게 유지되고 x1이 뒤죽박죽이됩니다.

누군가 내게 설명 할 수 있다면 정말 좋을 것입니다. 내가 실수를 고칠 수있는 방법.

편집 : 나는 사다리꼴 부분으로 실수를했다고 확신하지만, 나는 단지 어디에 있는지 모른다. f는 모든 루프마다 변경되지만 dx는 너무 작아서 변경되지 않습니다.

Answer 1

오류 중 하나는 다음과 같습니다. 나는 그것이 하나의 주장하지 않습니다

double precision function blub (x) 
double precision :: x 

blub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0) 

return 
end 

당신이 gfortran -Wall -std=f95 -c fun.f90로 컴파일하면 당신에게 벌어지고있는 이상한 뭔가를 경고 할 수

fun.f90:5.76: 

lub(x) = (1.0d0/sqrt(2.0d0 * (4.0d0 * atan(1.0d0))))*exp((-x)**2.0d0/2.0d0) 
                      1 
Warning: Obsolescent feature: Statement function at (1) 
fun.f90:1.30: 

double precision function blub (x) 
           1 
Warning: Return value of function 'blub' at (1) not set 

이

를 얻을. 원래 함수 내에서 올바른 값을 제공하는 대신 새로운 명령문 함수 (더 이상 사용되지 않는 Fortran 기능)를 작성했습니다. 당신은 사용해야

double precision function blub (x) 
    implicit none 
    double precision, intent(in) :: x 

    blub = (1/sqrt(8 * atan(1.0d0))) * exp((-x)**2/2) 
end function 

당신은 당신이 end 전에 넣어 stopreturn 끝을 잊을 수있다. 그들을 사용할 이유가 없습니다.

더 많은 디버깅 컴파일러 플래그 (예 : -g -fbacktrace -fcheck=all)를 사용하고 모든 함수와 서브 루틴을 모듈에 넣거나 내부로 만들 것을 권장합니다 (contains 사용).