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

[TERAB1T]

Очень нуждаюсь в помощи по написанию программы для вычисления кратных интегралов методом Монте-Карло Язык - любой, кроме С# и паскаля. Цель такова - вводим параметры (подинтегральную функцию (для реализации чего, как я понял, нужен интерпритатор), пределы интегрирования), n (для точности), получаем число-ответ.

@темы: Алгоритм, Монте-Карло, численное интегрирование, PHP

Comments

[mr Gray]

пхп этот должен запускаться при желании на любом компьютере, то есть без серверов и тд
Тут есть некоторые предложения на этот счет. Но лучше всего разместить программу на хосте и запускать по сети - тогда не придется копировать программу на все компьютеры, и, следовательно, не будет, например, проблем с ее обновлением. Ну а продемонстрировать программу преподавателю можно и через Zend, например.

По поводу интеграла по поверхности:
Если поверхность задается уравнением, то, задав случайный x, можно найти y из этого уравнения. И считать интеграл по найденным x и y. Только уравнение поверхности пользователь должен будет ввести, выразив y через x.

[TERAB1T]

По поводу интеграла по поверхности:
Если поверхность задается уравнением, то, задав случайный x, можно найти y из этого уравнения. И считать интеграл по найденным x и y. Только уравнение поверхности пользователь должен будет ввести, выразив y через x.

как это реализовать?

[mr Gray]

как это реализовать?

$g=lcg_value();
$x=$a+$g*($b-$a);
eval ("\$y = $func;");

Но при этом нужно учитывать, что $x может не попасть в область определения. Тогда вычисление $y вызовет ошибку, ее нужно перехватить через try и не учитывать эту итерацию.

[TERAB1T]

Не понимаю, чувствую себя тупым Что с границами интегрирования делать?

[mr Gray]

Сказали сделать так, чтоб область интегрирования могла быть и другой как это связано с Второе интегрирование у меня происходит в двойном интеграле строго в квадрате..
Что такое "второе интегрирование"? Что может быть другим, если задано строго?

Второе интегрирование у меня происходит в двойном интеграле строго в квадрате. Это область интегрирования? Если да, то здесь не нужно вычислять интеграл по поверхности. Это обычный двойной интеграл, где в качестве пределов по x и y берем координаты квадрата. Просто запрашиваете у пользователя координаты левого верхнего и правого нижнего углов квадрата, и интегрируете по ним.

PS: По поводу интегрирования по произвольной области - если я правильно понял условия задачи - оно здесь не понадобится. Но вообще для такого интеграла ИМХО нужно будет определять - лежит ли произвольная точка внутри области интегрирования, и если лежит - считать функцию. Наверное, как то так.

[TERAB1T]

Сказали сделать так, чтоб область интегрирования могла быть и другой как это связано с Второе интегрирование у меня происходит в двойном интеграле строго в квадрате..
Что такое "второе интегрирование"? Что может быть другим, если задано строго?

Насчет "второго интегрирования" я просто описАлся. Есть двойной интеграл, в нем конкретно заданы верхняя и нижняя границы интегрирования и x, и y, то есть мне остается только вписать их в соответсвующие поля, вписать функцию и нажать считать. Я же имею ввиду интегралы посложнее, пример: . Как сделать так, чтобы в этой программе считались такие интегралы, без лишних предварительных действий вне программы?

[TERAB1T]

Неужели никто не знает?

[mr Gray]

Неужели никто не знает?
Ваш преподаватель знает, и наверняка рассказывал про это на лекциях или при выдаче задания, если сейчас стоит такая задача.

