Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Iterative map in mathematica

  1. Sep 28, 2011 #1
    I have no programming experience and trying to get mathematica to do what I need it to do is frustrating. I have the following functions that I need to iterate. For notational purposes, k[t+1] is the value of K in the next period. w is a parameter.

    k[t+1] = -x[t] - y[t] + w
    x[t+1] = k[t]*sqrt(y[t]/(-y[t] + w)) - y[t]
    y[t+1] = k[t]*sqrt(x[t]/(-x[t] + w)) - x[t]

    As I care more about the behavior of x and y, this system can be reduced to:

    x[t+2] = (-x[t]-y[t]+w)*sqrt(y[t+1]/(-y[t+1] + w)) - y[t+1]
    y[t+2] = (-x[t]-y[t]+w)*sqrt(x[t+1]/(-x[t+1] + w)) - x[t+1]

    I have already found for what w the system becomes unstable by taking the Jacobian of the 3-system and seeing where it exceeds unity in absolute value. Now I just want to iterate the system and get some pretty graphs.
  2. jcsd
  3. Sep 29, 2011 #2
    NestList will generate an iterative map for you. You can look up the documentation for that to get an idea how it works in general. For your specific example problem this might do:

    k0 = .8; x0 = .7; y0 = .4; w = .3;
    k[v_] := v[[1]]; x[v_] := v[[2]]; y[v_] := v[[3]];
    f[t_] := {-x[t]-y[t]+w, k[t]*Sqrt[y[t]/(-y[t]+w)]-y[t], k[t]*Sqrt[x[t]/(-x[t]+w)]-x[t]}
    r = NestList[f, {k0, x0, y0}, 2]

    You should be able to check that those match the equations you gave.

    So f inside the NestList is your function to repeatedly apply. You didn't provide starting conditions so I used {k0,x0,y0} as the initial value to start with and the final 2 says to only do two steps, just so you can see how it starts. You can increase the 2 once you believe it is working.

    After you are done you will have a list of results, each containing {k,x,y} values.
    If all you want is a list of results, each containing {x,y} values then you could do
    and that would discard the k value from each result.

    I am not certain where you want to go with this or how much of what I've written will be understandable to someone with no Mathematica background.

    Study what I've written here, see if you can try it out. Carefully manually verify the calculations for a few steps with a few different starting values to try to expose any mistakes in this. If there are things you don't understand then provide more information about what you need explained and what you need this to do.

    It might be less error prone to "scrape" the code off the screen and paste it into Mathematica than to try to type it back in because Mathematica is very demanding about () versus {} versus [] and capitalization.

    If you get a list of real number pairs {x,y} then you can use ListPlot to directly provide your graphs from this.
    Last edited: Sep 29, 2011
  4. Sep 30, 2011 #3

    k0 = .8; x0 = .7; y0 = .4; w = .3;
    f[{k_, x_, y_}] := {-x-y+w, k*Sqrt[y/(-y+w)]-y, k*Sqrt[x/(-x+w)]-x}
    r = NestList[f, {k0, x0, y0}, 2]
    Last edited: Sep 30, 2011
  5. Oct 2, 2011 #4
    Sorry for getting back to you on this late, been busy!

    For the most part this is working as I would like. I put in:

    k0 = .6; x0 = .3; y0 = .1 ; w = 1;
    f[{k_, x_, y_}] := {-x - y + w, k*Sqrt[y/(-y + w)] - y,
    k*Sqrt[x/(-x + w)] - x}
    r = NestList[f, {k0, x0, y0}, 100]

    The numbers converge to precisely what my math has told me they should converge to, and the graph it generates is decent. However, it would be better if I could see the dots generated in time rather than just all the dots plotted at the same time. Difficult to see behavior when they are all plastered on there at the same time.

    X and Y denote production of 2 firms. I have also adapted this to handle 3 firms, so X, Y, and Z output (K is kinda a weird variable, but it helps to determine what the firms should do). The math tells me what it should converge to, but when I try to use the program to iterate the 3 firm case it outputs absurd answers, such as negative output and even complex valued output. The only thing I can reason for why its doing this is that it perhaps is not taking the values to enough significant figures.
  6. Oct 3, 2011 #5
    Hmm. For most of the similar models that I've been dealing with, it was possible to write the "reaction functions" (the functions I wrote in my first post) in such a way that you could actually have X in the next period determined by X in the last period. The test for stability of the continuously differentiable FUNCTION involved taking the Jacobian of the reaction functions and seeing if it passed unity in absolute value.

    You cannot put my functions together in such a way that X in the next period is based on X in the last period. So perhaps I have a continuously differentiable MAP? The test for stability of a continuously differentiable map is slightly different, and perhaps if I use this test I can see why the 2 firm case is stable and the 3 firm case is not stable.

    That is, if I am understanding the difference between a continuously differentiable function and continuously differentiable map.
    Last edited: Oct 3, 2011
  7. Oct 3, 2011 #6
    1. Your square roots and minus signs make it feasible and perhaps even easy to end up with complex values during your iterations. I can't say anything about whether your map is appropriate or not.

    2. Users occasionally ask if they can't somehow slow Mathematica down and watch it plot one point at a time. There has been nothing in the past to make that easy, but there are a handful of things you might try. There is the PlotJoined->True option that you can add as an argument to your ListPlot. That will connect the points in order with line segments. I don't know that will help you enough. My simple experimentation with your equations shows wild behavior of the iterations, even if I use ListPlot[Map[Re,yourlist],PlotJoined->True] to only see the real components.

    You might also explore Manipulate which is a new feature introduced in version 8 that might allow you to only plot points up to some boundary that you can select with a slider and see if you can figure out enough to be able to use that. That might involve the Manipulate containing your ListPlot containing a Take to extract the first part of of your list and the Manipulate control telling Take how many points to extract.

    The older way of doing that was to make a large collection of plots, each with one more point than the last and then export the whole collection as a movie.

    But doing things like this is perhaps more advanced than I expected you to be able to figure out how to implement.

    From what you have written I'd probably look first at making certain your function is correct and only after that has been verified would I start trying to wrestle Mathematica to the ground.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook