Метод Монте-Карло для вычисления интеграла


Рассмотренные методы называются детерминированными, то есть лишенными элемента случайности.
Методы Монте–Карло (ММК) – это численные методы решения математических задач с помощью моделирования
случайных величин. ММК позволяют успешно решать математические задачи, обусловленные вероятностными
процессами. Более того, при решении задач, не связанных с какими-либо вероятностями, можно искусственно
придумать вероятностную модель (и даже не одну), позволяющую решать эти задачи. Рассмотрим вычисление
определенного интеграла.
J = интеграл от а до b функции f(x)
Я рекомендую Вам два способа для вычисления интегралов методом Монте–Карло в Python
  • Статистический метод
  • При вычислении этого интеграла по формуле прямоугольников интервал разбивается на
    одинаковых интервалов, в серединах которых вычисляются значения подынтегральной функции.
    Вычисляя значения функции в случайных узлах, можно получить более точный результат:

    J = ((b-a)/N)*SUM(f(xj),j=1..N)

    где xj= a + r(b-a)

    Здесь r — случайное число, равномерно распределенное на интервале [0,1].
    Погрешность вычисления интеграла ММК значительно больше,
    чем у ранее рассмотренных детерминированных методов.
    Код в Python  В этом случае a = 1.2, b=2.2

    
    def Morte_Karlo():
        sum=0.0
        for i in range(n+1):
            r=ran.random()
            x=1.2 + r
            sum += function(x)
        return sum/n </em>
    
  • Обычный метод
  • J = (D*nf)/N

    где (xi,yi) — пар случайных чисел, равномерно распределенных в прямоугольнике [a,b]x[c,d]
    содержащем искомый интеграл (a<=xi<=b,0<=yi<=d,i=1,..,N), nf — количество точек, удовлетворяющих условию
    yi<=f(xi), D — площадь прямоугольника [a,b]x[c,d].
    Код в Python

    В этом случае a = 1.2, b=2.2.

    Выбираю 0 <= c <= f(a), d > max f(x)
    
    			def Morte_Karlo():
    			    c=0.68
                                d=9.98
    	                    D=(b-a)*(d-c)
    		            nf=0
    		            for i in range(n+1):
    		                x=ran.uniform(1.2,2.2)
    		                y=ran.uniform(c,d)
    		                if (y<=function(x)):
    		                    nf=nf+1
    		            return (D*nf)/n
    

0 комментариев

Только зарегистрированные и авторизованные пользователи могут оставлять комментарии.