С моей точки зрения, решение может выглядеть примерно так:
Чтобы решить любую задачу по программированию, ее сначала нужно формализовать. А тут непонятно, в каком виде поступают входные данные. Область может задаваться пересечением прямых, неравенством, несколькими неравенствами и т.д. Сперва определите - как должна задаваться область. Если неравенством, вида sqr(x)+sqr(y)<=25 (т.е. обязательно меньше или равно), то пользователь может выразить из этого неравенства x и ввести в поле для ввода выражение: sqrt(25-sqr(y)). На каждом шаге программы вы генерируете случайный y. Пусть на текущем шаге сгенерирован y1. Считаете x_max = sqrt(25-sqr(y1)). Теперь выбираете случайный x1, который меньше или равен x_max. По точкам x1 и y1 считаете подинтегральное выражение. Плюсуете его и переходите к следующему шагу и т.д. Еще нужно учесть, что при некоторых y1 выражение, введенное пользователем, может не посчитаться (например, если будет извлекаться корень из отрицательного числа). Тогда нужно проигнорировать этот шаг и повторить его, т.е. сгенерировать другой y1.

[TERAB1T]

Неужели никто не знает?
Ваш преподаватель знает, и наверняка рассказывал про это на лекциях или при выдаче задания, если сейчас стоит такая задача.

Преподаватель, который это задал - преподователь матана, а не информатики. Сказал - вот вы программисты, так что делайте как хотите

[TERAB1T]

Спасибо. Сейчас буду пробовать. А как такое сделать с 3+ кратными интегралами?

[mr Gray]

Преподаватель, который это задал - преподователь матана, а не информатики. Сказал - вот вы программисты, так что делайте как хотите
То-то я смотрю у вас вычислительная математика какая-то странная
Тогда советую выделить время, подойти к преподавателю и объяснить, что предложенные им условия задачи даны в слишком общем виде, и что нужно их уточнить.

А как такое сделать с 3+ кратными интегралами?
Если область интегрирования тоже задается как одно неравенство вида x1 <= F (x2, x3, x4 ...), то так же, но только на каждом шаге генерировать случайным образом уже несколько переменных: x2, x3, x4, ... Потом, используя эти переменные, считать x1_max = F (x2, x3, x4 ...). Потом генерировать x1 как случайное число, меньшее или равное x1_max. И т.д.

[TERAB1T]

mr Gray C этим случаем сейчас буду пробовать. А есть идеи насчет интеграла второго в примере, где S - треугольник?

[TERAB1T]

С моей точки зрения, решение может выглядеть примерно так:
Чтобы решить любую задачу по программированию, ее сначала нужно формализовать. А тут непонятно, в каком виде поступают входные данные. Область может задаваться пересечением прямых, неравенством, несколькими неравенствами и т.д. Сперва определите - как должна задаваться область. Если неравенством, вида sqr(x)+sqr(y)<=25 (т.е. обязательно меньше или равно), то пользователь может выразить из этого неравенства x и ввести в поле для ввода выражение: sqrt(25-sqr(y)). На каждом шаге программы вы генерируете случайный y. Пусть на текущем шаге сгенерирован y1. Считаете x_max = sqrt(25-sqr(y1)). Теперь выбираете случайный x1, который меньше или равен x_max. По точкам x1 и y1 считаете подинтегральное выражение. Плюсуете его и переходите к следующему шагу и т.д. Еще нужно учесть, что при некоторых y1 выражение, введенное пользователем, может не посчитаться (например, если будет извлекаться корень из отрицательного числа). Тогда нужно проигнорировать этот шаг и повторить его, т.е. сгенерировать другой y1.

Что в этом случае делать с границами интегрирования? Что туда вводить?

[mr Gray]

Что в этом случае делать с границами интегрирования? Что туда вводить?
Вместо границ интегрирования будет область интегрирования. То есть пользователь будет вводить не границы интегрирования а функцию F(x2, x3, x4 ...). Например, пользователь может ввести "sqrt(25-sqr($y))" это будет означать, что область интегрирования задана неравенством "sqr(x)+sqr(y)<=25" - круг.

