Take the 2-minute tour ×

How to draw the image of a circle $x^2+(y-1)^2<1/4$ under the action of a transformation of the phase flow for the equation $\dot{x}=y,\ \dot{y}=-\sin x$? Here $\dot{x}$ means $dx/dt$. Any help or suggestions will be appreciated!

The following code shows how to plot the trace of a point under that action which may be helpful.

splot = 
  StreamPlot[{y, -Sin[x]}, {x, -4, 4}, {y, -3, 3}, StreamColorFunction -> "Rainbow"];
      Evaluate[First[{x[t], y[t]} /. 
        NDSolve[{x'[t] == y[t], y'[t] == -Sin[x[t]], Thread[{x[0], y[0]} == point]}, 
          {x, y}, {t, 0, T}]]], 
      {t, 0, T}, 
      PlotStyle -> Red]], 
  {{T, 2}, 1, 20}, 
  {{point, {3, 0}}, Locator}, 
  SaveDefinitions -> True]
share|improve this question
Take a look at a related post Basins of Attraction. I guess the both answers can be helpful. –  Artes Oct 26 '13 at 12:40
@Artes Thanks! It's really helpful but I can not apply it in this question. –  Eden Harder Oct 26 '13 at 13:00

2 Answers 2

up vote 11 down vote accepted

One way is to create a Polygon and transform the vertices under the flow. I used NDSolve to solve the flow for initial points in a square containing the OP's disk. Then I made a listable Function that can be applied to a the vertices of the polygon.

Since only the flow depends on the time t, I used a GraphicsComplex so that the vertices come first as a list. That way I can apply flow to just that part and wrap it in Dynamic. That limits the updating (by the Kernel) to just the vertices (plus the redrawing by the Front End).

flow = Function[{t, a, b},
   Evaluate[{x[t, a, b], y[t, a, b]} /. 
     First @ NDSolve[{D[x[t, a, b], t] == y[t, a, b], 
        D[y[t, a, b], t] == -Sin[x[t, a, b]], 
        x[0, a, b] == a, y[0, a, b] == b},
        {x, y}, {t, 0, 20}, {a, -1/2, 1/2}, {b, 1/2, 3/2}]],
circlePts = With[{n = 100}, 
    Table[{0., 1.} + {Cos[t], Sin[t]}/2, {t, 0, 2 Pi - 2 Pi/n, 2 Pi/n}]];

    Dynamic @ flow[t, circlePts[[All, 1]], circlePts[[All, 2]]],
    {Opacity[0.6], Red, Polygon[Range@Length@circlePts]}
  AspectRatio -> Automatic
 {t, 0., 20.}

Manipulate animation

share|improve this answer

This can also be done using ParametricPlot, as I will show in code here below. Special thanks to Michael E2 for sorting out a bug and proposing polar coordinates. The alternative would have been to use RegionFunction.

Clear[x, y]

sol[px_, py_] := sol[px, py] = First@NDSolve[{
     x'[t] == y[t],
     Thread[{x[0], y[0]} == {px, py}]
     }, {x, y}, {t, 0, 20}]

pos[px_, py_, T_] := Through[({x, y} /. sol[px, py])[T]]

pp[t_] := pp[t] = ParametricPlot[
   pos[r Cos[theta], 1 + r Sin[theta], t], {r, 0, 0.5}, {theta, 0, 
    2 Pi}, PlotPoints -> {3, 30}, Mesh -> False


Array[pp, 20];

Manipulate[Show[splot, pp[t]], {t, 0, 20, 1}]

Observe that memoization is used to make this fast, that is why I have that seemingly useless line Array[pp, 20] in there. It ensures that once the manipulate is displayed it will be entirely smooth. It's not necessary and depends on one's purpose. I was initially worried about speed, but I am happy to report that with memoization in place it is not very slow at all.

share|improve this answer
Thanks very much! But why you write 'sol[px_, py_] := sol[px, py] =' and 'pp[t_] := pp[t] =' to define functions? –  Eden Harder Oct 27 '13 at 1:05
@EdenHarder This is a case of simple memoization to improve speed, read this for details: mathematica.stackexchange.com/a/25150/731 –  Pickett Oct 27 '13 at 1:12

Your Answer


