Usuario:Arkin/Retratos de fase

Los siguientes macros de Mathematica producen diversos tipos de retratos de fase para sistemas de ecuaciones diferenciales. Si deseas el archivo .nb que los contiene puedes pedírmelos y te los enviaré vía correo electrónico.

---

El código siguiente produce un retrato de fase para un sistema de ecuaciones diferenciales (de dos ecuaciones). Los argumentos son

• dxdt: la expresión que define a $dx/dt$.
• dydt: la expresión que define a $dy/dt$.
• x0: el valor de $x_{0}=x(0)$
• y0: el valor de $y_{0}=y(0)$
• t0: el valor inicial que debe tomar $t$
• tmax: el valor máximo que debe tomar $t$

PhasePortrait[dxdt_, dydt_, x0_, y0_, t0_, tmax_] := Module[{DXDT, DYDT, eq, sol, grph}, 
DXDT := dxdt /. {x -> x[t], y -> y[t]}; 
DYDT := dydt /. {x -> x[t], y -> y[t]}; 
eq := {x'[t] == DXDT, y'[t] == DYDT, x[t0] == x0, y[t0] == y0}; 
sol := NDSolve[eq, {x, y}, {t, t0, tmax}];
grph := ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, t0, tmax}, 
PlotRange -> All, PlotStyle -> {{Black, Thick}}
PlotPoints -> 500]; Show[grph]];

Por ejemplo, para el sistema de ecuaciones diferenciales de Van der Pol,

${dx \over dt}=y$

${dy \over dt}=-x+(1-x^{2})y$

podríamos escribir

PhasePortrait[y, -x + (1 - x^2) y, .5, .5, 0, 20]

para obtener

A diferencia del macro anterior, la función DirectedPhasePortrait, definida por el código que aparece a continuación, incluye flechas en la gráfica que indican la dirección en la que $x$ e $y$ cambian conforme $t$ crece. Puesto que la gráfica se produce de manera diferente, es necesario especificar un parámetro adicional step para el valor de $\Delta t$. Entre menor sea este, mas puntos serán calculados y mejor será la calidad de la gráfica, pero Mathematica tardará más tiempo en relizar la evaluación de la función.

DirectedPhasePortrait[dxdt_, dydt_, x0_, y0_, t0_, tmax_, step_] := 
 Module[{DXDT, DYDT, eq, sol, tablepoints, points, tablearrowpoints, 
   tablearrowpointsdt, arrowpoints, grph, arrows}, 
  DXDT := dxdt /. {x -> x[t], y -> y[t]}; 
  DYDT := dydt /. {x -> x[t], y -> y[t]}; 
  eq := {x'[t] == DXDT, y'[t] == DYDT, x[t0] == x0, y[t0] == y0}; 
  sol := NDSolve[eq, {x, y}, {t, t0, tmax}]; 
  tablepoints := 
   Table[Evaluate[{x[t], y[t]} /. sol], {t, t0, tmax, 
     ToExpression[step]}];
  points := Table[Extract[tablepoints, {n, 1}], {n, 1, tmax/step}];
  tablearrowpoints := 
   Table[Evaluate[{x[t], y[t]} /. sol], {t, t0, 
     tmax, (tmax - t0)/4}];
  tablearrowpointsdt := 
   Table[Evaluate[{x[t], y[t]} /. sol], {t, t0, 
     tmax, ((tmax - t0 + .001)/4)}];
  arrowpoints := 
   Table[{{Extract[tablearrowpoints, {n, 1, 1}], 
      Extract[tablearrowpoints, {n, 1, 2}]}, {Extract[
       tablearrowpointsdt, {n, 1, 1}], 
      Extract[tablearrowpointsdt, {n, 1, 2}]}}, {n, 1, 4}];
  grph := Graphics[{Thick, Line[points]}];
  arrows := 
   Graphics[{Thickness[0], Arrowheads[Large], Arrow /@ arrowpoints}];
  Show[{arrows, grph}, Axes -> True, AxesOrigin -> {0, 0}]
  ]

Un ejemplo de su uso sería

DirectedPhasePortrait[y, -x + (1 - x^2) y, .5, .5, 0, 20, .05]

que produce

El código sirve para crear espacios de fase. Es decir, se incluye una dimensión más correspondiente a los valores de $t$. La sintaxis es la misma que la de PhasePortrait excepto por el ultimo argumento opcional: RedFact. Ya que en ocasiones los valores máximos de $x$ e $y$ son, o muy grandes o muy pequeños en relación al valor máximo de $t$, a veces es necesario reducir o aumentar la escala del eje correspondiente a $t$. Este es precisamente el trabajo que realiza el argumento RedFact.

