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]
]
# ========================= 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), );
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):
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)