2013-05-17 2 views
0

Я пытаюсь построить бесконечную серию, беря только конечное количество точек. В моем случае достаточно 3 и 10 баллов.Python: Построение бесконечной серии только за счет конечных точек

Уравнение представляет собой линейный ряд Лагранжа в 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?

+0

много ошибок, для стартеров 1/2 ** (n-1) оценивается в 0, если n является целым числом. вам нужны квадратные скобки на встроенных петлях. google "pyhon tutorial" принесет тонны хорошей информации, чтобы вы начали .. – agentp

+0

@george Я добавил десятичные точки ко всему, чтобы исправить целочисленную проблему. – dustin

+0

@george Я не могу найти что-либо в линейных циклах. Он просто говорит, что '[]' are list и '()' являются кортежами, где список может быть изменен, а кортежи cant – dustin

ответ

0

Мне не нужно было использовать функцию пола из-за того, как деление работает на Python.

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

индексы в значении диапазона должны быть увеличены на единице для K счетчика, поэтому вместо того, чтобы:

для к в диапазоне (0, п/2, 1)

он должен : для k в диапазоне (0, n/2 + 1, 1)

в противном случае для диапазона n = 4 (0, 2, 1) возвращает 0,1, а не 0,1,2 по мере необходимости.

+0

Это не дает ответа на вопрос. Когда у вас будет достаточно [репутации] (http://stackoverflow.com/help/whats-reputation), вы сможете [прокомментировать любое сообщение] (http://stackoverflow.com/help/privileges/comment); вместо этого [предоставить ответы, которые не требуют разъяснений у аськи) (http://meta.stackexchange.com/questions/214173/why-do-need-50-reputation-to-comment-what-can- я-делать-вместо этого). – ddb

Смежные вопросы