PhasePortrait3D[dxdt_, dydt_, x0_, y0_, t0_, tmax_, RedFact_: 1] := 
Module[{DXDT, DYDT, eq, sol, grph, grphx, grphy, xtwall, ytwall, 
xywall, aux, auy, xmax, xmin, ymax, ymin, twall, grpht}, 
DXDT := dxdt /. {x -> x[t], y -> y[t]}; 
DYDT := dydt /. {x -> x[t], y -> y[t]};
eq := {x'[t] == DXDT, y'[t] == DYDT, x[t0] == x0, y[t0] == y0}; 
sol := NDSolve[eq, {x, y}, {t, t0, tmax}];
grph := ParametricPlot3D[Evaluate[{x[t], y[t], t/RedFact} /. sol], {t, t0, tmax}, 
PlotRange -> All, PlotStyle -> {{Thick, GrayLevel[.9]}}, ColorFunction -> 
Function[{x, y}, ColorData["AtlanticColors"][y]], ImageSize -> 500]; 
grphx := Table[ParametricPlot3D[Evaluate[{x[t], b, t/RedFact} /. sol], {t, t0, tmax}, 
PlotStyle -> {Thick, Opacity[.2]}, ColorFunction -> 
Function[{x, y}, ColorData["AtlanticColors"][x]]], {b, ymin - 3, ymin - 3}]; Table[ParametricPlot3D[
Evaluate[{b, y[t], t/RedFact} /. sol], {t, t0, tmax}, PlotStyle -> {Thick, Opacity[.2]}, 
ColorFunction -> Function[{x, y}, ColorData["AtlanticColors"][y]]], {b, xmin - 3, xmin - 3}];
auy := ParametricPlot[Evaluate[{t, y[t]} /. sol], {t, t0, tmax}];
aux := ParametricPlot[Evaluate[{t, x[t]} /. sol], {t, t0, tmax}];
xmax := Max[Last /@ Level[Cases[aux, _Line, Infinity], {-2}]];
xmin := Min[Last /@ Level[Cases[aux, _Line, Infinity], {-2}]];
ymax := Max[Last /@ Level[Cases[auy, _Line, Infinity], {-2}]];
ymin := Min[Last /@ Level[Cases[auy, _Line, Infinity], {-2}]];
PlotStyle -> {{Gray, Opacity[.05]}}], {b, ymin - 3, ymin - 3}]; 
ytwall := Table[ParametricPlot3D[{b, t, s}, {t, ymin - 3, ymax + 3}, {s, t0,tmax/RedFact}, Mesh -> False, 
PlotStyle -> {{Gray, Opacity[.05]}}], {b, xmin - 3, xmin - 3}];
twall := Table[ParametricPlot3D[{t, s, b}, {s, ymin - 3, ymax + 3}, {t, xmin - 3, xmax + 3}, Mesh -> False, 
PlotStyle -> {{Gray, Opacity[.05]}}], {b, t0, t0}];
grpht := Table[ParametricPlot3D[Evaluate[{x[t], y[t], t0} /. sol], {t, t0, tmax}, 
PlotStyle -> {Thick, Opacity[.2]}, 
ColorFunction -> Function[{x, y}, ColorData["AtlanticColors"][y]]], {b, xmin - 3, xmin - 3}];
Show[{grph, grphx, grphy, xtwall, ytwall, twall, grpht}, 
ViewPoint -> {2, 2, 2}, AxesLabel -> {x, y, t}]];

Por ejemplo, el sistema de ecuaciones diferenciales de Van der Pol anterior se vería así usando PhasePortrait3D en lugar de PhasePortrait

PhasePortrait3D[y, -x + (1 - x^2) y, .5, .5, 0, 20, 2]

El código siguiente sirve para exportar los puntos que usa Mathematica para graficar el retrato de fase del sistema de ecuaciones. Los argumentos usados son los mismos, excepto por los dos argumentos finales adicionales step y filename. El argumento step se usa para especificar la magnitud de $\Delta t$. Entre menor sea este, más puntos serán calculados. El argumento filename se usa para especificar el nombre del archivo, cuya extensión será .dat

ExportPhasePortrait[dxdt_, dydt_, x0_, y0_, t0_, tmax_, step_, filename_] := 
Module[{DXDT, DYDT, eq, sol, tablepoints, points}, 
DXDT := dxdt /. {x -> x[t], y -> y[t]}; 
DYDT := dydt /. {x -> x[t], y -> y[t]}; 
eq := {x'[t] == DXDT, y'[t] == DYDT, x[t0] == x0, y[t0] == y0}; 
sol := NDSolve[eq, {x, y}, {t, t0, tmax}]; 
tablepoints := Table[Evaluate[{x[t], y[t]} /. sol], {t, t0, tmax, ToExpression[step]}];
points := Table[Extract[tablepoints, {n, 1}], {n, 1, tmax/step}];
Grid[{{Graphics[Point /@ points, ImageSize -> 300, Axes -> True],
Export[ToString[filename] <> ".dat", points]}}, Frame -> All]];