Найти 1-е пересечение 2 PDF-функций с Mathematica - PullRequest
1 голос
/ 12 декабря 2011

В Mathematica 8.0.1.0 я использовал FindRoot[] для определения пересечения двух 2-х pdf-функций.

Но если у меня есть pdf-функции, которые пересекаются более чем в одной точке, и у меня есть верхняяпредел диапазона оси x за пределами второго пересечения, FindRoot[] возвращает только второе пересечение.

pdf1 = 1/x 0.5795367855565214` (E^(
  11.170058830053032` (-1.525439351903338` - Log[x]))
   Erfc[1.6962452696714152` (-0.5548887795964352` - Log[x])] + 
 E^(1.2932713057519` (2.60836043407439` + Log[x]))
   Erfc[1.6962452696714152` (2.720730943938539` + Log[x])]);

pdf2 = 1/x 0.4648445097126269` (E^(
  5.17560914275408` (-2.5500941338198615` - Log[x]))
   Erfc[1.7747318880142482` (-2.139288893723375` - Log[x])] + 
 E^(1.1332542415053757` (3.050849516581922` + Log[x]))
   Erfc[1.7747318880142482` (3.1407996592474956` + Log[x])]);

Plot[{pdf1, pdf2}, {x, 0, 0.5}, PlotRange -> All]   (* Shows 1st intersection *)
Plot[{pdf1, pdf2}, {x, 0.4, 0.5}, PlotRange -> All] (* Shows 2nd intersection *)

{x /. FindRoot[pdf1 == pdf2, {x, 0.00001, 0.5}],
x /. FindRoot[pdf1 == pdf2, {x, 0.00001, 0.4}]}

enter image description here enter image description here

Приведенные выше графики показывают проблему.При построении они пересекаются в двух точках:

{0.464719, 0.0452777}

соответственно.

Поскольку я не могу знать заранее, будет ли у меня второе пересечение, и я не знаю, где оно можетесли я упал на ось х, может кто-нибудь предложить способ, чтобы FindRoot[] возвращал только первое пересечение, а не второе?

Если нет, может кто-нибудь предложить другой способ сделать это?

Ответы [ 2 ]

4 голосов
/ 12 декабря 2011

С FindRoot[] вы можете получить только один корень для данной начальной точки.Перебор различных опций обременителен, и вы можете даже не получить желаемый результат для определенных крайних случаев, если не найдете правильный выбор начальной точки.

В этом случае что-то вроде NSolve или Reduce можетбыть лучшим вариантом.Если вы знаете, что ваши выражения затухают, используя разумную верхнюю границу для возможных значений x, вы можете использовать следующее, что довольно быстро и даст вам все корни.

NSolve[{pdf1 == pdf2, 0 < x < 1}, x] // Timing
Out[1]= {0.073495, {{x -> 0.0452777}, {x -> 0.464719}}}
1 голос
/ 12 декабря 2011

Как насчет следующего:

Сначала вы должны найти все корни за один шаг.Я делаю это с

roots=Reduce[pdf1==pdf2&&0.000001<x<0.5,x]

И тогда вы можете взять минимум (первый корень на оси x в вашем особом случае).

rootMin=Min[N[x/.{ToRules[roots]}]]
Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...