Введение

Александр Теофил Вандермонд (28 февраля 1735 - 1 января 1796) - французский музыкант и математик, известный благодаря своей работе в области высшей алгебры.

Главным увлечением Вандермонда длительное время была лишь музыка, но к 35-ти годам юный ученый обратился к математике. Первым делом он провел исследование симметрических функций и решения круговых полиномов, после чего выпустил три статьи: про задачу о ходе коня, про комбинаторику и про основы теории детерминантов. 

В честь Александра Теофила был назван специальный класс матриц - матрицы Вандермонда, о котором пойдет речь в данной статье. [1]

Матрица Вандермонда

В линейной алгебре матрица Вандермонда, представляет собой матрицу с членами геометрической прогрессии в каждой строке. 

Давайте рассмотрим узлы x_i, где x_i =\frac{i}{n} для i = 0, 1, …, n. 

Тогда первый столбец матрицы будет полностью заполнен 1, оно же x_i^0; второй столбец будет заполнен x_i, третий - x_i^2 и так далее. Общий вид матрицы:

\begin{bmatrix} 1 & x_0 & x_0^2 & ... & x_0^n\\ 1 & x_1 & x_1^2 & ...& x_1^n \\ 1 & x_2 & x_2^2 & ... & x_2^n\\ ... & ... & ... &... &... \\ ... &... & ... &... &x_n^n \end{bmatrix}

Применение - интерполяция

Задача интерполяции заключается в том, чтобы найти полином вида:

p(x)=a_0 + a_1*x + a_2*x^2 + … + a_n*x^n

который удовлетворяет условию:

p(x_0)=y_0, …, p(x_m)=y_m. 

Сделать это можно используя матрицу Вандермонда, следуя формуле:

V * a = y

где V - матрица Вандермонда, a = (a_0, …, a_n) - вектор коэффициентов и y = (y_0, …, y_n) = (p(x_0), …, p(x_n)) - вектор значений. [2]

\begin{bmatrix} 1 & x_0 & x_0^2 & ... & x_0^n\\ 1 & x_1 & x_1^2 & ...& x_1^n \\ 1 & x_2 & x_2^2 & ... & x_2^n\\ ... & ... & ... &... &... \\ ... &... & ... &... &x_n^n \end{bmatrix} \begin{bmatrix} a_0 \\ a_1  \\ ... \\ a_n  \end{bmatrix} =\begin{bmatrix} y_0 \\ y_1  \\ ... \\ y_n  \end{bmatrix}

Пример

Даны три точки: (1, 2); (3, 4); (2, 3), нужно найти полином вида

p(x)=a_0 + a_1x + a_2x^2 + … + a_n*x^n, проходящий через эти точки.

Для начала, построим матрицу Вандермонда. x_0=1, x_1=3, x_2=2

\begin{bmatrix} 1 & 1 & 1 \\ 1 & 3 & 9 \\ 1 & 2 & 4 \end{bmatrix}

Далее, расширим матрицу значениями вектора y

\begin{bmatrix} 1 & 1 & 1 &| & 2 \\ 1 & 3 & 9 & | & 4 \\ 1 & 2 & 4 & | & 3 \end{bmatrix}

Пользуясь правилами сложения и вычитания строк матрицы, упростим её до формы:

\begin{bmatrix} 1 & 0 & 0 & | & 1 \\ 0 & 1 & 0 & | & 1 \\ 0 & 0 & 1 & | & 0\end{bmatrix}

Далее, решаем систему методом подстановки, a_2=0, a_1=1, a_0=1.

Таким образом, наш полином равен p(x)=1 + x

Код

Мы создадим собственную функцию vandermonde_matrix,которая вычисляет матрицу

Вандермонда (n + 1) × (n + 1) с узлами x_i, где x_i =\frac{i}{n} для i = 0, 1, …, n.

Входные данные функции — степень n искомого полинома, а выходные данные — матрица.

import numpy as np

def vandermonde_matrix(n):
    A=np.zeros((n+1,n+1)) 
    for i in range (n+1):
        for j in range (n+1): 
            A[i][j]=(i/n)**j if n !=0 else 1
    return A

print(vandermonde_matrix(3)) #пример
результат прогонки кода (матрица 4х4)
результат прогонки кода (матрица 4х4)

Число обусловленности

