Лабораторная работа 2. Операторы повторения. Приближенное вычисление числа $\pi$

Вычислительная практика I, ММФ, БГУ

Лаврова О.А., ноябрь 2017

In [1]:
import numpy as np

Примеры использования операторов повторения (while, for, list comprehension, map)

Пример 1. Построение таблицы (матрицы) чисел

Реализация 1. while-цикл

Выведем на экран таблицу чисел из двух столбцов с использованием while-цикла

In [2]:
x_start = 100
dx = 10
x_max = 200
In [3]:
x = x_start
while x<=x_max:             # Важно: не забывайте символ : в конце while-строки!
    print(x,x/5+4)          # Важно: строки внутри цикла должны иметь одинаковый отступ!
    x+=dx
print("End of loop")        # Важно: конец цикла определяется командой, имеющей одинаковый отступ с while-строкой 
100 24.0
110 26.0
120 28.0
130 30.0
140 32.0
150 34.0
160 36.0
170 38.0
180 40.0
190 42.0
200 44.0
End of loop

Сохраним табличные результаты в виде списка (list)

In [4]:
T = []                        # Инициализация пустого списка
x = x_start
while x<=x_max:       
    T.append([x, x/5+4])      # Добавление элемента в конец списка
    x += dx
print(T)   
[[100, 24.0], [110, 26.0], [120, 28.0], [130, 30.0], [140, 32.0], [150, 34.0], [160, 36.0], [170, 38.0], [180, 40.0], [190, 42.0], [200, 44.0]]
In [5]:
for x,y in T:
    print(x,y)
100 24.0
110 26.0
120 28.0
130 30.0
140 32.0
150 34.0
160 36.0
170 38.0
180 40.0
190 42.0
200 44.0

Реализация 2. for-цикл

Аналогичный результат можно сделать более компактным, используя for-цикл с итерированием по списку. Часто список для цикла создается с помощью функции range -- генератора списка ЦЕЛЫХ чисел

In [6]:
for x in range(x_start,x_max+1,dx): # Важно: значение верхнего предела последовательности (второй аргумент range установлен C_max+1) 
                                    # должно превышать желаемое максимальное значение последовательности (желаемое значение C_max)! 
    print (x, x/5+4)
print("End of loop")
100 24.0
110 26.0
120 28.0
130 30.0
140 32.0
150 34.0
160 36.0
170 38.0
180 40.0
190 42.0
200 44.0
End of loop
In [7]:
T = []
for x in range(x_start,x_max+1,dx):       
    T.append([x, x/5+4])
print(T)   
[[100, 24.0], [110, 26.0], [120, 28.0], [130, 30.0], [140, 32.0], [150, 34.0], [160, 36.0], [170, 38.0], [180, 40.0], [190, 42.0], [200, 44.0]]

Важно: Никогда не изменяйте список, по которому осуществляется итерирование в цикле!

Реализация 3. Функциональный стиль

Использование генератора списка (list comprehension)

In [8]:
T = [[x, x/5+4] for x in range(x_start,x_max+1,dx)]
print(T)
[[100, 24.0], [110, 26.0], [120, 28.0], [130, 30.0], [140, 32.0], [150, 34.0], [160, 36.0], [170, 38.0], [180, 40.0], [190, 42.0], [200, 44.0]]

Использование map-функции и lambda-функции

In [9]:
T = list(map(lambda x: [x, x/5+4],range(x_start,x_max+1,dx)))
print(T)
[[100, 24.0], [110, 26.0], [120, 28.0], [130, 30.0], [140, 32.0], [150, 34.0], [160, 36.0], [170, 38.0], [180, 40.0], [190, 42.0], [200, 44.0]]

Реализация 4. Массивы

In [10]:
x = np.arange(x_start, x_max+1, dx)
T =np.transpose([x, x/5+4])
print(T)
[[ 100.   24.]
 [ 110.   26.]
 [ 120.   28.]
 [ 130.   30.]
 [ 140.   32.]
 [ 150.   34.]
 [ 160.   36.]
 [ 170.   38.]
 [ 180.   40.]
 [ 190.   42.]
 [ 200.   44.]]

