Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
5  
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? –  chris Oct 14 '12 at 19:47
    
@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? –  Seyhmus Güngören Oct 14 '12 at 20:01
    
well it's a progress; its now executable. –  chris Oct 14 '12 at 20:08
    
@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. –  Seyhmus Güngören Oct 14 '12 at 20:15
    
@SeyhmusGüngören What did FindRoot say? Was there any error message? This looks like a rather unfriendly equation. –  drN Oct 14 '12 at 20:30
show 4 more comments

1 Answer

up vote 7 down vote accepted

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

share|improve this answer
    
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)$ –  Seyhmus Güngören Oct 14 '12 at 21:05
    
There's also a root at a = b = 0.000622898. –  wxffles Oct 14 '12 at 21:07
    
@SeyhmusGüngören more to the point it seems all solutions have a==b –  chris Oct 14 '12 at 21:08
    
@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. –  Seyhmus Güngören Oct 14 '12 at 21:12
    
@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. –  Seyhmus Güngören Oct 14 '12 at 21:13
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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