Код Фортана для интеграции Монте-Карло в граничных точках a и b - PullRequest
Новая
       138

Код Фортана для интеграции Монте-Карло в граничных точках a и b

0 голосов
/ 01 марта 2019

Я понимаю, что моделирование методом Монте-Карло предназначено для оценки площади путем построения случайных точек и вычисления соотношения между точками за пределами кривой и внутри кривой.

Я хорошо рассчитал значение пи, предполагая, что радиус кривой равенбудьте едины.

Вот код 

program pi
implicit none

integer :: count, n, i
real :: r, x, y
count = 0
n=500
CALL RANDOM_SEED
DO i = 1, n
 CALL RANDOM_NUMBER(x)
 CALL RANDOM_NUMBER(y)
 IF (x*x + y*Y <1.0) count = count + 1
END DO
r = 4 * REAL(count)/n
print *, r
end program pi

Но, чтобы найти интеграцию, Учебник говорит применить ту же идею.Но я заблудился о том, как написать код, если я хочу найти интеграцию 

f(x)=sqrt(1+x**2) over a = 1 and b = 5

Раньше, когда радиус был равен единице, я предполагал, что точка попадает внутрь по условию x * 2 + y ** 2но как я могу решить выше одного? enter image description here

Любая помощь очень полезна

1 Ответ

0 голосов
/ 06 марта 2019

Сначала я напишу код, а затем объясню: 

Program integral
implicit none
real f
integer, parameter:: a=1, b=5, Nmc=10000000   !a the lower bound, b the upper bound, Nmc the size of the sampling (the higher, the more accurate the result)
real:: x, SUM=0

do i=1,Nmc                  !Starting MC sampling
  call RANDOM_NUMBER(x)     !generating random number x in range [0,1]
  x=a+x*(b-a)               !converting x to be in range [a,b]
  SUM=SUM+f(x)              !summing all values of f(x). EDIT: SUM is also an instrinsic function in Fortran so don't call your variable this, I named it so, to illustrate its purpose
enddo

print*, (b-a)*(SUM/Nmc)     !final result of your integral
end program integral

function f(x)           !defining your function
  implicit none
  real, intent(in):: x
  real:: f

  f=sqrt(1+x**2)
end function f

Итак, что происходит:

Интеграл integ можно записать как integ2.где:

g(x)

(это g (x) - равномерное распределение вероятности переменной x в [a, b]).И мы можем записать интеграл как:

inter3

, где this.

Итак, наконец, мы получаем, что интеграл должен быть:

final

Итак, все, что вам нужно сделать, это сгенерировать случайное число в диапазоне [a, b], а затем вычислить значение вашей функции для этого x.Затем сделайте это много раз (раз Nmc) и вычислите сумму.Затем просто разделите на Nmc, чтобы найти среднее, а затем умножьте на (ba).И это то, что делает код.

В интернете есть много материала для этого. вот один пример, который довольно хорошо это визуализирует

РЕДАКТИРОВАТЬ: Второй способ, такой же, как метод Пи: 

Nin=0                    !Number of points inside the function (under the curve)
do i=1,Nmc
  call random_number(x)
  call random_number(y)
  x=a+x*(b-a)
  y=f_min+y(f_max-f_min)
  if (f(x)<y) Nin=Nin+1
enddo
print*, (f_max-f_min)*(b-a)*(real(Nin)/Nmc)

Все это вы могли бы затемзаключите его во внешний цикл do, суммируя (f_max-f_min) (ba) (real (Nin) / Nmc), и в конце выведите его среднее значение.В этом примере вы по сути создаете вмещающий прямоугольник от a до b (измерение x) и от f_min до f_max (измерение y), а затем делаете выборку точек внутри этой области и подсчитываете точки в функции (Nin). Очевидно, вам нужно знать минимальное (f_min) и максимальное (f_max) значение вашей функции в диапазоне [a, b].В качестве альтернативы вы можете использовать произвольно низкие / высокие значения для вашего f_min f_max, но тогда вы будете тратить много очков и ваша ошибка будет больше.

...