В численном анализе число обусловленности функции показывает насколько может измениться значение функции при небольшом изменении аргумента. Число обусловленности отражает, насколько чувствительна функция к изменениям или ошибкам на входе и насколько ошибка на выходе является результатом ошибки на входе. [3]

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

Код

Для n = 1,2,...,100, мы построим матрицу Вандермонда (n + 1) × (n + 1);

вычислим число обусловленности матрицы и построим график зависимости числа обусловленности матрицы Вандермонда от n.

import numpy as np
import matplotlib.pyplot as plt


def vandermonde_matrix(n):
    A=np.zeros((n+1,n+1)) 
    for i in range (n+1):
        for j in range (n+1): 
            A[i][j]=(i/n)**j if n !=0 else 1
    return A

cond_vandermonde=np.zeros(101)

for i in range (1, 101):
    matrix=vandermonde_matrix(i)
    matrix_inv=np.linalg.inv(matrix)
    cond_vandermonde[i]=np.array(np.linalg.norm(matrix)*np.linalg.norm(matrix_inv))
    
n_values=np.zeros(101)

for i in range(1, 101):
    n_values[i]=np.array(i)
plt.semilogy(n_values, cond_vandermonde)
plt.ylabel('n')
plt.xlabel('Число Обусловленности')
plt.title('Зависимость числа обусловленности от n') 
plt.show()
результат прогонки кода
результат прогонки кода

Мы можем заметить, что когда увеличивается размер матрицы, число обусловленности также возрастает. Это связано с тем, что матрица Вандермонда содержит всё большие и большие разности элементов, что приводит к увеличению значений определителя. 

Из-за того, что определитель матрицы Вандермонда может сильно возрастать с увеличением размера матрицы, происходят численные нестабильности и потери точности вычислений. Благодаря чему мы можем подвести итог. 

Заключение

Метод интерполяции с использованием матрицы Вандермонда эффективен для небольшого количества точек, когда нужно построить матрицу небольшого размера. Однако, для больших данных или случаев с высокой чувствительностью к точности, могут потребоваться более устойчивые методы.

Список Литературы

  1. https://ru.m.wikipedia.org/wiki/Вандермонд,АлександрТеофил

  2. https://en.m.wikipedia.org/wiki/Vandermonde_matrix

  3. https://ru.m.wikipedia.org/wiki/Число_обусловленности

Комментарии (7)


  1. wataru
    17.08.2024 18:03
    +4

    Вас не смущает график числа обусловленности? Странное поведение: линейный рост, а потом какие-то колебания. Выглядят, как ошибка в вычисленях. Возможно, вызванная переполнением переменных.

    Заодно, могли бы вы объяснить что вот этот вот код в питоне делает?

    n_values[i] = np.array(i)


    1. kristina_ponomareva Автор
      17.08.2024 18:03

      Здравствуйте!

      Нет, не смущает. Колебания вызваны свойствами матрицы Вандермонда, которая часто плохо обусловлена, особенно при больших размерах матрицы. Это приводит к резкому увеличению числа обусловленности и его колебаниям!) ошибок в вычислениях нет.

      Что касается вопроса про код, эта строка берет число i, преобразует его в формат массива и добавляется в список n_values. Сейчас, когда вы написали этот комментарий, я поняла, что можно было не делать через array. Но ошибки нет:)

      Спасибо за Ваш комментарий!


      1. wataru
        17.08.2024 18:03
        +7

        Колебания вызваны свойствами матрицы Вандермонда

        Попробуйте поменять тип данных в матрице на np.float32 (по умолчанию там float64). Посмотрите, что произойдет с вашим графиком. Он станет колебаться раньше.

        У вас переполнение. Точности вещественных числел тупо не хватает для вычисления обратной матрицы.

        Попробуйте использовать пакет sympy. Или какой-нибудь Pygmp, правда обращение матрицы придется писать руками методом Гаусса, да и вычисление нормы тоже.

        Погуглите свойства этой самой матрицы. Уверен, для числа обусловленности есть какая-нибудь аж замкнутая формула.


    1. vadimr
      17.08.2024 18:03
      +2

      Так ясен пень, что если i делить на 10**17, то результат незначим.


      1. kristina_ponomareva Автор
        17.08.2024 18:03

        Есть методы интерполяции которые работают всегда, а есть которые нет) про это и рассказываю


        1. vadimr
          17.08.2024 18:03
          +1

          "Всегда" в вещественной арифметике не бывает.


          1. kristina_ponomareva Автор
            17.08.2024 18:03

            Это да