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












7














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.










share|improve this question









New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
















  • 1




    I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
    – Thies Heidecke
    yesterday






  • 1




    - perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
    – Thies Heidecke
    yesterday








  • 5




    The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
    – Lukas Lang
    yesterday










  • Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
    – user434180
    yesterday












  • @LukasLang Great idea, this seems like the most simple and best approach!
    – Thies Heidecke
    yesterday
















7














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.










share|improve this question









New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.
















  • 1




    I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
    – Thies Heidecke
    yesterday






  • 1




    - perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
    – Thies Heidecke
    yesterday








  • 5




    The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
    – Lukas Lang
    yesterday










  • Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
    – user434180
    yesterday












  • @LukasLang Great idea, this seems like the most simple and best approach!
    – Thies Heidecke
    yesterday














7












7








7


3





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.










share|improve this question









New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.











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.







numerical-integration parametric-functions numerical-value fixed-points






share|improve this question









New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.











share|improve this question









New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









share|improve this question




share|improve this question








edited yesterday





















New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.









asked yesterday









user434180

362




362




New contributor




user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.





New contributor





user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.






user434180 is a new contributor to this site. Take care in asking for clarification, commenting, and answering.
Check out our Code of Conduct.








  • 1




    I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
    – Thies Heidecke
    yesterday






  • 1




    - perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
    – Thies Heidecke
    yesterday








  • 5




    The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
    – Lukas Lang
    yesterday










  • Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
    – user434180
    yesterday












  • @LukasLang Great idea, this seems like the most simple and best approach!
    – Thies Heidecke
    yesterday














  • 1




    I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
    – Thies Heidecke
    yesterday






  • 1




    - perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
    – Thies Heidecke
    yesterday








  • 5




    The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
    – Lukas Lang
    yesterday










  • Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
    – user434180
    yesterday












  • @LukasLang Great idea, this seems like the most simple and best approach!
    – Thies Heidecke
    yesterday








1




1




I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
– Thies Heidecke
yesterday




I can think of multiple ways to go about this: -discretize your function by representing it by a vector and solve the discretized problem as an approximation. Then this might result in a finite dimensional Eigenvalue problem which can be solved with Eigensystem or NDEigensystem. -use Interpolation as function representation and sample and reinterpolate after each iteration. - use a variational approach to find the fixed point, perhaps the VariationalMethods package can help with that.
– Thies Heidecke
yesterday




1




1




- perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
– Thies Heidecke
yesterday






- perhaps the problem can be stated as an ordinary differential equation and either be directly solved by DSolve, numerically by NDSolve or iteratively by a Picard iteration. I'm not sure if all of those methods can successfully tackle your problem, but all those alleys could be explored in Mathematica. Hope this gives you some ideas! When you try something and need further help, update your question with Mathematica code and specific questions, so that people can help you with the details.
– Thies Heidecke
yesterday






5




5




The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
– Lukas Lang
yesterday




The following should work: replace the integral with $c$. Solve for $f(x)$. Compute the integral as a function of $c$. Set the integral equal to $c$ and solve for $c$.
– Lukas Lang
yesterday












Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
– user434180
yesterday






Thank you very much for these suggestions. I have tried the ODE approach. I don't think that it can be transformed into an ODE as I don't see a way to remove the integral. The ratio means it can't be transformed into a Fredholm equation, for which code already exists to numerically solve. I will try discretizing in line with the great answer below. And I am a little lost by your suggestion, Lukas, although I will try work through it as well.
– user434180
yesterday














@LukasLang Great idea, this seems like the most simple and best approach!
– Thies Heidecke
yesterday




@LukasLang Great idea, this seems like the most simple and best approach!
– Thies Heidecke
yesterday










5 Answers
5






active

oldest

votes


















6














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


verifying that plot also works for c>0



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


Function plot for different values of p



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! :)






