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

Я пытаюсь использовать правило Симпсона для вычисления интеграла 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, но он намного превосходит это значение. Я думал о том, чтобы сделать шаг меньше для точности, но это не проблема, так как он даже превзошел ожидаемое значение, когда я попробовал это.


fortran numerical-integration fortran90 simpsons-rule
person Insight    schedule 06.04.2019    source источник
comment
Вы обязательно должны установить h перед циклом. Это не главная проблема здесь, но это важно для будущего. Кроме того, вы действительно должны сообщить нам, какие именно результаты вы получаете и какие именно результаты вы ожидаете взамен.   -  person Vladimir F    schedule 06.04.2019
comment
Когда говорят о сходимости при численном интегрировании, обычно имеют в виду сходимость при h->0. Теперь у вас есть только h = 0,1, вы должны переходить к все более и более низким значениям. Кроме того, поскольку вы утверждаете, что являетесь новичком, обратите внимание, что real*8 не является стандартным Fortran stackoverflow.com/questions/838310/fortran-90-kind-parameter stackoverflow.com/questions/3170239/   -  person Vladimir F    schedule 07.04.2019
comment
Шаг устанавливается дважды, потому что мне нужно, чтобы он начинался с нуля, чтобы правильно оценить нечетные члены. Размер шага здесь не проблема, что-то еще не работает.   -  person Insight    schedule 07.04.2019
comment
Да, я делал уже говорил вам, что это не главная проблема здесь. Но, поверьте, размер шага h действительно должен быть установлен перед циклом и пусть зафиксирован. Пожалуйста, не ссылайтесь ни на какой pastebin. Поместите важную информацию прямо в свой вопрос. Вопросы здесь должны быть доступны для других даже тогда, когда внешние ресурсы перестанут работать в будущем.   -  person Vladimir F    schedule 07.04.2019

Ответы (1)


arrow_upward
4
arrow_downward

Вы должны проверить, что происходит, когда 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 является константой, но даже это можно легко изменить. Для правила Симпсона его можно легко оптимизировать. Вы, конечно, хотите только два вычисления функции за интервал. Если это школьное задание, и вам нужно рассматривать точки как нечетные-четные, а не как центральную границу, вы должны настроить это самостоятельно, это очень просто.

person Vladimir F    schedule 06.04.2019
comment
Спасибо за вашу помощь, но, похоже, это не работает. Я заменил значение 2 для interval_length, 10 для n и для функций, вычисляемых в точках, которые я только что написал, например, sqrt(1-(xa)**(2)) . Но он возвращает что-то близкое к 15 для значения s, и если я попытаюсь сделать меньшие шаги, изменив n, он вернет еще большие значения. - person Insight; 08.04.2019
comment
@Insight Я забыл разделить на n, это должно было быть довольно очевидно, я написал код без какого-либо тестирования, это ваша задача. - person Vladimir F; 08.04.2019