А есть идеи насчет интеграла второго в примере, где S - треугольник?
Пользователь может ввести три точки - вершины треугольника. На каждом шаге программа генерирует случайную точку, например (x1; y2). Проверяем - лежит ли эта точка внутри треугольника (например, считаем сумму площадей треугольников, которые получатся, если соединить сгенерированную точку отрезками с вершинами треугольника). Если площадь этих треугольников равна площади области интегрирования - то точка лежит в этой области. Тогда считаем интеграл, прибавляем его к сумме и переходим на следующую итерацию. Если точка не лежит в области интегрирования (исходном треугольнике), то повторяем эту итерацию, генерируя другую точку.

Возникнет проблема, связанная с тем, что мы не знаем, где расположен треугольник и какие случайные точки генерировать. Точки можно генерировать так: например, вершины исходного треугольника: (Ax; Ay), (Bx, By), (Cx, Cy). Выбираем минимальные и максимальные координаты:
Xmin = min(Ax, Bx, Cx);
Xmax = max(Ax, Bx, Cx);
И тогда можно на каждом шаге генерировать точку, у которой координата X находится между Xmin и Xmax.
С координатой 'Y' можно поступить также (хотя можно и немного по другому: считая пересечение сгенерированного X с треугольником и находя границы возможного Y, но это посложнее).

PS: пользователю придется посчитать точки пересечения прямых, чтобы узнать вершины треугольника. Иначе, в общем случае, программе придется решать систему из трех уравнений, вид которых не известен и может быть любым. Можно ввести ограничение: выражать во всех 3 прямых переменную Y. Но если прямая имеет вид X=10, то тут нечего выражать...

[TERAB1T]

Что в этом случае делать с границами интегрирования? Что туда вводить?
Вместо границ интегрирования будет область интегрирования. То есть пользователь будет вводить не границы интегрирования а функцию F(x2, x3, x4 ...). Например, пользователь может ввести "sqrt(25-sqr($y))" это будет означать, что область интегрирования задана неравенством "sqr(x)+sqr(y)<=25" - круг.

ну ок, но в формуле Монте-Карло самой-то есть еще b-a (причем дважды, учитывая что интеграл двойной), вот с ним я и запутался, не знаю, что туда вводить

[mr Gray]

Вот скрипт, вычисляющий двойной интеграл по поверхности методом Монте-Карло. Тестировал только на паре примеров, так что проверьте получше. Границы Xmin/max Ymin/max должны быть константами.
http://paste.org.ru/?7c2czf

PS: Пределы по X и Y (которые нужно поставить так, чтобы область определения поверхности интегрирования в них полностью помещалась) можно определять и программно с определенной точностью в определенном интервале. Хотя, конечно, нет никакой гарантии, что функция не определена за пределами этого интервала. Кроме того интервал может быть очень маленьким, меньше шага при его поиске. Так что я решил с этим не заморачиваться и предоставить выбор интервала пользователю.

[TERAB1T]

что-то у меня не получается, пробовал с примером, описанным выше (задание 3903), скорее всего что-то не так ввожу. В поле F(x, y) ввожу "(1/(sqrt(24+pow($x, 2)+pow($y, 2))))" в область интегрирования "pow($x, 2)+pow($y, 2)" "<=" "25" (пробовал выражать х через у, все равно ничего не получается), пределы: x - -5;5; y - -5;5; ответ около 13ти, хотя в у чебнике приближенный ответ 9,88. Что я делаю не так?

[mr Gray]

Задания отсюда решаются и ответ совпадает. Посчитайте интеграл 3903 вручную и напишите, что получится. Если все же ошибка в скрипте, то могло возникнуть переполнение.

[TERAB1T]

Скорее всего я что-то неправильно ввожу. Если не трудно, распишите что куда надо вводить например для решения первого примера по ссылке, что вы дали

[mr Gray]

распишите что куда надо вводить например для решения первого примера по ссылке, что вы дали
Для вычисления этого интеграла нужно вводить такие данные:

F(x, y) = $x*$y

Область интегрирования (укажите F1(x, y), знак и F2(x, y))
$x*$x+4*$y*$y <= 1

Количество итераций для метода 10000
Нижний предел для x -1
Верхний предел для x 1
Нижний предел для y -1
Верхний предел для y 1

