Неустойчивый метод Монте-Карло

[Geckelberryfinn]

Пытаюсь вычислить интеграл по треугольнику. Сначала вычислял его по квадратурным формулам Гаусса. Но потом, засомневавшись в правильности его работы, решил вычислить по тому же треугольнику интеграл от той же области методом Монте-Карло.
Так вот какая удивительная штуковина выходит: для простой функции вроде f(x,y)=x^2-y^2 метод работает хорошо. А вот для той функции, которая нужна мне (это полином второй степени - базисная функция в треугольнике), метод дает каждый раз разные результаты, причем разброс достаточно существенный. Увеличение количества выкидываемых точек не приводит к стабилизации.
Есть какие-нибудь идеи о том, в чем может скрываться причина такого странного поведения?
Вот код метода

N:=5000;
pintri:=1;   fmean:=0;
xmax:=max(x.x,max(y.x,z.x));
xmin:=min(x.x,min(y.x,z.x));
ymax:=max(x.y,max(y.y,z.y));
ymin:=min(x.y,min(y.y,z.y));
randomize;
while pintri<N do
begin
  cp.x:=xmin+random*(xmax-xmin);
  cp.y:=ymin+random*(ymax-ymin);
  inc(pintri);
  if InTri(cp,x,y,z) then //функция проверки принадлежности точки треугольнику - проверена.
    begin
      curv:=GetPolynom(cp.X,cp.Y,aa,bb,cc,dd,ee,ff,BlockNo)*
    (
        GetDerivative(cp.x,cp.y,BlockNo,a1,b1,c1,d1,e1,f1)
    );//GetDerivative - возвращает константу
      fmean:=fmean+curv;
    end;
end;
 Result:=(xmax-xmin)*(ymax-ymin)*(fmean)/N;

[Brukvalub]

Geckelberryfinn писал(а):

А вот для той функции, которая нужна мне (это полином второй степени - базисная функция в треугольнике)

Не по теме, но не могу скрыть удивления: полином второй степени по треугольнику на плоскости интегрируется аналитически ( то есть получается формула, вычисляющая истинное значение интеграла), зачем для точно решаемой задачи применять приближенные методы?

[Geckelberryfinn]

Это сделано для универсальности. Чтобы потом можно было просто увеличивать степень полинома, не вычисляя заново аналитически значение интеграла.

[незваный гость]

Вам не трудно поместить текст недостающих функций (и необходимые константы, вроде вершин треугольника)? Мне бы хотелось самому погонять…

Добавлено спустя 3 минуты 18 секунд:

Одна из сразу видимых проблем — это то, что Вы считаете и точки, не попавшие в треугольник. То есть, в вычислении у Вас может участвовать заметно меньше точек, чем Вы задаете.

[Geckelberryfinn]

Был вариант, когда я и не учитывал точки, не попавшие в треугольник. Тот же результат. Правда при количестве точек 5e+7 вроде более менее одинаковые результаты. Но это как то не спортивно...
недостающие функции:

type RealType=double;
type Dot=record
   x,y:RealType;

function GetPolynom(xx,yy,ga,gb,gc,gd,ge,gf:RealType;gblockno:integer;
                     IsBound:boolean=false): RealType;
begin

Result:=ga+gb*xx+gc*yy+gd*xx*yy+ge*sqr(xx)+gf*sqr(yy)

end;

function OnLine(TestPoint,X1,X2:Dot):boolean;
var ok,ob:RealType;
begin
ok:=(x1.Y-x2.y)/(X1.X-x2.x);
ob:=(x2.Y*x1.x-x1.Y*x2.x)/(x1.X-x2.x);
Result:=(abs(TestPoint.Y-(ok*TestPoint.X+ob))<1e-13) and

((testpoint.x>=min(x1.x,x2.x)) and
(testpoint.x<=max(x2.x,x1.x))) and ((testpoint.y>=min(x1.y,x2.Y)) and
(testpoint.y<=max(x2.y,x1.Y)));
end;

function Orientation(Dot1,bDot,eDot,cDot:Dot):Integer;
var oa,ob,oc:Dot;
    sa,sb:RealType;
begin

  oa.x:=eDot.x-bDot.X;
  oa.y:=eDot.y-bDot.y;
  ob.x:=Dot1.x-bDot.X;
  ob.y:=Dot1.y-bDot.y;
  oc.X:=cDot.X-bDot.X;
  oc.y:=cDot.y-bDot.y;

//a. x * b.y - b.x * a.y;
sa:=oa.x*ob.y-ob.x*oa.y;
sb:=oa.x*oc.y-oc.x*oa.y;
Result:=0;// otherwise
if (sa*sb>=0) then Result:=-1; //left

end;

function InTri(TestPoint,Tria,Trib,Tric:Dot):boolean;
begin

Result:=(Orientation(TestPoint,Tria,Trib,Tric)=-1) and
 (Orientation(TestPoint,Trib,Tric,Tria)=-1) and
 (Orientation(TestPoint,Tric,Tria,Trib)=-1) or
  OnLine(TestPoint,Tria,Trib) or OnLine(TestPoint,Tria,Tric) or
  OnLine(TestPoint,TriB,TriC)