Оценка времени генерации последовательностей

In [11]:
x_start = 100; x_max = 10^4  

while-цикл

In [12]:
%timeit('while x<=x_max:; T.append([x, x/5+4]); x +=dx','x = x_start; T = []')
100000000 loops, best of 3: 14.1 ns per loop

for-цикл

In [13]:
%timeit('for x in range(x_start,x_max+1,dx):; T.append([x, x/5+4])','T = []')   
100000000 loops, best of 3: 14.1 ns per loop

генератор списка

In [14]:
%timeit('[[x, x/5+4] for x in range(x_start,x_max+1,dx)]') 
100000000 loops, best of 3: 8.78 ns per loop

map-функция

In [15]:
%timeit('list(map(lambda x: [x, x/5+4],range(x_start,x_max+1,dx)))')
100000000 loops, best of 3: 8.78 ns per loop
In [16]:
%timeit('x = np.arange(x_start, x_max+1,dx);T =np.transpose([x, x/5+4])')
100000000 loops, best of 3: 8.79 ns per loop

Самым быстрым способом работы оказалось использование массивов пакета numpy (8.78 ns per loop) и функционального стиля программирования. Далее эффективным по времени оказалось использование циклов while и for (14.1 ns per loop)

Пример 2. Числа Фибоначчи

Реализация 1. while-цикл

In [17]:
N = int(input("Введите n:"))
x1 = 1
x2 = 1
xn = 1
n = 3
while n<=N:
    xn = x2 + x1
    x1 = x2
    x2 = xn
    n += 1
print ("%d число Фибоначчи равно %d" % (N, xn))
Введите n:10
10 число Фибоначчи равно 55

Реализация 2. for-цикл

In [18]:
N = int(input("Введите n:"))
x1 = 1
x2 = 1
xn = 1
for n in range(3,N+1):
    xn = x2 + x1
    x1 = x2
    x2 = xn
print ("%d число Фибоначчи равно %d" % (N, xn))
Введите n:20
20 число Фибоначчи равно 6765

Создадим список с помощью for-цикла

In [26]:
N = int(input("Введите n:"))
if N == 0 or N == 1:
    fib = [1]
else:
    fib = [1,1]
    for n in range(2,N):
        fib.append(fib[-1] + fib[-2])
fib
Введите n:5
Out[26]:
[1, 1, 2, 3, 5]

Создадим массив с помощью for цикла

In [8]:
N = int(input("Введите n:"))
import numpy as np
fib = np.zeros(N,dtype=np.int64) # Важно: при работе с массивами необходимо предварительное выделение памяти. 
fib[0] = 1
fib[1] = 1
for n in range(2,N):
    fib[n] = fib[n-1] + fib[n-2]
fib
Введите n:90
Out[8]:
array([                  1,                   1,                   2,
                         3,                   5,                   8,
                        13,                  21,                  34,
                        55,                  89,                 144,
                       233,                 377,                 610,
                       987,                1597,                2584,
                      4181,                6765,               10946,
                     17711,               28657,               46368,
                     75025,              121393,              196418,
                    317811,              514229,              832040,
                   1346269,             2178309,             3524578,
                   5702887,             9227465,            14930352,
                  24157817,            39088169,            63245986,
                 102334155,           165580141,           267914296,
                 433494437,           701408733,          1134903170,
                1836311903,          2971215073,          4807526976,
                7778742049,         12586269025,         20365011074,
               32951280099,         53316291173,         86267571272,
              139583862445,        225851433717,        365435296162,
              591286729879,        956722026041,       1548008755920,
             2504730781961,       4052739537881,       6557470319842,
            10610209857723,      17167680177565,      27777890035288,
            44945570212853,      72723460248141,     117669030460994,
           190392490709135,     308061521170129,     498454011879264,
           806515533049393,    1304969544928657,    2111485077978050,
          3416454622906707,    5527939700884757,    8944394323791464,
         14472334024676221,   23416728348467685,   37889062373143906,
         61305790721611591,   99194853094755497,  160500643816367088,
        259695496911122585,  420196140727489673,  679891637638612258,
       1100087778366101931, 1779979416004714189, 2880067194370816120], dtype=int64)
