2
$\begingroup$

I've animated a dynamics as follows in Python. However, I'd like to have the same in Mathematica (whose graphics output, in my view, is fancier than Matplotlib).

#!/usr/bin/env python
# -*- coding: utf-8 -*- 

from scipy.integrate import odeint
import matplotlib.pyplot as plt
import numpy as np
import matplotlib as mpl
import mpl_toolkits.mplot3d.axes3d as p3
from matplotlib import animation
plt.style.use('seaborn-pastel')

def f(p, t): 
      x1 = p[0]
      y1 = p[1]
      x2 = p[2]
      y2 = p[3]
      dx1dt = x2 - x1
      dy1dt = y2 - y1
      dx2dt = x1 - x2
      dy2dt = y1 - y2
      return np.array ([dx1dt, dy1dt, dx2dt, dy2dt])

ts = np.arange (0 , 20 , 0.1)
x10 = -5
y10 = -7
x20 = -5
y20 = 7

p = np.array ([x10 , y10, x20 , y20])

s = odeint (f , p , ts)
x1s = s [: , 0]
y1s = s [: , 1]
x2s = s [: , 2]
y2s = s [: , 3] 

fig = plt.figure()
ax = p3.Axes3D(fig)

def update_lines(num, dataLines, lines):
    for line, data in zip(lines, dataLines):
        line.set_data(data[0:2, :num])
        line.set_3d_properties(data[2, :num])
        line.set_marker("o")
        #line.set_label("p1")
        #plt.legend()
    return lines

data = np.array([[x1s,y1s,ts],[x2s,y2s,ts]])

lines = [ax.plot(dat[0, 0:1], dat[1, 0:1], dat[2, 0:1])[0] for dat in data]

ax.set_xlim(-10.1,10.1)
ax.set_ylim(-10.1,10.1)
ax.set_zlim(0,20)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("t")
ax.legend()

anim = animation.FuncAnimation(fig, update_lines, 50, fargs=(data, lines),
                                   interval=200, blit=True, repeat=True)

plt.rcParams['animation.ffmpeg_path'] ='C:\\...\\FFmpeg\\bin\\ffmpeg.exe'
FFwriter = animation.FFMpegWriter(fps = 100)
anim.save('test.mp4', writer = FFwriter)

plt.show()

As a newbie in using Mathematica, here is what I have planned to recreate the animation above.

1- Numerically integrating the system using DSolve[]

