# Composition of implicit functions

Let me explain what i have been doing so far. So I want to find an numerical solution or at least a grafic solution to this equation: R[x_] := (x/(x + w))*(1 - G0[1 - (x + w)*finf[x]]). Where:

w := 0.1

G0[x_] := Sum[P[k]*x^k, {k, 2, 250}]

P[x_] := (x^(-l))*Exp[-x/c]

l := 2

c := 25

finf[x]: = (x/(x + w))*(1 - G1[1 - (x + w)*finf[x]])

kmed = Sum[P[k]*k, {k, 2, 250}]*1.

G1[x_] := Sum[(k*P[k]*x^(k - 1))/(kmed*1.), {k, 2, 250}]

It´s important to note that I didn´t use the function finf as I just wrote it just for you to know the form of this implicit function. So, I tried Solve and NSolve finf and wouldn´t work so i came up with the idea of using a table with its values so i can at least get a discreate solution to R which lead to: d = Table[y = (x/(x + w))*(1 - G1[1 - (x + w)*y]), {x, 0, 1, 0.001}] At the same time I used ContourPlot to get this function but when I compare the results I get:

If I try using a smaller step, like 0.0001 or 0.00001 it gets a little better but still kind of bad when it shouldn´t and that part of the grafic is kind of important because it is when the solution stops being zero so I kind of need to make it better. Any advices or another way of doing this is very much apreciated (new in mathematica doing the best I can)

• Are R and finf missing parentheses? – Michael E2 Jan 5 at 18:02
• Yes, sorry, missed it when copying, already edited – Cbruno Jan 5 at 18:08

Please have a look at the distinction between immediate and delayed assignments.

First, define the constants:

w = 1/10;
c = 25;
l = 2;


I've modified your definitions to sum to infinity instead of to 250, assuming that that's what you had in mind:

P[x_] = (x^(-l))*Exp[-x/c];
G0[x_] = Sum[P[k]*x^k, {k, 2, ∞}] // FullSimplify
(*    -(x/E^(1/25)) + PolyLog[2, x/E^(1/25)]    *)

kmed = Sum[P[k]*k, {k, 2, ∞}] // FullSimplify
(*    -(1/E^(1/25)) - Log[1 - 1/E^(1/25)]    *)

G1[x_] = Sum[(k*P[k]*x^(k - 1))/(kmed), {k, 2, ∞}] // FullSimplify
(*    (x + E^(1/25) Log[1 - x/E^(1/25)])/(x + E^(1/25) x Log[1 - 1/E^(1/25)])    *)


Now we can define finf through numerical root-search:

finf[x_?NumericQ] := f /. NSolve[f == (x/(x + w))*(1 - G1[1 - (x + w)*f]), f, Reals]


Try it out: there are always two solutions,

finf[0.1]
(*    {-0.0249019, 0.}    *)

finf[0.2]
(*    {0., 0.257607}    *)

finf[0.3]
(*    {0., 0.43917}    *)


One of these solutions is always zero, because G1[1]==1. You may need extra code to isolate the nonzero solution.

Define the R function from here: it is also double-valued because finf has two solutions,

R[x_?NumericQ] := (x/(x + w))*(1 - G0[1 - (x + w)*finf[x]])


Try it out:

R[0.1]
(*    {0.236694, 0.242505}    *)

R[0.2]
(*    {0.32334, 0.413701}    *)

R[0.3]
(*    {0.363757, 0.54988}    *)


Make a plot of all solutions: you can see the two branches of solutions,

ListPlot[Join @@ Table[Thread@{x, R[x]}, {x, 0, 1, 1/1000}]]


One solution, coming from finf[x]==0, is (x/(x + w))*(1 - G0[1]); the other one is probably what you're after.