2
$\begingroup$

I have the following equations and I would like to solve them in Mathematica anyway. Numerical results are very much welcomed. However, when I use the command NSolve, I get the following error message:

The expression $b^{\frac{1}{2}\log (ab)}$ involves unknowns in more than one argument, so inverse functions cannot be used.

Here are the equations that I want to solve:

NSolve[
 {(Integrate[(E^(-x^2/2 + Log[a]^2/Log[a*b])*(E^(-(-1 + x)^2/2 + x^2/2))^(Log[a]/Log[a*b])*Log[E^(Log[a]^2/Log[a*b])*(E^(-(-1 + x)^2/2 + x^2/2))^(Log[a]/Log[a*b])])/
  Sqrt[2*Pi], {x, 1/2 + Log[a^(-1)], 1/2 + Log[b]}] + 
(a*Erfc[(1/2 + Log[b])/Sqrt[2]]*Log[a])/2)/((1 + Erf[(1/2 + Log[a^(-1)])/Sqrt[2]])/2 + 
(E^((Log[a]^2*(1 + 2*Log[a*b]))/(2*Log[a*b]^2))*(-Erf[(-2*Log[a] + (1 + 2*Log[a^(-1)])*Log[a*b])/(2*Sqrt[2]*Log[a*b])] + 
 Erf[(-2*Log[a] + (1 + 2*Log[b])*Log[a*b])/(2*Sqrt[2]*Log[a*b])]))/(2*a^(1/(2*Log[a*b]))) + (a*Erfc[(1/2 + Log[b])/Sqrt[2]])/2) - 
 Log[(1 + Erf[(1/2 + Log[a^(-1)])/Sqrt[2]])/2 + (E^((Log[a]^2*(1 + 2*Log[a*b]))/(2*Log[a*b]^2))*
  (-Erf[(-2*Log[a] + (1 + 2*Log[a^(-1)])*Log[a*b])/(2*Sqrt[2]*Log[a*b])] + Erf[(-2*Log[a] + (1 + 2*Log[b])*Log[a*b])/(2*Sqrt[2]*Log[a*b])]))/
  (2*a^(1/(2*Log[a*b]))) + (a*Erfc[(1/2 + Log[b])/Sqrt[2]])/2] == 0.1, 
(Integrate[(E^(-(-1 + x)^2/2 + Log[b]^2/Log[a*b])*Log[E^(Log[b]^2/Log[a*b])/
 (E^(-(-1 + x)^2/2 + x^2/2))^(Log[b]/Log[a*b])])/
((E^(-(-1 + x)^2/2 + x^2/2))^(Log[b]/Log[a*b])*Sqrt[2*Pi]), {x, 1/2 + Log[a^(-1)], 1/2 + Log[b]}] + 
 (b*(1 + Erf[(-1/2 + Log[a^(-1)])/Sqrt[2]])*Log[b])/2)/
   ((b*(1 + Erf[(-1/2 + Log[a^(-1)])/Sqrt[2]]))/2 + (E^((Log[b]^2*(1 + 2*Log[a*b]))/
 (2*Log[a*b]^2))*(Erf[(-1 + 2*Log[b]*(1 + Log[a*b]^(-1)))/(2*Sqrt[2])] - 
  Erf[(-1 + 2*Log[a^(-1)] + (2*Log[b])/Log[a*b])/(2*Sqrt[2])]))/
   (2*b^(1/(2*Log[a*b]))) + Erfc[(-1/2 + Log[b])/Sqrt[2]]/2) - 
  Log[(b*(1 + Erf[(-1/2 + Log[a^(-1)])/Sqrt[2]]))/2 + 
   (E^((Log[b]^2*(1 + 2*Log[a*b]))/(2*Log[a*b]^2))*(Erf[(-1 + 2*Log[b]*(1 + Log[a*b]^(-1)))/(2*Sqrt[2])] - 
    Erf[(-1 + 2*Log[a^(-1)] + (2*Log[b])/Log[a*b])/(2*Sqrt[2])]))/(2*b^(1/(2*Log[a*b]))) + Erfc[(-1/2 + Log[b])/Sqrt[2]]/2] == 0.1}, 
{a, b}]

