# 2 non-linear equations with 2 unknowns which I cannot solve

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.

• 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? 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? Oct 14 '12 at 20:01
• well it's a progress; its now executable. 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. 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. Oct 14 '12 at 20:30

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.

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

• 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)$ Oct 14 '12 at 21:05
• There's also a root at a = b = 0.000622898. Oct 14 '12 at 21:07
• @SeyhmusGüngören more to the point it seems all solutions have a==b 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. 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. Oct 14 '12 at 21:13