share|improve this answer































    5














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


    enter image description here






    share|improve this answer





















    • This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
      – user434180
      yesterday





















    5














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


    enter image description here



    @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{-1mpsqrt{5}}{2} $






    share|improve this answer























    • FixedPoint appears to be faster than NSolve for this.
      – Alan
      yesterday



















    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
    ]


    enter image description here






    share|improve this answer





























      0














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


      fig1






      share|improve this answer





















        Your Answer





        StackExchange.ifUsing("editor", function () {
        return StackExchange.using("mathjaxEditing", function () {
        StackExchange.MarkdownEditor.creationCallbacks.add(function (editor, postfix) {
        StackExchange.mathjaxEditing.prepareWmdForMathJax(editor, postfix, [["$", "$"], ["\\(","\\)"]]);
        });
        });
        }, "mathjax-editing");

        StackExchange.ready(function() {
        var channelOptions = {
        tags: "".split(" "),
        id: "387"
        };
        initTagRenderer("".split(" "), "".split(" "), channelOptions);

        StackExchange.using("externalEditor", function() {
        // Have to fire editor after snippets, if snippets enabled
        if (StackExchange.settings.snippets.snippetsEnabled) {
        StackExchange.using("snippets", function() {
        createEditor();
        });
        }
        else {
        createEditor();
        }
        });

        function createEditor() {
        StackExchange.prepareEditor({
        heartbeatType: 'answer',
        autoActivateHeartbeat: false,
        convertImagesToLinks: false,
        noModals: true,
        showLowRepImageUploadWarning: true,
        reputationToPostImages: null,
        bindNavPrevention: true,
        postfix: "",
        imageUploader: {
        brandingHtml: "Powered by u003ca class="icon-imgur-white" href="https://imgur.com/"u003eu003c/au003e",
        contentPolicyHtml: "User contributions licensed under u003ca href="https://creativecommons.org/licenses/by-sa/3.0/"u003ecc by-sa 3.0 with attribution requiredu003c/au003e u003ca href="https://stackoverflow.com/legal/content-policy"u003e(content policy)u003c/au003e",
        allowUrls: true
        },
        onDemand: true,
        discardSelector: ".discard-answer"
        ,immediatelyShowMarkdownHelp:true
        });


        }
        });






        user434180 is a new contributor. Be nice, and check out our Code of Conduct.










        draft saved

        draft discarded


















        StackExchange.ready(
        function () {
        StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f188770%2ffunctional-fixed-points-ie-fixed-point-of-mapping-from-function-space-c0-1-to%23new-answer', 'question_page');
        }
        );

        Post as a guest















        Required, but never shown

























        5 Answers
        5






        active

        oldest

        votes








        5 Answers
        5






        active

        oldest

        votes









        active

        oldest

        votes






        active

        oldest

        votes









        6














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


        verifying that plot also works for c>0



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


        Function plot for different values of p



        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! :)






        share|improve this answer




























          6














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


          verifying that plot also works for c>0



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


          Function plot for different values of p



          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! :)






          share|improve this answer


























            6












            6








            6






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


            verifying that plot also works for c>0



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


            Function plot for different values of p



            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! :)






            share|improve this answer














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


            verifying that plot also works for c>0



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


            Function plot for different values of p



            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! :)







            share|improve this answer














            share|improve this answer



            share|improve this answer








            edited 14 hours ago

























            answered 23 hours ago









            Thies Heidecke

            6,9462438




            6,9462438























                5














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


                enter image description here






                share|improve this answer





















                • This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                  – user434180
                  yesterday


















                5














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


                enter image description here






                share|improve this answer





















                • This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                  – user434180
                  yesterday
















                5












                5








                5






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


                enter image description here






                share|improve this answer












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


                enter image description here







                share|improve this answer












                share|improve this answer



                share|improve this answer










                answered yesterday









                Ulrich Neumann

                7,580516




                7,580516












                • This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                  – user434180
                  yesterday




















                • This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                  – user434180
                  yesterday


















                This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                – user434180
                yesterday






                This is a very helpful addition. I was about to comment on Henrik's answer that it requires the integral to be analytically solvable by mathematica. Although the one I posted in the question is solvable, the actual ones I care about are generally not. So this is very useful. Thanks!
                – user434180
                yesterday













                5














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


                enter image description here



                @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{-1mpsqrt{5}}{2} $






                share|improve this answer























                • FixedPoint appears to be faster than NSolve for this.
                  – Alan
                  yesterday
















                5














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


                enter image description here



                @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{-1mpsqrt{5}}{2} $






                share|improve this answer























                • FixedPoint appears to be faster than NSolve for this.
                  – Alan
                  yesterday














                5












                5








                5






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


                enter image description here



                @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{-1mpsqrt{5}}{2} $






                share|improve this answer














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


                enter image description here



                @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{-1mpsqrt{5}}{2} $







                share|improve this answer














                share|improve this answer



                share|improve this answer








                edited 11 hours ago

























                answered yesterday









                Okkes Dulgerci

                4,1151816




                4,1151816












                • FixedPoint appears to be faster than NSolve for this.
                  – Alan
                  yesterday


















                • FixedPoint appears to be faster than NSolve for this.
                  – Alan
                  yesterday
















                FixedPoint appears to be faster than NSolve for this.
                – Alan
                yesterday




                FixedPoint appears to be faster than NSolve for this.
                – Alan
                yesterday











                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
                ]


                enter image description here






                share|improve this 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
                  ]


                  enter image description here






                  share|improve this answer
























                    4












                    4








                    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
                    ]


                    enter image description here






                    share|improve this answer












                    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
                    ]


                    enter image description here







                    share|improve this answer












                    share|improve this answer



                    share|improve this answer










                    answered yesterday









                    Henrik Schumacher

                    49.1k467139




                    49.1k467139























                        0














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


                        fig1






                        share|improve this answer


























                          0














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


                          fig1






                          share|improve this answer
























                            0












                            0








                            0






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


                            fig1






                            share|improve this answer












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


                            fig1







                            share|improve this answer












                            share|improve this answer



                            share|improve this answer










                            answered 22 hours ago









                            Alex Trounev

                            6,1601419




                            6,1601419






















                                user434180 is a new contributor. Be nice, and check out our Code of Conduct.










                                draft saved

                                draft discarded


















                                user434180 is a new contributor. Be nice, and check out our Code of Conduct.













                                user434180 is a new contributor. Be nice, and check out our Code of Conduct.












                                user434180 is a new contributor. Be nice, and check out our Code of Conduct.
















                                Thanks for contributing an answer to Mathematica Stack Exchange!


                                • Please be sure to answer the question. Provide details and share your research!

                                But avoid



                                • Asking for help, clarification, or responding to other answers.

                                • Making statements based on opinion; back them up with references or personal experience.


                                Use MathJax to format equations. MathJax reference.


                                To learn more, see our tips on writing great answers.





                                Some of your past answers have not been well-received, and you're in danger of being blocked from answering.


                                Please pay close attention to the following guidance:


                                • Please be sure to answer the question. Provide details and share your research!

                                But avoid



                                • Asking for help, clarification, or responding to other answers.

                                • Making statements based on opinion; back them up with references or personal experience.


                                To learn more, see our tips on writing great answers.




                                draft saved


                                draft discarded














                                StackExchange.ready(
                                function () {
                                StackExchange.openid.initPostLogin('.new-post-login', 'https%3a%2f%2fmathematica.stackexchange.com%2fquestions%2f188770%2ffunctional-fixed-points-ie-fixed-point-of-mapping-from-function-space-c0-1-to%23new-answer', 'question_page');
                                }
                                );

                                Post as a guest















                                Required, but never shown





















































                                Required, but never shown














                                Required, but never shown












                                Required, but never shown







                                Required, but never shown

































                                Required, but never shown














                                Required, but never shown












                                Required, but never shown







                                Required, but never shown







                                Popular posts from this blog

                                1300-talet

                                1300-talet

                                Display a custom attribute below product name in the front-end Magento 1.9.3.8