如何将极坐标形式的x和y值从这些耦合的ODE转换为笛卡尔形式并绘制图形？

Question

我已经编写了这个代码来模拟弹簧摆的运动

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()

其中x和y是两个变量，我试图将x和y转换为极坐标n（x轴）和u（y轴），然后在图形上将图n和u转换为x在x上-axis和u在y轴上。但是，当我绘制上面的代码时，它给了我：

相反，我应该得到一个与此类似的图像：

代码的第一部分 - 从“def deriv（z，t）：到sol：odeint（deriv ...”是生成x和y的值的地方，然后我可以将它们转换为直角坐标和如何更改我的代码来执行此操作？我是Python的新手，因此我可能不了解某些术语。谢谢！

Answer 1

第一个解决方案应该给你预期的结果，但是在执行ode时有一个错误。

传递给odeint的函数应该返回一个包含一阶微分方程系统解的数组。

在你的情况下，你正在解决的是

相反，你应该解决

为此，请将代码更改为此

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()

代码的第二部分看起来像是在尝试更改坐标。我不确定你为什么要再次尝试解决颂歌，而不仅仅是这样做。

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))

截至目前，您正在做的是解决这个系统：

这导致x = e ^ t和y = e ^ t并且n'=（0.4 + e ^ t）* sin（e ^ t）u'=（0.4 + e ^ t）* cos（e ^ t）。

没有过多的细节，你可以看到，这将导致吸引子作为n的导数，你将开始以指数速率更快地切换符号，并导致n和u折叠到你的情节所示的吸引子。

如果你真的想要解决另一个微分方程，我需要看到它以便进一步帮助你

如果您进行转换并将时间设置为1000，则会发生这种情况：