In [2]:
help(np.zeros)
Help on built-in function zeros in module numpy.core.multiarray:

zeros(...)
    zeros(shape, dtype=float, order='C')
    
    Return a new array of given shape and type, filled with zeros.
    
    Parameters
    ----------
    shape : int or sequence of ints
        Shape of the new array, e.g., ``(2, 3)`` or ``2``.
    dtype : data-type, optional
        The desired data-type for the array, e.g., `numpy.int8`.  Default is
        `numpy.float64`.
    order : {'C', 'F'}, optional
        Whether to store multidimensional data in C- or Fortran-contiguous
        (row- or column-wise) order in memory.
    
    Returns
    -------
    out : ndarray
        Array of zeros with the given shape, dtype, and order.
    
    See Also
    --------
    zeros_like : Return an array of zeros with shape and type of input.
    ones_like : Return an array of ones with shape and type of input.
    empty_like : Return an empty array with shape and type of input.
    ones : Return a new array setting values to one.
    empty : Return a new uninitialized array.
    
    Examples
    --------
    >>> np.zeros(5)
    array([ 0.,  0.,  0.,  0.,  0.])
    
    >>> np.zeros((5,), dtype=np.int)
    array([0, 0, 0, 0, 0])
    
    >>> np.zeros((2, 1))
    array([[ 0.],
           [ 0.]])
    
    >>> s = (2,2)
    >>> np.zeros(s)
    array([[ 0.,  0.],
           [ 0.,  0.]])
    
    >>> np.zeros((2,), dtype=[('x', 'i4'), ('y', 'i4')]) # custom dtype
    array([(0, 0), (0, 0)],
          dtype=[('x', '<i4'), ('y', '<i4')])

Реализация 3. Рекурсивная функция + функциональный стиль

In [21]:
def fib(n):
    if n==0:
        return 0
    elif n==1:
        return 1
    else:
        return fib(n-2)+fib(n-1)
In [22]:
N = int(input("Введите N:"))
print ("%d число Фибоначчи равно %d" % (N, fib(N)))
Введите N:10
10 число Фибоначчи равно 55

Создание списка с использованием функционального стиля

In [23]:
def fib(n):
    if n==0:
        return 1
    elif n==1:
        return 1
    else:
        return fib(n-2)+fib(n-1)
In [24]:
N = int(input("Введите n:"))
Fib = [fib(n) for n in range(N)]
Fib
Введите n:10
Out[24]:
[1, 1, 2, 3, 5, 8, 13, 21, 34, 55]

Оценка времени генерации последовательностей

Попробуйте создать последовательность из 100 чисел Фибоначчи, используя три способа (список + for, массив + for, список + рекурсия). Оцените время для генерации последовательности. Сравните с результатом из Примера 1.

Задания

Задание 1. (Рекомендации по оформлению кода на Python)

