Ctrl+K

Домашнее задание

Домашнее задание#

Задача 1. Вводное задание.#

Вектора х1 и у1 из архива дают значения некоторой функции \(f(x) .\) Постройте на графике, используя библиотечные функции:

(a) Интерполяционный Лагранжев многочлен.

(b) Интерполяционный кубический сплайн (CubicSpline).

(c) Монотонный кубический интерполянт (PchipInterpolator).

(d) Аппроксимацию полиномами, используя функции np. polyfit и np.polyval.

Чтобы прочитать из архива данные, используйте следующий код:

import numpy as np
with np.load('data_interp.npz') as data:
  x1, y1 = data['x1'], data['y1']

Задача 2. Интерполяция полиномом Лагранжа.#

Допишите класс, который конструирует интерполяционный полином Лагранжа, проходящий через точки, заданные как xk и yk.

import numpy as np

class LagrangeInterpolator:
    """Lagrange interpolating polynomial.

    Given a set of pairs ``(x_k, y_k)``, construct
    a Lagrange polynomial ``f(x)``, such that

    .. math::

        f(x_k) = y_k   for k =0, ..., n-1

    Parameters
    ----------
    xk : array_like, shape(n,)
        Abscissas
    yk : array_like, shape(n,)
        Ordinates

    Attributes
    ----------
    __call__

    """
    def __init__(self, xk, yk):
        self.xk = np.asarray(xk, dtype=float)
        self.yk = np.asarray(yk, dtype=float)

    def __call__(self, x):
        """Evaluate the interpolator at a given point.

        Parameters
        ----------
        x : float

        Returns
        -------
        the value of the interpolator at ``x``.
        """
        # YOUR CODE HERE
        raise NotImplementedError()

def runge_func(x, a=25):
    return 1.0 / (1.0 + a*x**2)

xx = np.linspace(-2, 2, 21)
yy = runge_func(xx)

lagr = LagrangeInterpolator(xx, yy)

from numpy.testing import assert_allclose

assert_allclose(yy,
                [lagr(xval) for xval in xx],
                atol=1e-14)

Задача 3. Феномен Рунге#

Рассмотрим функцию Рунге, \(y = 1/(1 + 25x^2)\). Интерполируйте эту функцию на интервале \(x\in [-2, 2]\), используя полином Лагранжа с \(m=3, 5, 7, 11\). Используйте равномерную сетку. Нарисуйте результат интерполяции вместе с исходной функцией на одном графике.

# YOUR CODE AND COMMENTS HERE

Допишите функцию, которая возвращает узлы Чебышева.

def cheb_nodes(n, a=-1, b=1):
    r"""Chebyshev nodes of degree $n$ on $[a, b]$
    """
    # YOUR CODE HERE
    raise NotImplementedError()

nodes_11 = cheb_nodes(11)
nodes_11 = np.asarray(nodes_11)
assert (nodes_11[1:] > nodes_11[:-1]).all()

from scipy.special import roots_chebyt
nodes, weights = roots_chebyt(5)

assert_allclose(cheb_nodes(5),
                nodes, atol=1e-14)

assert_allclose(cheb_nodes(5, a=-1, b=3),
                nodes*2 + 1, atol=1e-14)

Повторите интерполяцию Лагранжа функции Рунге с помощью узлов Чебышева. Постройте интерполянты. Также постройте интерполяцию кубическим сплайном тех же данных (scipy.interpolate.CubicSpline). Сравните величину явления Рунге для равномерной сетки и сетки Чебышева. Демонстрирует ли интерполяция сплайнами феномен Рунге?

# YOUR CODE AND COMMENTS HERE

Другим известным трудным тестом для интерполяции является следующая периодическая функция:

\(u(x)=\frac{\sqrt{\varepsilon(2+\varepsilon)}}{2 \pi(1+\varepsilon-\cos x)}, \quad-\pi \leq x \leq \pi, \quad \varepsilon=0.21 \)

Проведите аналогичные вычисления и постройте графики для неё.

# YOUR CODE AND COMMENTS HERE

Задача 4. UnivariateSpline.#

Прочитать документацию к UnivariateSpline и определить оптимальное значение для параметра сглаживания s, учитывая, что шум был задан в виде \(u_i = \dots + 0.5 \xi_i, \quad \xi_i \sim \mathcal N(0, 1)\). Экспериментально проверить оптимальность этого значения для какой-либо функции

Задача 5. Производная и шум#

Цель - найти производную следующей таблично-заданной функции с шумом. Всего есть 1000 узлов. Генерировать новые точки нельзя.

\[ u(t)=\frac{200}{1+\frac{t}{200}}\left(0.5+0.5 \cdot \cos \left(0.04 \cdot t\right)\right)+\text{noise} \]

  1. Сделайте это через конечные разности, используя минимальное даное расстояние между точками. Вспоминаем семинар 1.

  2. Сделайте это через оптимальный шаг численного дифференциирования. То есть, при оценке производной надо брать несоседние точки исходя из погрешностей \(u_i\). Какой шаг надо брать по индексам точек?

  3. Сделайте это через UnivariateSpline, используя метод .derivative() - он аналитически дифференциирует объект сплайна и выдаёт также объект сплайна (т.е. можем обращаться к нему как к функции). Не забывайте корректно определить параметр s.

Сравните все 3 способа численного взятия производной с аналитической производной как графически, так и используя какой-нибудь Loss.

import numpy as np

t = np.linspace(0, 1000, 1000)
u = 200 / (1 + t/200) * (0.5 + 0.5 * np.cos(0.04 * t)) + 0.5 * np.random.randn(len(t))

Задача 6*. Аппроксимация узлов Чебышева.#

Вектора х5 и у5 из архива задают некоторую функцию на равномерной сетке на отрезке \([-1,1] .\)

Составьте соответствующий интерполяционный многочлен и постройте его график. Ожидаемо, получившийся интерполянт демонстрирует явление Рунге. Теперь проделайте замену переменной \(x \rightarrow s\) по формуле \(x=g(s)\), где

\[ g(s)=\frac{1}{53089}\left(40320 s+6720 s^{3}+3024 s^{5}+1800 s^{7}+1225 s^{9}\right) \]

функция, приближённо отображающая чебышевскую сетку в равномерную. Да, вам надо численно реализовать обратное отображение. В новой переменной \(s\), составьте интерполяционный многочлен \(P(s) .\) Постройте график получившейся интерполяции исходной функции \(P(s(x))\). Делает ли такая процедура явление Рунге менее выраженным?

Содержание