# dealing with discontinuity at integration

My problem is about the solution of $2$ non-linear equations which have a single parameter. To be more precise, For given two densities; for example:

f0[y_]: = 1/(E^((1 + y)^2/2) Sqrt[2*Pi])


and

f1[y_]:= 1/(E^((-1 + y)^2/2) Sqrt[2*Pi])


define

l[y_] := f1[y]/f0[y] == E^(-(1/2) (-1 + y)^2 + 1/2 (1 + y)^2)

ll[y_] := Log[y]/2


then

h1[cu_] := Integrate[f0[y], {y, -∞, ll[cu]}] + 1/cu Integrate[f1[y], {y, ll[cu], ∞}]


and

h2[cl_] := cl*Integrate[f0[y], {y, -∞, ll[cl]}] + Integrate[f1[y], {y, ll[cl], ∞}]


have unique solutions which can be obtained via

FindRoot[h1[cu] == 1/(1 - 0.001), {cu, 0.1}]


and

FindRoot[h2[cl] == 1/(1 - 0.001), {cl, 0.1}]


the number $0.001$ above can be any number between $0$ and $1$. For the example above mathematica gives me the following results:

$${cu -> 68.9457}\quad {cl -> 0.0145042}$$

My Question: Now I would like to replace $f_0$ with $g_0$ and $f_1$ by $g_1$ for

g0[y_] :=
1/(4 (1/2 + m00/l00)^2) f0[y] UnitStep[-y + yll] +
(l00*Sqrt[f0[y]] + l11 Sqrt[f1[y]])^2/
(l00 + l11 + 2 (m00 + m11 - 1))^2 (UnitStep[y - yll] - UnitStep[y - yuu]) +
1/(4 (1/2 + (m00 - 1)/l00)^2) f0[y] UnitStep[y - yuu]


and

g1[y_] :=
1/(4 (1/2 + (m11 - 1)/l11)^2)f1[y] UnitStep[-y + yll] +
(l00 Sqrt[f0[y]] + l11 Sqrt[f1[y]])^2/
(l00 + l11 + 2 (m00 + m11 - 1))^2 (UnitStep[y - yll] - UnitStep[y - yuu]) +
1/(4 (1/2 + m11/l11)^2) f1[y] UnitStep[y - yuu]


with parameters :

l11 = 5.924124553902465

m11 = 0.2351824501432836

yuu = 0.375085

yll = -0.375085

m00 = m11

l00 = l11


 lll[y_] :=
0.5 Log[0.47228603296849714 y] (UnitStep[y] - UnitStep[y - 1]) +
0.5 Log[2.117360942720707 y] UnitStep[y - 1]


Now I am searching solutions to the following equations as above

h3[cu_] := Integrate[g0[y], {y, -∞, lll[cu]}] + 1/cu Integrate[g1[y], {y, lll[cu], ∞}]

h4[cu_]:= cl Integrate[g0[y], {y, -∞, lll[cl]}] + Integrate[g1[y], {y, lll[cl], ∞}]


via

FindRoot[h3[cu] == 1/(1 - 0.001), {cu, 1}]


and

FindRoot[h4[cl] == 1/(1 - 0.001), {cl, 0.1}]


However i get an error for the equation related to $h_3$

FindRoot::njnum: The Jacobian is not a matrix of numbers at {cu} = {1.}. >>

For the eqaution related to $h_4$ I am able to get a result:

$${cl -> 0.0225872}$$

lll[y] is not continious at $y=1$, I am binding the problem to this. I am expecting to get a result like $1/cl =44.2729$.

Because in the first example I also have the relation $cu=1/cl$.

How can I deal with the singularity at $y=0$ at the integration?

EDIT: Here is the gaphic of $h3$ my estimation seems correct. I get $1.001$ at around $44$

-
Please don't enter Mathematica code in LaTex format. Please read the formatting guidelines and format your code accordingly. It is essentially impossible for any to work with your code when it is posted in the way you have done it here. –  m_goldberg Mar 12 '13 at 2:23
@m_goldberg Yes I am aware of this. Please dont think that I did it because I like it or it would seem better or it was of much fun. I tried to get the codes in the inputform. I am using //InputForm command and right click copy as "inputform". With both of them I couldnt get something useful. It was starting with InputForm.... Thats the reason I posted the ones which I couldnt post in inputform as Latex. Could you please help me about what to do? I can definitely send whatever I get via //InputForm to you. –  Seyhmus Güngören Mar 12 '13 at 9:44
FindRoot[h1[cu] == 1/(1 - 0.001), {cu, 0.1}] \\ InputForm Syntax::tsntxi: "\!(* StyleBox[\"\\\"\\\\\\\"\\\"\", \"MT\"])\!(* StyleBox[ RowBox[{ RowBox[{\"FindRoot\", \"[\", RowBox[{ RowBox[{ RowBox[{\"h1\", \"[\", \"cu\", \"]\"}], \"==\", RowBox[{\"1\", \"/\", RowBox[{\"(\", RowBox[{\"1\", \"-\", \"0.001\"}], \")\"}]}]}], \",\", RowBox[{\"{\", RowBox[{\"cu\", \",\", \"0.1\"}], \"}\"}]}], \"]\"}], \" \", \"InputForm\"}], \"MT\"])\!(* StyleBox[\"\"\", \"MT\"])" is incomplete; more input is needed. –  Seyhmus Güngören Mar 12 '13 at 9:48
For the integral equations I get : Syntax::tsntxi: "g1[y][DifferentialD]y InputForm" is incomplete; more input is needed. Syntax::sntxi: Incomplete expression; more input is needed. –  Seyhmus Güngören Mar 12 '13 at 9:50
for the substitutions you need to do something like f-> Function[x,Sin[x]]  –  chris Mar 12 '13 at 9:55