Qfyilyi

I am looking for some tips or guidance as to what machinery in mathematica can help me get at this problem numerically. I am looking for fixed points of a mapping, but the objects in question are themselves functions. Hence I am looking for a fixed point function.

The setup (simplified version):
suppose we restrict our search to continuous functions $f: [0,1]rightarrow [0,1]$. $p$ is a known parameter. I am looking for a fixed point (function) such that, for all $xin [0,1]$, $f(x)$ solves

$$f(x) = frac{x^p}{x^p + int_0^1 f(x) x^p , dx }.$$

It's not as simple as finding lots of fixed points for each $x$ in isolation, as the value of the expression at a single $x$ depends on the entire function $f$. Any help to try and solve this type of thing numerically would be much appreciated.

I believe this is @Lucas suggestions in the comment.

    ClearAll[p, c]

    p = 2;

    f[x_] := x^p/(x^p + c)    



c = c /. Quiet@FindRoot[NIntegrate[f[x] x^p, {x, 0, 1}] == c, {c, 1}]

0.227879

fixedPoints = NSolve[f[x] == x, x]

{{x -> 0.64873}, {x -> 0.35127}, {x -> 0.}}

 Plot[{f[x], x}, {x, 0, 1}, AspectRatio -> 1, Frame -> True, 

 GridLines -> {Flatten@Values@fixedPoints, 

   Flatten@Values@fixedPoints}]

This uses a discretization by piecewise-linear functions.

n = 1000;

x = Subdivide[0., 1., n - 1];

p = 2;

(* quadrature weights for trapezoidal rule *)

ω = ConstantArray[1./(n - 1), n];

ω[[1]] = ω[[n]] = 0.5/(n - 1);

Applying fixed point iteration; I use Dot and ω to compute the integral:

data = FixedPointList[

 f [Function] x^p/(x^p + (x^p f).ω), 

 ConstantArray[0.5, n]

 ];

Checking the $L^infty$ error:

Max[Abs[step[data[[-1]]] - data[[-1]]]]

2.22045*10^-16

Plotting the iterates:

ListLinePlot[

 data,

 PlotLegends -> Automatic

 ]

Just as an addition to @Okkes and @Ulrich's answer, we can also start with a symbolic solution of the integral for every p:

Integrate[x^p/(x^p + c) x^p, {x, 0, 1}, Assumptions -> p > 0 && c < -1]

Hypergeometric2F1[1, 2 + 1/p, 3 + 1/p, -(1/c)]/(c + 2 c p)

Here we had to cheat a bit with the assumption c < -1 to get the solution without condition, but we can see, that this is also a valid solution in the region c > 0 which is of interest to us (here for p==2):

Plot[(Hypergeometric2F1[1, 2 + 1/p, 3 + 1/p, -(1/c)]/(c + 2 c p)) /. p -> 2, {c, -2, 2}]

Next, we can use this solution to construct a function, which can numerically compute the constant c for arbitrary p (like @Okkes showed):

croot[p_?NumericQ] := Re[c /. FindRoot[

  Hypergeometric2F1[1, 2 + 1/p, 3 + 1/p, -(1/c)]/(c + 2 c p) == c,

  {c, 3/10}]

]

and then our solution will be

solution[p_] := x^p/(x^p + croot[p])

Now we can plot this for a range of p values:

Plot3D[solution[p], {x, 0, 1}, {p, 0.01, 4}, AxesLabel -> {"x", "p"}, MeshFunctions -> {#2 &}]

An interesting observation is, that f[x] seems to tend to a constant value as p tends to zero. With our symbolic solution from earlier we can determine that this value

Limit[Hypergeometric2F1[1, 2 + 1/p, 3 + 1/p, -(1/c)]/(c + 2 c p) - c, p -> 0, Direction -> -1] == 0

% && c > 0 // Solve

Limit[x^p/(x^p + c) /. First[%], p -> 0, Direction -> -1]

% // N

-((-1 + c + c^2)/(1 + c)) == 0

{{c -> 1/2 (-1 + Sqrt[5])}}

-(1/2) + Sqrt[5]/2

0.618034

is the golden ratio! :)

supplement

The list of fixpoint-functions can be obtained strictly numerical (variable p) using Nintegrate:

int[c_?NumericQ, p_?NumericQ] :=NIntegrate[x^(2 p)/(x^p + c), {x, 0, 1}, Method -> "LocalAdaptive"   ]    



f[Infinity] =Table[ x^ p /(x^p + c) /.NMinimize[{1, c == int[c, p]}, c][[2]] , {p, 1, 5}]

(*{x/(0.323829 + x), x^2/(0.227879 + x^2), x^3/(0.1782 + x^3)

, x^4/(0.147254 + x^4), x^5/(0.125923 + x^5)}*)



Plot[f[Infinity], {x, 0, 1}, PlotRange -> {0, 1}]

Another method

k = 20; int[0] = NIntegrate[x^p, {x, 0, 1}]; 

f[i_, x_] := x^p/(x^p + int[i])

 Table[

 Do[int[i] = NIntegrate[f[i - 1, x]*x^p, {x, 0, 1}]; kp = i; 

  If[Abs[int[i] - int[i - 1]] > 10^-6, Continue, Break], {i, 1, 

   k}]; x^p/(x^p + int[kp]), {p, 2, 7}]



(* {x^2/(0.227879 + x^2), x^3/(0.1782 + x^3), x^4/(

 0.147254 + x^4), x^5/(0.125923 + x^5), x^6/(0.110239 + x^6), x^7/(

 0.0981784 + x^7)}*)

For p=7

{ListPlot[Table[{i, int[i]}, {i, 0, kp}], PlotRange -> All], 

 Plot[Evaluate[Table[f[i, x], {i, 1, kp}]], {x, 0, 1}], 

 Plot[f[kp, x], {x, 0, 1}]}