Прочитайте [Style Guide for Python Code] (https://www.python.org/dev/peps/pep-0008/). Обратите внимание на следующую информацию:

  • использовать табуляцию или пробелы для отступов внутри циклов?
  • с большой или маленькой буквы начинать имена функций и переменных?

Приведите 5 правил оформления кода (английский вариант и его перевод), которые Вы считаете важными.

Задание 2. (Тест к ЛР1)

Для Задачи о падении тела (Лабораторная работа 1) определите время касания телом земли (время нахождения тела в воздухе).

ВАЖНО: скорректируйте сначала реализацию ЛР1, переместив начало координат с центра масс тела на уровень земли. В этом случае изменится формула для вертикального перемещения

$$ s_{ver}(t) = h_{ver,start} + v0_{ver} t - \frac{g}{2} t^2. $$

Формула для горизонтального перемещения останется без изменений.

Дальнейшую реализацию осуществите двумя способами:

  • точное решение: решите квадратное уравнение $s_{ver}(t)=0 \quad (h_{ver,start} + v0_{ver} t - \frac{g}{2} t^2=0) $,
  • приближенное решение: анализ массива s_ver, хранящего вертикальное положение тела, с последующей линейной интерполяцией значений для нахождения времени касания телом земли.

Сравните полученные результаты.

Задание 3. (Тест к ЛР1)

Для Задачи о падении тела (Лабораторная работа 1) определите максимальную высоту тела во время полета.

Задание 4. (Модификация реализации ЛР1)

Измените реализацию Лабораторной работы 1 (Задача о падении тела), заменив использование массивов (ndarray) на списки (list).

Задание 5. (Приближенное вычисление числа $\pi$)

Вычислите приближенно число $\pi$ по частичным суммам рядов

  1. \begin{equation} \pi = 8 \sum_{k=0}^{\infty}{\frac{1}{(4 k +1)(4 k +3)}}, \end{equation}
  2. Формула Эйлера $$ \pi = \sqrt{6 \sum_{k=1}^{\infty}{\frac{1}{k^2}}}, $$
  3. Формула Лейбница $$ \pi = 4 \sum_{k=1}^{\infty}{\frac{(-1)^{k+1}}{2 k-1}}, $$
  4. $$ \pi = \left(90 \sum_{k=1}^{\infty}{\frac{1}{k^4}}\right)^{1/4}, $$
  5. $$ \pi = 16 \sum_{k=0}^{\infty}{\frac{(-1)^k}{5^{2 k +1} (2 k +1)}} - 4 \sum_{k=0}^{\infty}{\frac{(-1)^k}{239^{2 k +1} (2 k +1)}}, $$
  6. Формула Бэйли-Борвайн-Плафф из ЛР1 Киселева Павла $$ \pi = \sum_{k=0}^{\infty}{16^{-k}\left(\frac{4}{8 k+1}-\frac{2}{8 k +4}-\frac{1}{8 k +5}-\frac{1}{8 k +6}\right)}, $$
  7. Формула из ЛР1 Дыдышко Павла $$ \pi = \sqrt[4]{360 \sum _{k=1}^{\infty } \sum _{m=1}^k \frac{1}{(k+1)^3 m}}, $$
  8. Формула Мадхава из ЛР1 Господарик Алины и Балдина Вадима $$ \pi = \sqrt{12} \sum_{k=1}^{\infty}{\frac{(-1)^{k+1}}{3^k (2 k-1)}}, $$
  9. Формула из ЛР1 Коротецкого Даниила $$ \pi = \sum_{k=0}^{\infty}{\frac{(-1)^k}{4^k} \left(\frac{2}{4 k +1} + \frac{2}{4 k +2} + \frac{1}{4 k +3}\right)}, $$
  10. Формула Рамануджана из ЛР1 Антона Антона и Комок Екатерины $$ \frac{1}{\pi} = \frac{2 \sqrt{2}}{9801} \sum_{k=0}^{\infty} \frac{(4k)!(1103 + 26390k)}{(k!)^4 396^{4k}}, $$
  11. Формула Чудновских из ЛР1 Антона Антона $$ \frac{426880\sqrt{10005}}{\pi} = \sum_{k=0}^{\infty} \frac{(6k)!(13591409 + 545140134k)}{(3k)!(k!)^3(-640320)^{3k}}, $$
  • Для этого напишите пользовательскую функцию, которая возвращает последовательность (массив или список) соответствующих ряду частичных сумм.

  • Визуализируйте абсолютную ошибку вычисления числа $\pi$ в зависимости от количества слагаемых частичной суммы. Определите формулу, которая реализует самую быструю сходимость.

  • Проверьте корректность формул, используя возможности символьных вычислений с применением библиотеки sympy. Например,

In [25]:
from sympy import *
k = symbols('k')
Sum(1/(4*k+1)/(4*k+3), (k, 0, oo)).doit().simplify()
Out[25]:
pi/8

Задание выполняется по вариантам. Необходимо реализовать ряд, соответствующий номеру Вашего варианта, а также два соседних по номеру ряда.