Численная интеграция с использованием Fortran 90 - PullRequest
Новая
       61

Численная интеграция с использованием Fortran 90

1 голос
/ 06 апреля 2019

Я пытаюсь использовать правило Симпсона для оценки интеграла sqrt (1-x ^ 2) в интервале от -1 до 1. Однако сумма, представленная переменной "s", в коде, который я разработано, не сходится вообще к пи за 2. Я абсолютный новичок в фортране и программировании в целом, поэтому, пожалуйста, потерпите меня. Что я делаю не так? 

PROGRAM integral

REAL*8 :: x
REAL*8 :: h
REAL*8 :: fodd
REAL*8 :: feven
REAL*8 :: simpson
REAL*8 :: s

x = -1
s = 0
simpson = 0
h = 0

   DO WHILE (x<=1)

     fodd = sqrt(1-(x+(2*h+0.1))**(2))
     feven = sqrt(1-(x+2*h)**(2))
     simpson = 4*fodd + 2*feven
     s = s + simpson*(h/3)
     WRITE(*,*) x,h, fodd, feven, simpson, s
     h = 0.1
     x = x + h
     END DO

END PROGRAM

Вот ссылка на вывод, который он генерирует: https://pastebin.com/mW06Z6Lq

Поскольку этот интеграл составляет лишь половину площади круга с радиусом 1, он должен сходиться к пи над 2, но он значительно превосходит это значение. Я думал об уменьшении шага для точности, но это не проблема, так как он превысил ожидаемое значение еще больше, когда я попробовал это.

1 Ответ

4 голосов
/ 07 апреля 2019

Вы должны проверить, что происходит, когда вы получаете с x близко к 1. Вы, конечно, не можете использовать DO WHILE (x<=1), потому что когда x==1, тогда x+2*h выше 1, и вы вычисляете квадратный корень из отрицательного числа.

Я предлагаю вообще не использовать DO WHILE.Просто установите количество делений, вычислите шаг, используя размер шага в качестве длины интервала / число делений, а затем используйте 

sum = 0
h = interval_length / n
x0 = -1

do i = 1, n
  xa = (i-1) * h + x0 !start of the subinterval
  xb = i * h + x0 !end of the subinterval    
  xab = (i-1) * h + h/2 + x0 !centre of the subinterval

  !compute the function values here, you can reuse f(xb) from the
  !last step as the current f(xa)


  !the integration formula you need comes here, adjust as needed
  sum = sum + fxa + 4 * fxab + fxb
end do

! final normalization, adjust to follow the integration formula above
sum = sum * h / 6

Обратите внимание, что вышеприведенное гнездо цикла написано очень общим способом, а неспецифичные для правила Симпсона.Я просто предположил, что h постоянен, но даже это можно легко изменить.По правилу Симпсона его можно легко оптимизировать.Вы, конечно, хотите только две оценки функции за интервал.Если это школьное задание, и вам нужно считать точки нечетно-четными, а не центральными, вы должны сами это скорректировать, это очень просто.

...