f0[x_] := 1/Sqrt[2*Pi*1^2]*E^(-(x + 0)^2/(2*1^2))
f1[x_] := 1/Sqrt[2*Pi*1^2]*E^(-(x - 1)^2/(2*1^2))

yl = -1.070751889428462
yu = 2.0707518894284505

a = 1/l[yl]
b = l[yu]

l[x] = f1[x]/f0[x]

Solve[-(1/2) (-1 + x)^2 + x^2/2 - Log[y] == 0, x]   

--> here I must find the inverse function of $l$ but I couldn't do that either, so then I did it half-manually, which is:

ll[x_] := 1/2 + Log[x]

NSolve[{(Integrate[a*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])]*Log[a], {x, ll[b], Infinity}] + Integrate[E^(Log[a]^2/Log[a*b])*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])]*
        l[x]^(Log[a]/Log[a*b])*Log[E^(Log[a]^2/Log[a*b])*l[x]^(Log[a]/Log[a*b])], {x, ll[a^(-1)], ll[b]}])/
     (Integrate[InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])], {x, -Infinity, ll[a^(-1)]}] + Integrate[a*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])], {x, ll[b], Infinity}] + 
      Integrate[E^(Log[a]^2/Log[a*b])*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])]*l[x]^(Log[a]/Log[a*b]), {x, ll[a^(-1)], ll[b]}]) - 
    Log[Integrate[InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])], {x, -Infinity, ll[a^(-1)]}] + Integrate[a*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])], {x, ll[b], Infinity}] + 
      Integrate[E^(Log[a]^2/Log[a*b])*InputForm[1/(E^(x^2/2)*Sqrt[2*Pi])]*l[x]^(Log[a]/Log[a*b]), {x, ll[a^(-1)], ll[b]}]] == 0.1, 
  (Integrate[b*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])]*Log[b], {x, -Infinity, ll[a^(-1)]}] + 
      Integrate[(E^(Log[b]^2/Log[a*b])*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])]*Log[E^(Log[b]^2/Log[a*b])/l[x]^(Log[b]/Log[a*b])])/l[x]^(Log[b]/Log[a*b]), 
       {x, ll[a^(-1)], ll[b]}])/(Integrate[InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])], {x, ll[b], Infinity}] + 
      Integrate[b*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])], {x, -Infinity, ll[a^(-1)]}] + 
      Integrate[(E^(Log[b]^2/Log[a*b])*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])])/l[x]^(Log[b]/Log[a*b]), {x, ll[a^(-1)], ll[b]}]) - 
    Log[Integrate[InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])], {x, ll[b], Infinity}] + Integrate[b*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])], 
       {x, -Infinity, ll[a^(-1)]}] + Integrate[(E^(Log[b]^2/Log[a*b])*InputForm[1/(E^((-1 + x)^2/2)*Sqrt[2*Pi])])/l[x]^(Log[b]/Log[a*b]), {x, ll[a^(-1)], ll[b]}]] == 
   0.1}, {a, b}]

Please let me know if there is still something missing.

MORE NOTES:

1- $y_u>y_l$ should be satisfied

2- $a,b>0$ and $ab>1$

3- $l[x]$ monotone increasing. This means one can choose some densities $f0[x_]$ and $f1[x_]$ such that $l$ is monotone increasing. Because we need its inverse!

According to these notes, is the solution unique? or are there still more roots if I change $f0[x_]$ and $f1[x_]$?

Is there anyway to get a solution? I also tried FindRoot and SolveAlways which gave no result at all. Thank you very much in advance.

