1
$\begingroup$

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:

The soft one is from ContourPlot

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)

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

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}]]

enter image description here

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

$\endgroup$
0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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