How do I convert the x and y values in polar form from these coupled ODEs to to cartesian form and graph them?

Question

I have written this code to model the motion of a spring pendulum

import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi, array
import matplotlib.pyplot as plt

def deriv(z, t):
    
    x, y, dxdt, dydt = z

    dx2dt2=(0.415+x)*(dydt)**2-50/1.006*x+9.81*cos(y)
    dy2dt2=(-9.81*1.006*sin(y)-2*(dxdt)*(dydt))/(0.415+x)
    
    return np.array([x,y, dx2dt2, dy2dt2])

init = array([0,pi/18,0,0])
time = np.linspace(0.0,10.0,1000)
sol = odeint(deriv,init,time)

def plot(h,t):
    n,u,x,y=h
    n=(0.4+x)*sin(y)
    u=(0.4+x)*cos(y)
    return np.array([n,u,x,y])

init2 = array([0.069459271,0.393923101,0,pi/18])

time2 = np.linspace(0.0,10.0,1000)
sol2 = odeint(plot,init2,time2)

plt.xlabel("x")
plt.ylabel("y")
plt.plot(sol2[:,0], sol2[:, 1], label = 'hi')
plt.legend()
plt.show()

where x and y are two variables, and I'm trying to convert x and y to the polar coordinates n (x-axis) and u (y-axis) and then graph n and u on a graph where n is on the x-axis and u is on the y-axis. However, when I graph the code above it gives me:

Instead, I should be getting an image somewhat similar to this:

The first part of the code - from "def deriv(z,t): to sol:odeint(deriv..." is where the values of x and y are generated, and using that I can then turn them into rectangular coordinates and graph them. How do I change my code to do this? I'm new to Python, so I might not understand some of the terminology. Thank you!

Answer 1

The first solution should give you the expected result, but there is a mistake in the implementation of the ode.

The function you pass to odeint should return an array containing the solutions of a 1st-order differential equations system.

In your case what you are solving is

While instead you should be solving

In order to do so change your code to this

import numpy as np
from scipy.integrate import odeint
from numpy import sin, cos, pi, array
import matplotlib.pyplot as plt


def deriv(z, t):

    x, y, dxdt, dydt = z

    dx2dt2 = (0.415 + x) * (dydt)**2 - 50 / 1.006 * x + 9.81 * cos(y)
    dy2dt2 = (-9.81 * 1.006 * sin(y) - 2 * (dxdt) * (dydt)) / (0.415 + x)

    return np.array([dxdt, dydt, dx2dt2, dy2dt2])


init = array([0, pi / 18, 0, 0])
time = np.linspace(0.0, 10.0, 1000)
sol = odeint(deriv, init, time)

plt.plot(sol[:, 0], sol[:, 1], label='hi')
plt.show()

The second part of the code looks like you are trying to do a change of coordinate. I'm not sure why you try to solve the ode again instead of just doing this.

x = sol[:,0]
y = sol[:,1]

def plot(h):
    x, y = h
    n = (0.4 + x) * sin(y)
    u = (0.4 + x) * cos(y)
    return np.array([n, u])

n,u = plot( (x,y))

As of now, what you are doing there is solving this system:

Which leads to x=e^t and y=e^t and n' = (0.4 + e^t) * sin(e^t) u' = (0.4 + e^t) * cos(e^t).

Without going too much into the details, with some intuition you could see that this will lead to an attractor as the derivative of n and u will start to switch sign faster and with greater magnitude at an exponential rate, leading to n and u collapsing onto an attractor as shown by your plot.

If you are actually trying to solve another differential equation I would need to see it in order to help you further

This is what happen if you do the transformation and set the time to 1000: