Numerical solution of two equations in two functions

I would like to obtain a solution to the following two equations, where g0[y] and g1[y] are the unknown functions:

-0.23094021403064358 + (-2 g0[y] + g1[y])/(1.043747093724331 g0[y] + 1.0437470937243307 g1[y])
- (1.043747093724331 (-g0[y]^2 + g0[y] g1[y]))/(1.043747093724331 g0[y] + 1.0437470937243307 g1[y])^2
- 0.9580865000847749 Log[E^(1/2 (1 + y)^2) Sqrt[2 \[Pi]] g0[y]] == 0

-0.23094021403064324 + (-g0[y] + 2 g1[y])/(1.043747093724331 g0[y] + 1.0437470937243307 g1[y])
- (1.0437470937243307 (-g0[y] g1[y] + g1[y]^2))/(1.043747093724331 g0[y] + 1.0437470937243307 g1[y])^2
- 0.9580865000847752 Log[E^(1/2 (-1 + y)^2) Sqrt[2 \[Pi]] g1[y]] == 0


I dont need an analytical solution, a numerical solution would suffice. I tried using Solve but I didn't manage to obtain any solutions.

• Consider using FindRoot. After obtaining roots for a particular value of y, use those roots as initial guesses for a nearby value of y, etc. – bbgodfrey Dec 1 '14 at 5:05

Just to elaborate the advice of bbgodfrey. I take it, you do not need to have that precise numbers. The 2-digit precision is enough? Then here are your equations:

  eq1 = -0.23094021403064358 + (-2 g0[y] +
g1[y])/(1.043747093724331 g0[y] +
1.0437470937243307 g1[y]) - (1.043747093724331 (-g0[y]^2 +
g0[y] g1[y]))/(1.043747093724331 g0[y] +
1.0437470937243307 g1[y])^2 -
0.9580865000847749 Log[E^(1/2 (1 + y)^2) Sqrt[2 \[Pi]] g0[y]] == 0

eq2 = -0.23094021403064324 + (-g0[y] +
2 g1[y])/(1.043747093724331 g0[y] +
1.0437470937243307 g1[
y]) - (1.0437470937243307 (-g0[y] g1[y] +
g1[y]^2))/(1.043747093724331 g0[y] +
1.0437470937243307 g1[y])^2 -
0.9580865000847752 Log[E^(1/2 (-1 + y)^2) Sqrt[2 \[Pi]] g1[y]] ==0


This will make them less ugly and rationalize them. The latter will help Mma to work faster:

 rnd[expr_, n_Integer] := expr /. x_Real :> Round[x, .1/10^(n - 1)];
eq1A = Rationalize[Map[rnd[#, 2] &, eq1, Infinity] /. f_[y] -> f]

(*   -(23/100) + (-2 g0 + g1)/((26 g0)/25 + (26 g1)/25) - (
26 (-g0^2 + g0 g1))/(25 ((26 g0)/25 + (26 g1)/25)^2) -
24/25 Log[E^(1/2 (1 + y)^2) g0 Sqrt[2 \[Pi]]] == 0    *)


and

 eq2A = eq1A =
Rationalize[Map[rnd[#, 2] &, eq2, Infinity] /. f_[y] -> f]

(* -(23/100) + (-g0 + 2 g1)/((26 g0)/25 + (26 g1)/25) - (
26 (-g0 g1 + g1^2))/(25 ((26 g0)/25 + (26 g1)/25)^2) -
24/25 Log[E^(1/2 (-1 + y)^2) g1 Sqrt[2 \[Pi]]] == 0     *)


These equations can be solved with the FindRoot for a sequence of yvalues and the solution is formed as a list with the structure {y,{g0,g1}}:

  lst = Table[{y, (FindRoot[{-(23/100) + (-2 g0 + g1)/((26 g0)/25 + (
26 g1)/25) - (26 (-g0^2 + g0 g1))/(
25 ((26 g0)/25 + (26 g1)/25)^2) -
24/25 Log[E^(1/2 (1 + y)^2) g0 Sqrt[2 \[Pi]]] ==
0, -(23/100) + (-g0 + 2 g1)/((26 g0)/25 + (26 g1)/25) - (
26 (-g0 g1 + g1^2))/(25 ((26 g0)/25 + (26 g1)/25)^2) -
24/25 Log[E^(1/2 (-1 + y)^2) g1 Sqrt[2 \[Pi]]] == 0}, {{g0,
1}, {g1, 2}}] /. {x_ -> a_, xx_ -> b_} -> {a, b}) //
Chop}, {y, 0, 2, 0.1}]


Evaluate it and you will get the list in question. Now you may do whatever you need with this list. For example, one may plot it:

 ListPlot[{(lst /. {x_, {y_, z_}} -> {x,
y}), (lst /. {x_, {y_, z_}} -> {x, z})},
PlotStyle -> {Blue, Red},
AxesLabel -> {Style["y", 16, Italic],
Row[{Style["\!$$\*SubscriptBox[\(g$$, $$0$$]\)", 16, Italic, Blue],
Style["  \!$$\*SubscriptBox[\(g$$, $$1$$]\)", 16, Italic,
Red]}]}]


This should look like the following:

Have fun!

• What is the range of the plot? Is it possible to get a figure for y in -6 to 6? – Seyhmus Güngören Dec 1 '14 at 15:46
• @Seyhmus Güngören The PlotRange is up to you. I did not fix it, but you can do it as you like. You should understand that I do not know, what are you after with this problem, so I just indicate the direction in which you could come. Concerning going in y from -6 to 6 just try, and have a look at the lst output. I remember that at y>2 the result became complex, That is, no real solution. However, I have arbitrarily chosen the initial guess: g0=1; g1=2 . It may be that with some other initial values you will obtain solutions. I recommend to ParametricPlot3D first your functions. – Alexei Boulbitch Dec 1 '14 at 15:55