end;

GetDerivative - возвращает константу.
Тестировал на треугольнике
(0,54030230587, 0,84147098481)
(0,4782535689, 0,74549535193)
(0,48028810795, 0,87711078739)
значения коэффициентов полинома
aa=144,40721255
bb=-152,03135368
cc=-256,00702271
dd=134,64505084
ee=39,979898071
ff=113,36503213

[VLarin]

Geckelberryfinn писал(а):

Это сделано для универсальности. Чтобы потом можно было просто увеличивать степень полинома, не вычисляя заново аналитически значение интеграла.

Тоже не по теме:))
Функции формы выше 6 не используются в программах МКЭ ввиду сложности.
А для первых 6 порядков получить интегралы по КЭ для составления матрицы и по границам КЭ для учета ГУ 2 и 3 рода трудности нету, один раз вычислили, записали в код и забыли:)

Квадратурные формулы Гаусса и результаты для Вашего полинома:
Вычисление по 3 точкам:
1. Координаты площади L1 = (0.5, 0, 0), L2 = (0, 0.5, 0), L3 = (0, 0, 0.5), вес 1/3 каждый
A = 1.2883e-009
2. Координаты площади L1 = (2/3, 1/6, 1/6), L2 = (1/6, 2/3, 1/6), L3 = (1/6, 1/6, 2/3), вес 1/3
A = 1.2883e-009

[незваный гость]

1) Вы неправильно (ошибаетесь на 1) считаете количество повторений цикла в основной функции.

2) Я не вижу нестабильности. То, что наблюдаю я — стремление интеграла к 0 (по крайней мере, в пределах точности вычисления). Естественно, что если интеграл равен 0, то результат будет колебаться вокруг ответа в + и -. Однако амплитуда этих колебаний уменьшается у меня примерно обратно пропорционально количеству итераций.

3) Проверьте качество датчика случайных чисел. Поскольку Вы всегда генерируете точки парами, для Вас это очень существенно. Один из простых способов быстро улучшить датчик — алгорифм Маклорена-Марсальи (см. Кнут, т.2)

4) Принадлежность треугольнику проще проверять по равенству сумм удвоенных площадей (удвоенная площадь суть абсолютная величина векторного произведения): если $D$ лежит в треугольнике $\triangle ABC$ или на его границе, то $S_{ABC} = S_{ABD} + S_{BCD} + S_{CAD}$.

[VLarin]

Тоже бросилась в глаза слишком сложная проверка принадлежности треугольнику.
Зачем вычислять 3 площади и еще проверять отдельно принадлежность границам?:)

Проверка равенства трех площадей сравнением с площадью S_{ABC} - здесь нужно быть аккуратным, т.к. используются вещественные числа.
Самое простое и быстрое на мой взгляд - это вычислить последовательно 3 векторных произведения (площади) и проверить, есть ли хоть одно отрицательное:

if ((x - x2)*(y1 - y2) - (y - y2)*(x1 - x2) < 0)
   return 0;
if ((x - x3)*(y2 - y3) - (y - y3)*(x2 - x3) < 0)
   return 0;
if ((x - x1)*(y3 - y1) - (y - y1)*(x3 - x1) < 0)
   return 0;
return -1;

[Brukvalub]

незваный гость писал(а):

4) Принадлежность треугольнику проще проверять по равенству сумм удвоенных площадей (удвоенная площадь суть абсолютная величина векторного произведения): если $D$ лежит в треугольнике $\triangle ABC$ или на его границе, то $S_{ABC} = S_{ABD} + S_{BCD} + S_{CAD}$.

Треугольник лежит целиком по одну сторону от каждой из своих сторон, то есть является пересечением 3 полуплоскостей, так не проще ли выписать ур-ния сторон и проверять выполнение системы трех нестрогих линейных неравенств?:shock:

[незваный гость]