DSolve[{x1'[t] == -x1[t] + x2[t], 
  y1'[t] == -y1[t] + y2[t], 
  x2'[t] == x1[t] - x2[t], 
  y2'[t] == y1[t] - y2[t], -5, -7, -5, 7}, x1[t], y1[t],
  x2[t], y2[t], t]

2- Animating the integration result using Animate[Plot3D[...]]

But, I don't know how to set variation ranges corresponding to my variables. Can you please give a clue?


Edit:

The desired script does numerically integrate the following system of equations

dx1dt = x2 - x1
dy1dt = y2 - y1
dx2dt = x1 - x2
dy2dt = y1 - y2

According to the initial conditions below

x10 = -5
y10 = -7
x20 = -5
y20 = 7

Then the overall dynamics of the system is animated by varying the time t between $0$ to $20$ units.

$\endgroup$
  • 2
    $\begingroup$ Can you please explain what your script is doing and what is the desired result for all those who do not speak python. $\endgroup$ – yarchik Jul 21 at 12:42
  • $\begingroup$ @yarchik: Please check the edit $\endgroup$ – Roboticist Jul 21 at 12:46
  • $\begingroup$ I still do not understand. You solve ODE for 4 functions. They all depend on t. How do you make a 3D plot out of them? 3D plot needs one 1 function that depends on 2 parameters. $\endgroup$ – yarchik Jul 21 at 12:59
  • $\begingroup$ @yarchik: The result is actually the trajectories of the simulated agents. What is actually drawn is a set of tuples (x,y,t) in a 3D space. $\endgroup$ – Roboticist Jul 21 at 13:02
5
$\begingroup$

You can use NDSolveValue and Animate:

ClearAll[x1, x2, y1, y2]
ndsv = NDSolveValue[{x1'[t] == -x1[t] + x2[t], 
    y1'[t] == -y1[t] + y2[t], 
    x2'[t] == x1[t] - x2[t], 
    y2'[t] == y1[t] - y2[t], 
    x1[0] == -5, x2[0] == -7, y1[0] == -5,  y2[0] == 7}, 
   {x1, x2, y1, y2}, 
   {t, 0, 20}, 
   "ExtrapolationHandler" -> {Undefined &, "WarningMessage" -> False}];


Animate[Plot[Evaluate@Through@ndsv[t], {t, 0, tmax}], {tmax, .1, 20}]

enter image description here

frames = Table[Plot[Evaluate@Through@ndsv[t], {t, 0, tmax}], {tmax, .1, 20, .1}];

Export["frms.gif", frames]

enter image description here

frames2 = Table[Plot[Evaluate@Through@ndsv[t], {t, 0, tmax}, 
    PlotRange -> {{0, 20}, {-10, 10}}], {tmax, .1, 20, .1}];

Export["frms2.gif", frames2]

enter image description here

Update: For 3D visualization, you can use ParametricPlot3D with {{x1[t], y1[t], t}, {x2[2], y2[t], t}} as the first argument:

Animate[ParametricPlot3D[{{ndsv[[1]]@t, ndsv[[3]]@t, t}, 
    {ndsv[[2]]@t, ndsv[[4]]@t, t}}, {t, 0, tmax}, 
  BoxRatios -> 1,  BaseStyle -> Thick, PlotRange -> {{-10, 10}, {-10, 10}, {0, 20}}], 
 {tmax, .1, 20}]

enter image description here

$\endgroup$
  • $\begingroup$ Thanks, but the result needs to be 3D, say, $(x,y,t)$ points. $\endgroup$ – Roboticist Jul 21 at 13:06
  • $\begingroup$ @Roboticist, for the error message add ClearAll[x1,x2, y1,y2] before the code. $\endgroup$ – kglr Jul 21 at 13:09
  • 1
    $\begingroup$ @Roboticist, for 3D visualization please see the update. $\endgroup$ – kglr Jul 21 at 13:19
  • $\begingroup$ I do appreciate your help! $\endgroup$ – Roboticist Jul 21 at 13:19
  • 1
    $\begingroup$ @Roboticist, my pleasure. Welcome to mma.se. $\endgroup$ – kglr Jul 21 at 13:20
2
$\begingroup$

The OP mentions "numerically integrating...using DSolve[]." DSolve[] does a symbolic integration, which I will show, and I will also show numerical integration using NDSolve[], which can used to construct a faster (more responsive) Animate[]/Manipulate[] visualization (because all of the points are pre-calculuated and stored).

Symbolic solution:

ClearAll[x1, x2, y1, y2, t]
dsol = DSolveValue[  (* integrate *)
  {x1'[t] == -x1[t] + x2[t], y1'[t] == -y1[t] + y2[t], 
   x2'[t] == x1[t] - x2[t], y2'[t] == y1[t] - y2[t], x1[0] == -5, 
   x2[0] == -7, y1[0] == -5, y2[0] == 7},
  {x1, y1, x2, y2}, {t, 0, 20}]
(*
  {Function[{t}, -6 + E^(-2 t)], Function[{t}, 1 - 6 E^(-2 t)], 
   Function[{t}, -6 - E^(-2 t)], Function[{t}, 1 + 6 E^(-2 t)]}
*)

P1 = Function[t, Evaluate[Flatten@{t, Through[dsol[[{1, 2}]][t]]}]]
P2 = Function[t, Evaluate[Flatten@{t, Through[dsol[[{3, 4}]][t]]}]]
(* 
  Function[t, {t, -6 + E^(-2 t), 1 - 6 E^(-2 t)}]
  Function[t, {t, -6 - E^(-2 t), 1 + 6 E^(-2 t)}]
*)

Animate[
 ParametricPlot3D[{P1[tt], P2[tt]}, {tt, 0, t + $MachineEpsilon}, 
  PlotRange -> {{0, 20}, {-11, -1}, {-5, 7}}, 
  PlotRangePadding -> Scaled[0.05],
  (* Manipulate/Animate localize t which gives it a ugly symbol name *)
  (* The string processed by ToExpression won't be mapped to the local var. *)
  (* One could simply use strings, but 
  (*   this way, you get the nice italic look of TraditionalForm *)
  AxesLabel -> (ToExpression[#, StandardForm, HoldForm] & /@ {"t", "x", "y"})],
 {t, 0, 20}
 ]

Numerical solution:

The trick here is to construct the integral curves from the values computed by NDSolve[] (NDSolveValue[] is just syntactic sugar for NDSolve[]). The evolution of the integral curve can be done by selecting an initial segment of each curve up to the current animation time.

ndsol = NDSolveValue[{
    x1'[t] == -x1[t] + x2[t], y1'[t] == -y1[t] + y2[t], 
    x2'[t] == x1[t] - x2[t], y2'[t] == y1[t] - y2[t], x1[0] == -5, 
    x2[0] == -7, y1[0] == -5, y2[0] == 7,
    tt'[t] == 1, tt[0] == 0},  (* dummy int. to get time steps more easily *)
   {x1, y1, x2, y2, tt}, {t, 0, 20}];

(* get points from values computed by NDSolve[] *)
(* xy coords. of pt. 1 are parts {1, 2}, of pt. 2 are parts {3, 4}; time is part 5 *)
values = Through[ndsol["ValuesOnGrid"]];
points = {Transpose@values[[{5, 1, 2}]], Transpose@values[[{5, 3, 4}]]};

(* make a nicely scaled PlotRange - I cheated and also used this above *)
plotrng = Apply[{#5, MinMax[{#1, #3}], MinMax[{#2, #4}]} &, MinMax /@ values];
With[{m = Max[-Subtract @@@ plotrng]},
 plotrng =    (* expand zero-length ranges -- not needed in this specific case *)
  plotrng /. {a1_Real, a2_Real} /; a2 == a1 :> {a1 - 0.5 m, a2 + 0.5 m};
 plotrng = Rescale[#, Mean[#] - Subtract @@ # {-1, 1}, 
     Mean[#] + Max[0.5 m, -Subtract @@ #] {-1, 1}] & /@ plotrng
 ]
(*  {{0., 20.}, {-11., -1.}, {-5., 7.}}  *)

Animate[
 Graphics3D[
  {Thick, Riffle[
    ColorData[97] /@ Range@2,
    Line /@ points[[All, ;; LengthWhile[values[[5]], # <= t &]]] (* segments of curves *)
    ]},
  PlotRange -> plotrng, PlotRangePadding -> Scaled[0.05],
  Boxed -> True, Axes -> True,
  AxesLabel -> (ToExpression[#, StandardForm, HoldForm] & /@ {"t", "x", "y"})
  ],
 {t, Sequence @@ plotrng[[1]]}
 ]

Both animations produce an output like the one below. The second is more responsive because there is very little computation compared to ParametricPlot3D.

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.