I am trying to solve numerically with NDsolve a pretty standard system of non-linear ODE's of the form y'=F(y), where y is a list of n points in the plane, and F is a certain function defined inside the program. All goes well, until I reach NDsolve. My guess is that I have an elementary syntax error in passing values of F to NDsolve, steming from my lack of experience with working with functions in Mathematica (my programming background is in C...). Thanks ahead for any suggestions. I paste the program below.
Polygon flow
This is a numerical integration program that shows a "polygon flow" in the plane. That is, p={p[1],..., p[n]} is a list of n points, where each p[i] is in R^2, moving by a certain rule, given by the ODE p'=pflow[p]. The function pflow[p] requires a small program to evaluate it. After solving this ODE, with some initial condition p0,
the program shows how the n points move in the plane.
Preliminary definitions
In[1]:= Clear[n, canon, p0, long, e, pPlot, reflect, pPlot, pflow, A]
(* Define the number of vertices of the polygon; should be odd *)
n = 3;
(* Define the cannonical basis of R^2 *)
canon = Table[IdentityMatrix[2][[i]], {i, 1, 2}];
(* Define some initial polygon *)
p0 = Table[{2 Cos[i 2 Pi/n], Sin[i 2 Pi/n]} // N, {i, 1, n}];
(* This function "lengthen " p, by adding the 1st vertex as the nth+1 vertex \
*)
long[p_] := Join[p, {p[[1]]}];
(* Define the list of edges (unit vectors) of the poligon *)
e[p_] := Module[{edge, lp},
lp = long[p];
Table[
edge = lp[[i + 1]] - lp[[i]];
edge/Norm[edge], {i, 1, n}]
]
(* This generates a drawing of the polygon *)
pPlot[p_] := Graphics[Line[long[p]]]
(* This reflects a vector v about the edge e *)
reflect[v_, e_] := -v + 2 (v.e) e
(* This is the composition of reflections about the edges of the polygon *)
A[e_] := Module[{v},
Table[
For[
j = 1; v = canon[[i]],
j < n + 1, j++,
v = reflect[v, e[[j]]]
]; v,
{i, 1, 2}
]
]
(* This defines a vector field on the space of polygons *)
pflow[p_] := Module[
{edges, R, u},
edges = e[p];
R = A[edges];
(* R is a composition of an odd number of reflections, so is a reflection,
hence has a unique fixed direction.
We pick a unit vector in this direction at the 1st vertex,
then spread it around the polygon by reflections wrt the edges *)
u = First@NullSpace[R - IdentityMatrix[2]];
Table[
If[
i == 1,
u = u/Norm[u],
u = reflect[u, edges[[i - 1]]]
],
{i, 1, n}
]
]
Integrating the ODE p'=pflow[p]
In[11]:= (* Integration time *)
T = 10;
(*Integrating the ODE*)
s = NDSolve[{y'[t] == pflow[y[t]], y[0] == p0}, y, {t, 0, T}];
yy := y /. s;
Animate[
Graphics[
Line[long[yy[t]]],
PlotRange -> 3],
{t, 0, T},
AnimationRate -> .1
]
