Functional fixed points (ie fixed point of mapping from function space C[0,1] to itself)

Question

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 $x\in [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.

Answer 1

Just as an addition to @Okkes and @Ulrich's answer following the idea lined out by @LukasLang, 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! :)

Answer 2

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

Answer 3

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

@Thies's observation can be done analytically.

Set p=0. Then, $f(x) = \frac{1}{1+ c }$ where $c=\int_0^1 f(x) \, dx$

$c=\int_0^1 \frac{1}{1+ c } \, dx$

$c= \frac{1}{1+ c }x|_0^1 $

$c= \frac{1}{1+ c } $

$c^2+c-1= 0 $

$c= \frac{-1\mp\sqrt{5}}{2} $

Answer 4

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
 ]

Answer 5

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