Это, конечно, верно, вопрос лишь, по какую (ответ: по ту же, что и оставшаяся вершина). Если Вы попытаетесь сделать это всё алгорифмически, и сравните объем вычислений (плюс, чем меньше if'ов, тем лучше), Вы увидите, что метод площадей не так уж плох. Он, помимо всего прочего, позволяет обойтись вообще без делений.

Добавлено спустя 12 минут 47 секунд:

VLarin писал(а):

Проверка равенства трех площадей сравнением с площадью S_{ABC} - здесь нужно быть аккуратным, т.к. используются вещественные числа.

Разумеется.

VLarin писал(а):

Самое простое и быстрое на мой взгляд - это вычислить последовательно 3 векторных произведения (площади) и проверить, есть ли хоть одно отрицательное:

Чего-то у меня не получилось . Я рассматривал треугольник (0,0), (4,0), (0,3), и точку (1,1). То есть, если я рассматривал, как написал, то точка внутри. А если (0,0), (0,3), (4,0) — то нет.

Кроме того, при этом подходе гораздо больше if'ов. Командопровод просто визжит от негодования.

Добавлено спустя 3 минуты 47 секунд:

я писал(а):

То есть, в вычислении у Вас может участвовать заметно меньше точек, чем Вы задаете.

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

Добавлено спустя 6 минут 47 секунд:

P.S. VLarin, Brukvalub:
Если действительно делать Монте-Карло (то есть, проверку большого количества точек к одному и тому же треугольнику), то, вероятно, имеет смысл пойти путем предвычисления (в той или иной форме) параметров треугольника с целью оптимизации вычисления принадлежности. С этой точки зрения ваши подходы эквиваленты (на мой взгляд), а подход с площадями — хуже.

[Someone]

незваный гость писал(а):

я писал(а):

То есть, в вычислении у Вас может участвовать заметно меньше точек, чем Вы задаете.

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

Чтобы генерировать поменьше лишних точек, можно воспользоваться барицентрическими координатами. Пусть вершины треугольника будут $M_1(x_1,y_1)$, $M_2(x_2,y_2)$, $M_3(x_3,y_3)$. Чтобы получить точку в треугольнике, генерируем два псевдослучайных числа $p$ и $q$, равномерно распределённых на $[0;1)$. Если $p+q\leqslant1$, вычисляем координаты точки по формулам $x=px_1+qx_2+(1-p-q)x_3$ и $y=py_1+qy_2+(1-p-q)y_3$. Если же $p+q>1$, то заменяем $p$ на $1-p$, $q$ на $1-q$ и вычисляем координаты точки по тем же формулам. И не надо ничего проверять.

[незваный гость]

Я чего-то не могу понять с ходу, будет ли работать такой подход: $M' = M_1 p + M_2 (1-p)$, $m = M' q + M_3 (1-q)$. Вроде он еще проще…

[Zai]

Интегрировать по треугольнику проще используя формулы из Зенкевича( границы прямолинейные и система координат расположена в центре тяжести).
(Зенкевич О. Метод конечных элементов в технике. -М.: Мир, 1975. - 318 с)

Приложение 3 стр. 229

[Geckelberryfinn]

VLarin
У меня треугольник билинейный. Поэтому там координаты и коэффициенты чуть-чуть другие.
Определение принадлежности точки треугольнику по площадям я тоже использовал (первый алгоритм, который пришел мне в голову). Однако, из-за "плавающей" арифметики он иногда игнорил вершины. А в вершинах базисные функции равны 1 или 0 и не хотелось терять эти значения. Приходилось проверять равенство площадей не строго а с неким допуском.... Что не придавало изящности алгоритму.
P.S. Однако, всем спасибо за советы и замечания. Мне они пригодятся. У меня большого опыта в МКЭ нету, поэтому важны любые заметки на эту тему.

[VLarin]

Geckelberryfinn писал(а):

VLarin
У меня треугольник билинейный. Поэтому там координаты и коэффициенты чуть-чуть другие.

Это в каком смысле чуть-чуть другие?:)))
У Вас я как понял обычный треугольный КЭ 2 порядка. Возникает вопрос зачем Вы записали полином в явном виде, а не использовали функции формы:)

Geckelberryfinn писал(а):

Определение принадлежности точки треугольнику по площадям я тоже использовал (первый алгоритм, который пришел мне в голову). Однако, из-за "плавающей" арифметики он иногда игнорил вершины. А в вершинах базисные функции равны 1 или 0 и не хотелось терять эти значения. Приходилось проверять равенство площадей не строго а с неким допуском.... Что не придавало изящности алгоритму.

Я про это и говорил, что с проверкой по сумме площадей здесь нужно быть аккуратным.
Еще раз повторюсь, что на мой взгляд здесь самое лучшее - это обход треугольника против часовой и проверка по векторным произведениям. Если нужно учесть допуск, то проверяем не просто на отрицательное значение, а с некоторой delta <= 1e-6.
Таким образом, если произведение AB x AX < -1e-6*|AB| - то точка X лежит ниже стороны AB, не входит в треугольник ABC, и расстояние до нее больше, чем допуск.
Если AB x AX >= -1e-6*|AB| и AB x AX <= 0 - точка с некоторым допуском лежит на грани AB.
Опыт показывает, что для double здесь хорошо работает допуск 1е-6 - 1е-8.

Более того в такой проверке меньше вычислений, чем в вычислении площадей.
Если одна из 3 проверок не прошла, остальные считать не нужно:)

Таким образом, в Вашем случае функция InTri будет содержать -

if ((x - x2)*(y1 - y2) - (y - y2)*(x1 - x2) < 1e-6*L12) 
   return 0; 
if ((x - x3)*(y2 - y3) - (y - y3)*(x2 - x3) < 1e-6*L23) 
   return 0; 
if ((x - x1)*(y3 - y1) - (y - y1)*(x3 - x1) < 1e-6*L31) 
   return 0; 
return -1; 

Geckelberryfinn писал(а):

P.S. Однако, всем спасибо за советы и замечания. Мне они пригодятся. У меня большого опыта в МКЭ нету, поэтому важны любые заметки на эту тему.

Спрашивайте еще:))