0

我试图通过只取有限数量的点来绘制一个无限系列。就我而言，3分和10分就足够了。

e该方程是偏心率中的拉格朗日幂级数。

E = Me + \sum_{n = 1}^{\infty}a_n e ** n

a_n在哪里

a_n = (1 / 2 ** (n - 1) * \sum_{k = 0}^{\lfloor n/2\rfloor} (-1) ** k / 
      ((n - 2 * k)! * k!) * (n - 2 * k) ** (n - 1) * np.sin((n - 2 * k) * Me))

\lfloor n/2\rfloor乳胶的地板功能也是如此n/2。

自变量是E且依赖Me的，因此该函数不会像正常遇到此类函数那样编写，但我看不到一种明确解决的方法，Me以便我们可以编写Me(E)

所以到目前为止我所做的是（见下文）这是错误的，因为它不起作用。我该怎么做才能使代码和绘图正常工作？ 

import numpy as np
import pylab as py
import math
from scipy.misc import factorial as fact

Me = np.linspace(0, 2 * np.pi, 50000.0)
e = 0.65
a = [1.0 / 2.0 ** (math.floor(n / 2.0) - 1.0) *
     sum([(-1.0) ** math.floor(n / 2.0) /
          (fact(math.floor(n / 2.0) - k) * fact(k)) *
          (math.floor(n / 2.0) - 2.0 * k) ** (math.floor(n / 2.0) - 1.0) *
          np.sin((math.floor(n / 2.0) - 2.0 * k) * Me)
          for k in range(1, 4, 1)])
          for n in range (1, 4, 1)]

print a

def E2(x):
    return Me + sum(a[n] * e ** n for n in range(1, 4, 1)) - x

fig = py.figure()
ax = fig.add_subplot(111)
ax.plot(Me, E2(Me))
py.xlim((0, 2 * np.pi))
py.ylim((0, 2 * np.pi))
py.show()

有了这个程序，我得到

In [2]: /usr/bin/ipython:17: RuntimeWarning: divide by zero encountered in double_\
scalars
/usr/bin/ipython:17: RuntimeWarning: invalid value encountered in multiply
/usr/bin/ipython:17: RuntimeWarning: invalid value encountered in add
[array([ nan,  inf,  inf, ..., -inf, -inf, -inf]), array([ nan,  inf,  inf, ..., -\
inf, -inf, -inf]), array([ nan,  inf,  inf, ..., -inf, -inf, -inf])]

无穷大根本不应该是一个值，所以我不确定它是如何得出的。

最后的错误是list of index out of range

/home/dustin/Documents/School/UVM/Engineering/OrbitalMechanics/lagrangeseries.py i\
n <genexpr>((n,))
     17
     18 def E2(x):
---> 19     return Me + sum(a[n] * e ** n for n in range(1, 4, 1)) - x
     20
     21 fig = py.figure()

IndexError: list index out of range

这怎么超出范围？一切都是从 1 到 3 的总和？

4

2 回答 2

0

由于 Python 中除法的工作方式，我不需要使用 floor 函数。

import numpy as np
import pylab as py
from scipy.misc import factorial as fact

e = 0.65


def E(M):
    return (M + sum((1.0 / 2.0 ** (n - 1) *
                     sum((-1) ** (k) / (fact(n - k) * fact(k)) *
                         (n - 2 * k) ** (n - 1) *
                         np.sin((n - 2 * k) * M)
                         for k in range(0, n / 2, 1))) * e ** n
                         for n in range(1, 4, 1)))


M = np.linspace(0, 2 * np.pi, 50000.0)

fig = py.figure()
ax = fig.add_subplot(111)
ax.plot(E(M), M)
py.xlim((0, 2 * np.pi))
py.ylim((0, 2 * np.pi))
py.show()
于 2013-05-27T04:09:00.610 回答
0

对于 k 计数器，范围值中的索引应该增加一，而不是：

对于范围内的 k (0, n / 2, 1)

它应该是：对于 k in range(0, n / 2+1, 1)

否则对于 n =4 range(0, 2, 1) 根据需要返回 0,1 而不是 0,1,2。

于 2016-08-03T08:55:37.877 回答