MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Spline Toolbox

Задание граничных условий в функции csape

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

Пусть для - точки разрыва, для - значения интерполируемой функции. Функция csape допускает задание следующих стандартных граничных условий для сплайна , кроме упомянутого:

1) задать или ;

2) задать или , в том числе и для получения естественного сплайна;

3) задать периодические условия: , ;

4) задать условия отсутствия узла (not-a-knot condition), которые соответствуют непрерывности в точках и , т.е. на отрезках и полином один и тот же, соответственно на отрезках и полином тоже один;

5) задать более сложные граничные условия вида: или ;

Для задания граничных условий есть несколько способов. Простые граничные условия указываются при помощи дополнительного входного аргумента - строки:

  • pp = csape(ksi, y, 'periodic') - построение сплайна в кусочно-полиномиальной форме с периодическими граничными условиями;
  • pp = csape(ksi, y, 'variational') - построение естественного сплайна в кусочно-полиномиальной форме;
  • pp = csape(ksi, y, 'complete') или pp = csape(ksi, y, 'clamped') - то же, что и по умолчанию pp = csape(ksi, y), т.е. построение сплайна, первая производная которого в ) совпадает с производной интерполяционного полинома третьей степени, построенного по ) (см. предыдущий раздел).
  • pp = csape(ksi, y, 'not-a-knot') - построение сплайна с условиями отсутствия узла.

    Следующие примеры демонстрируют построение сплайнов с заданными граничными условиями и вычисление производных и скачков при помощи функций fnder и fnjmp для проверки получающихся сплайнов.

    Задание точек разрыва и значений функции:
    >> ksi = 0:pi/10:2*pi;
    >> L = length(ksi)-1;
    >> y = sin(ksi);
    
    Пример периодического сплайна:
    >> pp = csape(ksi, y, 'periodic');
    >> pp1 = fnder(pp);
    >> pp2 = fnder(pp, 2);
    >> fnjmp(pp1, ksi(1)) - fnjmp(pp1, ksi(L+1))
    ans =
         0
    >> fnjmp(pp2, ksi(1)) - fnjmp(pp2, ksi(L+1))
    ans =
         0
    
    Пример сплайна с условиями отсутствия узла:
    >> pp = csape(ksi, y, 'not-a-knot');
    >> pp3 = fnder(pp, 3);
    >> fnjmp(pp3, ksi(2))
    ans =
     -2.1316e-014
    >> fnjmp(pp3, ksi(L))
    ans =
      1.3323e-015
    
    Пример естественного сплайна 
    >> pp = csape(ksi, y, 'variational');
    >> pp2 = fnder(pp, 2);
    >> fnjmp(pp2, ksi(1))
    ans =
         0
    >> fnjmp(pp2, ksi(L+1))
    ans =
         0 
    

    Третьим входным аргументом функции csape может быть не только текстовая строка, но и вектор conds из двух элементов, содержащий числа 0, 1 или 2:

    pp = csape(ksi, y, conds)

    Такой способ предоставляет более широкие возможности для постановки граничных условий: во-первых, граничные условия могут быть разных типов и зависят от значений conds(1) и conds(2). Первый элемент conds(1) соответствует условию на левой границе отрезка, а conds(2) - на правой. Во-вторых, можно задавать значения первой или второй производных сплайна . Эти дополнительные значения должны быть указаны в первом и последнем элементе вектора - второго входного аргумента функции csape. Его элементы со второго по предпоследний будут интерпретироваться как значения интерполируемой функции.

    Будем считать, что вектор y содержит значения интерполируемой функции: y(i) - значение функции в ksi(i). Рассмотрим смысл значений элементов вектора conds

    conds = [0 0] - периодические граничные условия, при этом дополнительных элементов в векторе y не требуется (то же самое, что pp = csape(ksi, y, 'periodic')).

    conds = [1 1] - на обеих границах отрезка заданы значения первой производной сплайна. Возможны два случая:

  • pp = csape(ksi, [A y B], [1 1]) - считается, что и ;
  • pp = csape(ksi, y, [1 1]) - по умолчанию значение в точке (соответственно ) совпадает со значением первой производной интерполяционного полинома в этой точке, построенного по узлам (соответственно ).

    conds = [2 2] - на обеих границах отрезка заданы значения второй производной сплайна. Возможны два случая:

  • pp = csape(ksi, [A y B], [2 2]) - считается, что и ;
  • pp = csape(ksi, y, [2 2]) - по умолчанию .

    conds = [1 2] - на левой границе отрезка задано значение первой производной сплайна, на правой - второй. Возможны два случая:

  • pp = csape(ksi, [A y B], [1 2]) - считается, что и ;
  • pp = csape(ksi, y, [1 2]) - по умолчанию совпадает со значением первой производной интерполяционного полинома в этой точке, построенного по узлам , а .

    conds = [2 1] - аналогично предыдущему случаю, только на левой границе отрезка задано значение второй производной сплайна, на правой - первой.

    Рассмотрим теперь случай нестандартных граничных условий. В справочной системе по Spline Toolbox приведен пример построения сплайна при помощи csape, удовлетворяющего условию:

    для заданных a, b и c. Идея заключается в отыскании сплайна в виде , где - интерполяционный сплайн для исходных данных, - некоторый сплайн, значения которого в точках равны нулю, а коэффициент выбирается так, чтобы выполнялось заданное граничное условие.

    Несложно предложить простой подход, позволяющий построить сплайн, который удовлетворяет двум граничным условиям следующего вида:

    Будем искать сплайн в виде линейной комбинации , где - интерполяционный сплайн для исходных данных, и - некоторые сплайны, значения которых в точках равны нулю (см. пример на рис. 1, где точки разрыва отмечены круглыми маркерами).

    Рис. 1. Пример сплайнов для

    Коэффициенты и выберем, исходя из заданных граничных условий. Очевидно, для нахождения и получается система из двух уравнений:

    Преобразуем ее и получаем:

    Если ее определитель не равен нулю, то нужный сплайн построен.

    Рассмотрим пример построения сплайна , удовлетворяющего граничным условиям

    для данных:

    , , ,

    y = [0, 0.1, 0.4, 1, 0.3, -0.1, -1, 0.5, 0.1, 0.5, 1].

    % Задаем данные
    ksi = 0:pi/5:2*pi;
    L = length(ksi)-1;
    y = [0 0.1 0.4 1 0.3 -0.1 -1 --0.5 0.1 0.5 1];
    
    % Вводим коэффициенты граничных условий
    a1 = 1.2;    b1 = 0.4;   c1 = 0.7;
    a2 = -0.3;   b2 = 1.8;   c2 = -3.9;
    
    % Строим сплайн для аппроксимации данных с граничными условиями по умолчанию
    % (производная сплайна в граничных точках совпадает с производной 
    % кубического интерполяционного полинома, проходящего через граничную точку 
    % и три точки, расположенные рядом)
    pp1 = csape(ksi, y);
    
    % Строим два вспомогательных сплайна (см. выше про задание граничных условий)
    pp2 = csape(ksi, [1 zeros(size(y)) 0], [1 1]);
    pp3 = csape(ksi, [0 zeros(size(y)) 1], [1 1]);
    
    % Находим первую и вторую производные от сплайнов pp1, pp2, pp3 
    % (здесь можно сделать более эффективно, см. примечание после программы)
    pp11 = fnder(pp1, 1);
    pp12 = fnder(pp1, 2);
    pp21 = fnder(pp2, 1);
    pp22 = fnder(pp2, 2);
    pp31 = fnder(pp3, 1);
    pp32 = fnder(pp3, 2);
    
    % Вычисляем элементы матрицы системы
    M = zeros(2);
    M(1, 1) = a1*fnval(pp21, ksi(1)) + b1*fnval(pp22, ksi(1));
    M(1, 2) = a1*fnval(pp31, ksi(1)) + b1*fnval(pp32, ksi(1));
    M(2, 1) = a2*fnval(pp21, ksi(L+1)) + b2*fnval(pp22, ksi(L+1));
    M(2, 2) =a2*fnval(pp31, ksi(L+1)) + b2*fnval(pp32, ksi(L+1));
    
    % Вычисляем элементы вектора правой части системы
    c = zeros(2,1);
    c(1) = c1 - a1*fnval(pp11, ksi(1)) - b1*fnval(pp12, ksi(1));
    c(2) = c2 - a2*fnval(pp11, ksi(L+1)) - b2*fnval(pp12, ksi(L+1));
    
    % Решаем систему
    d = M\c;
    
    % Берем линейную комбинацию pp1, pp2 и pp3, используя функцию fncmb 
    % для выполнения арифметических операций над сплайнами
    % Записываем искомый сплайн в P.
    P = fncmb(pp2, d(1), pp3, d(2));
    P = fncmb(pp1, '+', P);
    
    % Проверяем, удовлетворяются ли граничные условия
    P1 = fnder(P, 1);
    P2 = fnder(P, 2);
    a1*fnval(P1, ksi(1)) + b1*fnval(P2, ksi(1)) - c1
    a2*fnval(P1, ksi(L+1)) + b2*fnval(P2, ksi(L+1)) - c2
    
    Примечание.
    Не обязательно считать производные для всех сплайнов pp1, pp2 и pp3. Достаточно было выделить первую и
    последнюю полиномиальную части при помощи функции fnbrk и найти производные только для них:
    
    ...
    % Выделяем первую и последнюю полиномиальные части в pp1, pp2, pp3
    pp1L = fnbrk(pp1, 1);
    pp1R = fnbrk(pp1, L);
    pp2L = fnbrk(pp2, 1);
    pp2R = fnbrk(pp2, L);
    pp3L = fnbrk(pp3, 1);
    pp3R = fnbrk(pp3, L);
    
    % Находим первую и вторую производные этих частей. 
    pp11L = fnder(pp1L,1);
    pp12L = fnder(pp1L,2);
    pp11R = fnder(pp1R,1);
    pp12R = fnder(pp1R,2);
    
    pp21L = fnder(pp2L,1);
    pp22L = fnder(pp2L,2);
    pp21R = fnder(pp2R,1);
    pp22R = fnder(pp2R,2);
    
    pp31L = fnder(pp3L,1);
    pp32L = fnder(pp3L,2);
    pp31R = fnder(pp3R,1);
    pp32R = fnder(pp3R,2);
     
    % Формируем матрицу системы уравнений
    M=zeros(2);
    M(1,1) = a1*fnval(pp21L, ksi(1)) + b1*fnval(pp22L, ksi(1));
    M(1,2) = a1*fnval(pp31L, ksi(1)) + b1*fnval(pp32L, ksi(1));
    M(2,1) = a2*fnval(pp21R, ksi(L+1)) + b2*fnval(pp22R, ksi(L+1));
    M(2,2) = a2*fnval(pp31R, ksi(L+1)) + b2*fnval(pp32R, ksi(L+1));
    
    % Формируем вектор системы уравнений
    c = zeros(2, 1);
    c(1) = c1 - a1*fnval(pp11L, ksi(1)) - b1*fnval(pp12L, ksi(1));
    c(2) = c2 - a2*fnval(pp11R, ksi(L+1)) - b2*fnval(pp12R, ksi(L+1));
    ...
    Все остальное аналогично.
    

    Функция csapi для интерполяции с условиями отсутствия узла (not-a-knot)

    Напомним, что условия отсутствия узла (not-a-knot condition) заключаются в непрерывности третьей производной сплайна в точках и , т.е. на отрезках и полином один и тот же. Соответственно, на отрезках и полином тоже один. Поскольку интерполяция с условиями отсутствия узла достаточно широко распространена, то для нее предусмотрена отдельная функция csapi, которая строит интерполяционный сплайн в кусочно-полиномиальной форме. Обращения

    pp = csapi(ksi, y)    и      pp = csape(ksi, y, 'not-a-knot')
    
    приводят к одинаковому результату (рассмотренная выше функция csape вызывает csapi, если заданы граничные условия отсутствия узла). Проинтерполируем функцию при помощи csape и csapi, а затем найдем разность полученных сплайнов и убедимся, что она тождественно равна нулю:

    ksi = 0:0.1:10;
    y = sin(x);
    % интерполируем функцией csapi
    ppi = csapi(ksi, y);
    % интерполируем функцией csape с условиями отсутствия узла
    ppe = csape(ksi, y, 'not-a-knot');
    % находим разность сплайнов и сравниваем коэффициенты
    pp0 = fncmb(ppi, '-', ppe);
    max(max(abs(pp0.coefs)))
    
    ans =
         0
    
    Вызов функции csapi аналогичен следующему обращению к функции spline, входящей в MATLAB:
    ppi = csapi(ksi, y);
    pps = spline(ksi, y);
    
    и приводит к тому же самому сплайну (с точностью до ошибок округления при вычислениях):
    pp0 = fncmb(ppi, '-', pps);
    max(max(abs(pp0.coefs)))
    ans =
      1.1519e-014
    
    Интерфейс функций csapi и spline одинаков, они умеют не только выводить структуру с информацией о кусочно-полиномиальной форме сплайна, но и значения сплайна в заданных точках. Для нахождения значений сплайна следует вызвать csapi с третьим входным аргументом - массивом точек xx:
    yy = csapi(ksi, y, xx)
    
    Иногда вычисление значений сплайна в достаточно большом числе точек позволяет построить его более точно, чем функция fnplot, входящая в Spline Toolbox, например:
    ksi = 0:0.1:10;
    y = sin(x.^2);
    figure
    ppi = csapi(ksi, y);
    subplot(2,1,1)
    fnplt(ppi)
    xx = ksi(1):(ksi(end)-ksi(1))/1000:ksi(end);
    yy = csapi(ksi, y, xx);
    subplot(2,1,2)
    plot(xx, yy)
    
    приводят к двум графикам одного и того же сплайна, нижний из которых более точен.



    Примечание.
    Для вычисления сплайна в точках xx не обязательно было второй раз вызывать csapi и заново конструировать сплайн. Можно воспользоваться функцией fnval, которая вычисляет значения сплайна в заданных точках. В приведенном выше примере построение точного графика сплайна может быть реализовано при помощи следующих команд
    figure
    ppi = csapi(ksi, y);
    yy = fnval(ppi, xx);
    plot(xx, yy)
    
    Аналогичным образом вычисляются значения сплайнов, построенных другими функциями Spline Toolbox.


  • Поиск по сайту:

    Система Orphus

    Яндекс.Метрика