pythonmatplotlibplotwolfram-mathematicacode-conversion

Convert a parameteric plot from mathematica to python


I have a small parametric plot that I am trying to convert to python from Mathematica. The problem is that my python script does not match this, I am struggling to get the list plot to work as it would be used in Mathematica.

How would you write this in python using matplotlib?

With[{r = 2.1, dt = 0.5, tfinal = 7.0},
 x = 1;
 y = 0;
 xylist = {{x, y}};
 Do[
  x = x - r y dt;
  y = y + r x dt;
  AppendTo[xylist, {x, y}];
  , {i, 0, tfinal, dt}];
 p1 = ParametricPlot[{Cos[t], Sin[t]}, {t, 0, 2 \[Pi]}, 
   PlotStyle -> Dashed];
 p2 = ListPlot[xylist, Joined -> True, AspectRatio -> 1];
 Show[p1, p2, PlotRange -> All]
 ]

enter image description here

# ========================= Python  ========================= #
# Parameters

r = 2.1
tfinal = 7.0
dt = 0.5
n = int(tfinal/dt)

# Containers
tspan = np.linspace(0, tfinal, n)
xtraj = np.zeros(n+1)
ytraj = np.zeros(n+1)


x0 = 1 # Initialize X
y0 = 0 # Initialize Y

value = odeint(f1, [x0,y0],tspan) # ODE values

# Euler scheme to calculate trajectory

# Initialize
xtraj[0] = x0
ytraj[0] = y0
for i in range(n):
    xtraj[i+1] = xtraj[i] - r*ytraj[i]*dt 
    ytraj[i+1] = ytraj[i] + r*xtraj[i]*dt 

fig, ax1 = plt.subplots(figsize=(9, 9))

# Plot ODE result portion

plt.plot(xtraj ,ytraj ,label='simulation $\Delta t = 0.1$, $r=1$, $T = 7$ ', color='#96CDCD', marker='o', linestyle='dashed', alpha=0.7)

theta = np.linspace( 0 , 2 * np.pi , 150 )
radius = 1
 
a = radius * np.cos( theta ) + 0
b = radius * np.sin( theta ) + 0
 
ax1.plot( a, b , label='actual, r=1', c='k')
ax1.set_aspect( 1 )
ax1.set(xlim=(-1.5, 1.5), ylim=(-1.5, 1.5), );


Solution

  • The issue is that in your Mathematica, the value of x is being updated each iteration through the "do" loop. Replace

    ytraj[i+1] = ytraj[i] + r*xtraj[i]*dt
    

    with

    ytraj[i+1] = ytraj[i] + r*xtraj[i+1]*dt
    

    and you'll see the same results (up to plot formatting, like color and linestyle, of course):

    enter image description here

    And note that as it stands, your python is a correctly implemented Euler integrator; it's your Mathematica that isn't quite 'Eulerian' (it's close to an Euler-Cromer method, though not quite)