Ответ: 0.12604340369705, а это примерно 1/8

[TERAB1T]

а почему <= если в задании просто =?

[mr Gray]

а почему <= если в задании просто =?
Предполагается, что <=. Просто в задании указали линию, которая очерчивает область, по которой нужно интегрировать.

$x*$x+4*$y*$y = 1 это окружность, то есть линия. Ее площадь = 0 и, следовательно, интеграл тоже = 0.
$x*$x+4*$y*$y <= 1 это круг, имеющий площадь pi*r*r, и по нему уже можно интегрировать.

[TERAB1T]

Пределы по X и Y (которые нужно поставить так, чтобы область определения поверхности интегрирования в них полностью помещалась) можно определять и программно с определенной точностью в определенном интервале. Хотя, конечно, нет никакой гарантии, что функция не определена за пределами этого интервала. Кроме того интервал может быть очень маленьким, меньше шага при его поиске. Так что я решил с этим не заморачиваться и предоставить выбор интервала пользователю.

Можно все же про это поподробнее? у меня в тестах все косяки именно с определением границ

[mr Gray]

Можно все же про это поподробнее? у меня в тестах все косяки именно с определением границ
Суть в том, что нужно найти максимальные и минимальные X и Y, для которых выполняется условие: F1(x, y) [<=, >=] F2(x, y). Эти X и Y и будут границами.
Искать можно разными методами. Вот самый простой (но и самый неэффективный):
1) Разбиваем плоскость XOY на сетку, с определенным шагом h (например, h=0.1). Сетка будет не бесконечной, поэтому ее размер может быть, например, 200 на 200 (с центром в начале координат)
2) Для каждой точки сетки смотрим, выполняется ли условие F1(x, y) [<=, >=] F2(x, y).
3) Ищем Xmin, Xmax - минимальную и максимальную координаты X точек сетки, для которых выполняется условие F1(x, y) [<=, >=] F2(x, y).
4) Для Y делаем то же самое.
5) В результате придется затратить (200/0.1)*(200/0.1) = 4 млн итераций.

PS: можно использовать и другие методы, например – для поиска первой точки использовать крупную сетку, а потом уменьшать h и др. Но тогда это будет уже целая отдельная задача.

[TERAB1T]

как это реализовать?

[TERAB1T]

mr Gray Такой вопрос -разбираю код вашей программы, сделал по аналогии тройной интеграл, все также работает, но сам я так и не понял по какому принципу. То есть я не увидел в коде, где использовалось бы хоть как-то выражение, которое определяет границы, можете прояснить?

[mr Gray]

TERAB1T
Используется только метод Монте-Карло, без всяких посторонних ухищрений.
Нам не нужны границы, нам нужно знать площадь основания (фигуры, заключенной в эти границы).
Например, у нас есть фигура, основание которой лежит в плоскости XOY и задано в виде неравенства F1(x, y) [<=, >=] F2(x, y). Бока фигуры параллельны оси Z а сверху ее ограничивает подынтегральная функция. Тогда программа выбирает прямоугольник на плоскости XOY и начитает создавать в ней случайные точки: какой процент точек попадет в основание фигуры, такой процент площади прямоугольника является площадью фигуры. Зная площадь фигуры основания и среднюю высоту фигуры, перемножаем их и получаем объем - интеграл.

PS: Если еще не сделали, автоматическое определение границ сделаю вечером. Пока много дел.

[mr Gray]

Вот скрипт, автоматически вычисляющий область определения и считающий интеграл:
paste.org.ru/?a38o3i

PS: 58-ой комментарий... Иду на рекорд...

[TERAB1T]

Спасибо, осталось последнее - вычислять площадь фигуры, ограниченной некоторыми прямыми.

[TERAB1T]

mr Gray Спасибо вам за помощь, очень помогли, больше впринципе ничего не надо, этого оказалось достаточно. Теперь задание приобретает немного другой оборот, поэтому создаю новую тему.