$\endgroup$
9
  • 5
    $\begingroup$ Welcome to mathematica.stackexchange. Could you please type //InputForm at the end of the cell you copied, evaluate it and replace what you wrote with the corresponding result? $\endgroup$
    – chris
    Commented Oct 14, 2012 at 19:47
  • $\begingroup$ @chris thanks alot for welcoming. Let me try it. Should I do it in mathematica? I guess yes. I am doing it. How about now? $\endgroup$ Commented Oct 14, 2012 at 20:01
  • $\begingroup$ well it's a progress; its now executable. $\endgroup$
    – chris
    Commented Oct 14, 2012 at 20:08
  • $\begingroup$ @chris great! Mathematica, I just downloaded yesterday and learned a bit. It seems a very strong tool to deal with equations. I would never ever use MATLAB for this reason. $\endgroup$ Commented Oct 14, 2012 at 20:15
  • $\begingroup$ @SeyhmusGüngören What did FindRoot say? Was there any error message? This looks like a rather unfriendly equation. $\endgroup$
    – dearN
    Commented Oct 14, 2012 at 20:30

1 Answer 1

7
$\begingroup$

It seems this problem falls into the ?NumericQ category

If you type

 eqn = <your input>[[1]];
 F1[a_, b_] = eqn[[1, 1]] /. Integrate -> Int;
 F2[a_, b_] = eqn[[2, 1]] /. Integrate -> Int;
 F11[a_?NumericQ, b_?NumericQ] := F1[a, b] /. Int -> NIntegrate;
 F22[a_?NumericQ, b_?NumericQ] := F2[a, b] /. Int -> NIntegrate; 

 ContourPlot[{F11[a, b] == 0.1, F22[a, b] == 0.1}, {a, 0.1, 0.9}, {b, 0.1, 0.9}]

you get this contour plot.

Mathematica graphics

which suggests a solution exists near {0.4,0.4}. The command

 FindRoot[{F11[X, Y] == 0.1, F22[X, Y] == 0.1}, {{X, 0.4}, {Y, 0.4}}]

finds it

(* {X->0.38513,Y->0.38513} *)

If you take more seriously the hint that all solutions should have a==b

Plot[{0.1 + x*0, F11[x, x]}, {x, 0.01, 15}]

returns

Mathematica graphics

$\endgroup$
5
  • $\begingroup$ Woow it is great. Meanwhile @drN and chris, I tried findroot and here is the result : {$a$ -> $4.810263620768754$, $b$ -> $4.810263620768698$}. Very interesting because it seems the question has more solutions! I checked $0.38513$ and $4.810263620768754$. For $0.38513$ we have $y_l=1.45$ and $y_u=-0.45$. For $4.810263620768754$ we have $y_l=-1.070751889428462$ and $y_u=2.0707518894284505$. In my problem $y_u>y_l$ should be satisfied. Please let me also introduce you what $y_l$ and $y_u$ are. $a = 1/l[yl]$ and $b = l[yu]$ and $l[x]=E^(-(-1 + x)^2/2 + x^2/2)$ $\endgroup$ Commented Oct 14, 2012 at 21:05
  • $\begingroup$ There's also a root at a = b = 0.000622898. $\endgroup$
    – wxffles
    Commented Oct 14, 2012 at 21:07
  • $\begingroup$ @SeyhmusGüngören more to the point it seems all solutions have a==b $\endgroup$
    – chris
    Commented Oct 14, 2012 at 21:08
  • $\begingroup$ @wxffles for that case $y_l=7.88113$ and $y_u=-6.88113$. Again $y_l>y_u$. I am interested in if a solution is unique, with the constraint $y_u>y_l$. I think I could better post the equations. $\endgroup$ Commented Oct 14, 2012 at 21:12
  • $\begingroup$ @chris it is due to the symmetry of the densities. It is an expected result. Please give me some time and I post some more equations. Thanks again. $\endgroup$ Commented Oct 14, 2012 at 21:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.