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

Spline Toolbox

Вычисление коэффициентов кубических полиномов Эрмита и кубических сплайнов

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

Вычисление коэффициентов полиномов Эрмита

При интерполировании полиномами Эрмита требуется знать не только значения интерполируемой функции в узлах, но и значения ее производной. Кусочные полиномы строятся так, чтобы их значения совпадали со значениями интерполируемой функции в узлах и значения их первых производных так же совпадали со значениями первой производной интерполируемой функции в узлах интерполирования.

Итак, нам задано

где xk - упорядоченные по возрастанию узлы интерполирования, а ƒ(x) - интерполируемая функция. Требуется найти n - 1 кубический полином

определенный на отрезке [xk, xk+1] такой, что выполняются следующие четыре условия (как раз четыре условия нужны для определения четырех коэффициентов каждого из полиномов):

Ниже на рисунке приведен пример для пяти узлов интерполирования и, соответственно, четырех полиномов.

Поскольку

то сразу же получаем, что

Оставшиеся два условия приводят нас к системе из двух линейных алгебраических уравнений относительно двух неизвестных коэффициентов ck и dk:

Решая эту систему относительно ck и dk получим, что

Здесь введено обозначение

Эта величина называется разделенной разностью 1-го порядка функции ƒ(x).

Таким образом, выражения для коэффициентов ak, bk, ck и dk каждого из кубических полиномов Pk(x) найдены.

Для изучения поведения полиномов Эрмита может оказаться полезной функция Hpoly, приведенная ниже, которая по заданным узлам (вектор x), значениям функции (вектор y) и ее производной (вектор y1) в узлах строит график кусочного полинома и выводит рядом с каждым полиномом его коэффициенты, начиная с коэффициента при нулевой степени. Для вычисления разделенных разностей использована функция diff, которая из заданного вектора делает вектор, состоящий из разностей его соседних элементов, т.е. если, например исходный вектор x=[2 3 5 9], то diff(x) вернет 1, 2, 4.

function Hpoly(x,y,y1)
% определение числа узлов
n=length(x);
% нахождение всех коэффициентов a(k) и b(k)
a = y;
b = y1;
% вычисление разделенных разностей 
DivDif = diff(y)./diff(x);
% вычисление всех коэффициентов c(k) и d(k)
for k = 1:n-1
    c(k) = (3*DivDif(k) - 2*y1(k) - y1(k+1))/(x(k+1)-x(k));
    d(k) = (y1(k) - 2*DivDif(k) + y1(k+1))/(x(k+1)-x(k))^2;
end
% создание графического окна и осей
figure('Color','w')
axes
% вывод данных маркерами
plot(x,y,'ro')
hold on
% построение графиков каждого из полиномов
for k=1:n-1
    % вычисление k-го полинома в 30 точках между x(k) и x(k+1)
    xx=linspace(x(k),x(k+1),30);
    yy=a(k) + b(k)*(xx-x(k)) + c(k)*(xx-x(k)).^2 + d(k)*(xx-x(k)).^3;
    % выбор случайного цвета
    rgbcolor = [rand rand rand];
    % построение графика k-го полинома
    plot(xx,yy,'Color', rgbcolor) 
    % вывод коэффициентов k-го полинома
    text('Position',[xx(15) yy(15)],...
        'String',num2str([a(k); b(k); c(k); d(k)]), 'Color', rgbcolor,...
        'VerticalAlignment','middle','HorizontalAlignment','center')
end

Например, для интерполяции функции на отрезке с шагом кубическими полиномами Эрмита достаточно выполнить следующие команды

>> x=1:10;
>> y=sin(x)+sin(2*x);
>> y1=cos(x)+2*cos(2*x);
>> Hpoly(x,y,y1)

В результате получаем


Интерполяция функции sinx + sin2x кубическими полиномами Эрмита

Примечание про pchip. В состав MATLAB входит функция pchip, предназначенная для построения кусочных полиномов Эрмита. Однако в ней задаются только узлы xk и значения функции ƒ(xk). Значения производных ƒ'(xk) функция pchip назначает сама, вычисляя их описанным ниже способом.

Введем обозначение для разностной производной функции ƒ(x).

В каждом внутреннем узле xk для k = 2,3,...,n - 1 значения производной ƒ'(xk) определяют в зависимости от значений Δk и Δk-1 следующим образом. Возможны несколько случаев.

  1. Если знаки величин Δk и Δk-1 различны, или одна из этих величин равна нулю, то полагают

  2. Если знаки величин Δk и Δk-1 одинаковы и xk+1 - xk = xk - xk-1, то прибегают к усреднению

  3. Если знаки величин Δk и Δk-1 одинаковы и xk+1 - xk ≠ xk - xk-1, то выполняют такое усреднение

Таким образом, производная кусочного полинома в узлах интерполяции равна нулю или является результатом некоторого осреднения левой и правой разностных производных Δk-1 и Δk в k-ом узле интерполяции и наклон кусочного полинома в узле xk определяется наклонами двух отрезков, соединяющих точки (xk-1, ƒ(xk-1)), (xk, ƒ(xk)) и (xk, ƒ(xk)), (xk+1, ƒ(xk+1)) (см. рисунок ниже).

В граничной точке x1 для определения ƒ'(x1) используется трехточечная формула:

Значение ƒ' в правой граничной точке xn определяется по похожей формуле.

Вычисление коэффициентов кубического сплайна

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

Получение системы уравнений относительно тангенсов углов наклона и формул для остальных коэффициентов

Нам задана табличная функция

где xk - упорядоченные по возрастанию абсциссы узлов интерполирования. Задача состоит в нахождении кубического сплайна, обозначим его s(x) , который образован (n-1)-m кубическими

которых определен на отрезке [xk,xk+1] и которые должны удовлетворяют нижеприведенным условиям в узлах интерполирования, состоящим в прохождении графика сплайна через заданные точки и наличия у сплайна определенной гладкости.

1. Условия прохождения через заданные точки

2. Условия непрерывности первой производной во внутренних точках

3. Условия непрерывности второй производной во внутренних точках

Кроме этого, наложим на сплайн следующие граничные условия (для однозначного определения его коэффициентов).

1. На левой границе отрезка s(x1)=0, т.е

2. На правой границе отрезка s(xn)=0,

При выводе выражений для коэффициентов полиномов Эрмита (см. раздел Вычисление коэффициентов полиномов Эрмита) мы знали заданные значения производной во всех точках k=1,2,K, n и для коэффициентов bk (тангенсов углов наклона сплайна в узлах интерполирования) тогда получались простые выражения.

Теперь значения первых производных функции ƒ(x) в узлах интерполирования заранее неизвестны и мы поступим следующим образом: выберем коэффициенты bk таким образом, чтобы обеспечить непрерывность второй производной кубического сплайна во всех внутренних узлах на отрезке интерполирования.

Итак, для k=1,2,...,n-1:

ak=ƒ(xk)

bkследует определить

где

является разделенной разностью первого порядка интерполируемой функции ƒ(x), которая построена по точкам xk, xk+1.

Из условий непрерывности второй производной сплайна во внутренних узлах отрезка интерполирования

получаем, что

Подставим в полученные уравнения приведенные выше выражения для коэффициентов Ck и dk полиномов Pk(x) через неизвестные пока коэффициенты bk. Тогда получаем систему линейных алгебраических уравнений относительно коэффициентов bk, т.е. относительно тангенсов углов наклона кубического сплайна во всех (как внутренних, так и граничных) узлах отрезка интерполирования.

Здесь k=1,2,...,n-2.

Преобразуем систему уравнений, вычислив коэффициенты при bk, bk+1 и bk+2 , получим тогда

где введены следующие обозначения

Заметим, что поскольку мы разыскиваем сплайн с нулевыми наклонами в граничных точках, то b1= 0 и bn= 0 , но данные коэффициенты входят в первое и последнее уравнение системы. Поэтому система из (n-2)-ух уравнений относительно неизвестных коэффициентов bk для k=2,3,...,n-1 выглядит следующим образом:

Матрица системы линейных уравнений является трехдиагональной, все остальные ее коэффициенты, за исключением коэффициентов, стоящих на главной и побочных диагоналях, равны нулю. Решение системы с трехдиагональной матрицей осуществляется за число операций, пропорциональное числу неизвестных (например, методом прогонки). Решив эту систему мы находим коэффициенты b2,b3, .... bn-1 и подставляя их, а так же b1=bn= 0 в выражения для Ck и dk при к=1,2,...,n-1 находим все коэффициенты полиномов Pk(x) составляющих кубический сплайн s(x).

Для понимания приведенных выше выкладок может оказаться полезной функция mycubic10 (ее текст с комментариями приведен в разделе Функция mycubic10) Функция mycubic10 выполняет следующие действия.

1.По заданным упорядоченным по возрастанию абсциссам узлов интерполирования (в массивеx x) и значениям табличной функции в них (в массиве y) функция mycubic10 вычисляет значения кубического сплайна s(x) в заданных точках (в массиве xx), принадлежащих отрезку интерполирования [x1,xn]. Причем строится такой кубический сплайн, что значения его первой производной в первом и последнем узлах интерполирования равны нулю. Вычисленные значения возвращаются в массиве yy. . Причем строится такой кубический сплайн, что значения его первой производной в первом и последнем узлах интерполирования равны нулю. Вычисленные значения возвращаются в массиве yy.

2. Функция mycubic10 возвращает в двумерном массиве C значения коэффициентов каждого из полиномов Pk(x) составляющих кубический сплайн s(x) при этом каждая строка массива C содержит значения коэффициентов соответствующего полинома, начиная с коэффициента при старшей степени.

Функция mycubic10 содержит подфункцию CalcCoeffs, которая вычисляет значения коэффициентов ak, bk, ck и dk полиномов Pk(x), образующих кубический сплайн, выполняя следующие действия:

1) формирование трехдиагональной матрицы (в разреженном формате) и вектора правой части системы линейных алгебраических уравнений относительно тангенсов углов наклона сплайна b2, b3,.... bn-1 во внутренних узлах отрезка интерполирования;

2) решение этой системы линейных уравнений;

3) вычисление значений всех остальных коэффициентов a1, a2, .... , an-1, c1, c2, .... ,cn-1, d1, d2, .... , dn-1.

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

x = x(:)

Преобразование в вектор-столбцы происходит потому, что массивы в MATLAB хранятся по столбцам, а обращение x(:) служит для получения последовательно хранящихся элементов массива в виде одного столбца.

Найденные коэффициенты каждого полинома записываются в матрицу, являющуюся вторым выходным аргументом функции mycubic10. После этого происходит вычисление значений сплайна в заданных точках xx. Для вычисления значения сплайна в некоторой точке необходимо определить, какому из отрезков [xk,xk+1] данная точка для того, чтобы вычислять значение соответствующего полинома Pk(x) в этих точках и записать их на нужные позиции в вектор yy.

Определение индексов тех элементов массива xx, значения которых принадлежат текущему отрезку [xk,xk+1] , в функции mycubic10 осуществляется при помощи функция find. Входным аргументом функции find является логическое выражение, содержащее массив (& означает логическое «и»)

(xx >= x(k)) & (xx<= x(k+1))

а возвращает она массив индексов ind тех элементов массива xx, для которых данное логическое условие истинно.

Примечание

Если говорить более точно, то функция find работает следующим образом. Результатом логического выражения, в которое входит массив, является вектор из нулей и единиц, причем нули соответствуют тем позициям массива, для которых данное выражение ложно, а единицы — тем позициям массива, для которых выражение истинно, например

>> xx = [5.5 4.4 1.1 3.3 6.6 2.2];
>> z = xx<=4

z =
     0     0     1     1     0     1

Функция find ищет позиции всех ненулевых элементов массива и возвращает их.

>> ind = find(z)

ind =
     3     4     6

Далее сразу для всех элементов вектора xx с найденными при помощи функции find индексами вычисляются значения соответствующего полинома Pk(x) и результат записывается в элементы вектора yy с индексами, хранящимися в векторе ind. При этом применяются поэлементные операции.

Примечание

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

Остановимся более подробно на подфункции CalcCoeffs, принимающей на вход вектор-столбцы x и y с абсциссами и ординатами узлов интерполирования, соответственно, и возвращающей вектор-столбцы коэффициентов полиномов, составляющих кубический сплайн. Первый вектор-столбец d содержит коэффициенты полиномов при старшей степени, столбец c — при второй степени, столбец b — при первой степени, а столбец a — при нулевой.

После записи в переменную n количества узлов, в подфункции CalcCoeffs вычисляются расстояния между соседними узлами и они записываются в массив h, длина которого на единицу меньше длины массива x. Для вычисления расстояний между узлами используется функция diff. Таким образом

h(k) = x(k+1) –x(k)

Далее заполняются массивы alpha, beta и gamma, содержащие, соответственно, коэффициенты поддиагонали, основной диагонали и наддиагонали матрицы системы уравнений. При этом используются поэлементные операции и соответствующие части массива h. Так, например, заполнение массива beta по формуле

производится при помощи следующего оператора

beta = 2* (1./h(1:n-2) + 1./h(2:n-1));

Т.е. получается, что beta(k) = 2 * ( 1/h(k) + 1/h(k+1) ).

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

A = diag(beta) + diag(alpha, -1) + diag(gamma, 1);
Примечание.

Функция diag формирует диагональную матрицу из элементов вектора, указанного в ее входном аргументе. Функцию diag можно также вызывать со вторым аргументом, который задает номер заполняемой наддиагонали (1,2, …) или поддиагонали (-1, -2, …). Например

>> T = diag([3 3 3 3 3 3 3]) + diag([4 4 4 4 4 4], -1) + diag([2 2 2 2 2 2], 1)

T =
     3     2     0     0     0     0     0
     4     3     2     0     0     0     0
     0     4     3     2     0     0     0
     0     0     4     3     2     0     0
     0     0     0     4     3     2     0
     0     0     0     0     4     3     2
     0     0     0     0     0     4     3

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

Примечание

Для хранения разреженных матриц в MATLAB используется специальный тип данных double (sparse). Работе с разреженными матрицами посвящен раздел справочной системы MATLAB: Mathematics: Sparse Matrices.

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

D3 = [-1 -6  -11
          -2 -7  -12
          -3 -8  -13
          -4 -9  -14
          -5 -10 -15];

и формируем трехдиагональную матрицу 5 на 5 при помощи функции spdiags (аргумент [-1 0 1] содержит номер диагонали, на которую надо поставить соответствующий столбец матрицы D3, причем главная диагональ имеет номер 0, наддиагонали нумеруются положительными числами, поддиагонали — отрицательными)

  
A = spdiags(D3, [-1 0 1], 5, 5)

то получаем тогда представление трехдиагональной матрицы в следующем формате

A =
   (1,1)       -6
   (2,1)       -1
   (1,2)      -12
   (2,2)       -7
   (3,2)       -2
   (2,3)      -13
   (3,3)       -8
   (4,3)       -3
   (3,4)      -14
   (4,4)       -9
   (5,4)       -4
   (4,5)      -15
   (5,5)      -10

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

>> full(A)

ans =
    -6   -12     0     0     0
    -1    -7   -13     0     0
     0    -2    -8   -14     0
     0     0    -3    -9   -15
     0     0     0    -4   -10

После формирования трехдиагональной матрицы A заполняется вектор правой части delta системы линейных уравнений относительно тангенсов углов наклона сплайна в узлах интерполяции и находится решение системы при помощи оператора обратной косой черты. Если матрица системы уравнений трехдиагональная, то в MATLAB используется исключение по Гауссу, выполняемое только для ненулевых элементов.

Далее в функции CalcCoeffs вычисляются остальные коэффициенты полиномов, составляющих кубический сплайн.

Функция mycubic10

Данный раздел содержит основную функцию mycubic10 вместе с подфункцией CalcCoeffs, которые должны быть записаны в файл mycubic10.m

function [yy, C] = mycubic10(x, y, xx)
% Вычисление коэффициентов кубического сплайна и его значений в заданных точках

% Входные аргументы
%    x, y – узлы интерполяции
%    xx – значения аргумента, для которых надо вычислить кубический сплайн
%  Выходные аргументы  
%    yy – значения сплайна в xx
%    С – матрица коэффициентов кубических полиномов
%    С = [d(1)      c(1)       b(1)       a1
%            d(2)      c(2)       b(2)       a2 
%             ...
%            d(n-1)   c(n-1)   b(n-1)    a(n-1) ];

% находим число узлов интерполяции
n = length(x);

% преобразуем массивы узлов в вектор-столбцы
x = x(:);
y = y(:);

% находим коэффициенты полиномов, образующих сплайн
[d, c, b, a]  = CalcCoeffs(x, y);
% записываем найденные коэффициенты в матрицу
C = [d c b a];

% создаем нулевой вектор для записи значений сплайна
yy = zeros(size(xx));

% вычисляем значения сплайна
% в цикле проходим по отрезкам между узлами интерполирования
for k = 1:n-1
    % находим индексы элементов массива xx, значения которых 
    % лежат между x(k) и x(k+1)
    ind = find((xx >= x(k)) & (xx<= x(k+1)));
    % вычисляем значения соответствующего полинома
    yy(ind) = a(k) + b(k)*(xx(ind)-x(k)) + ...
        c(k)*(xx(ind)-x(k)).^2 + d(k)*(xx(ind)-x(k)).^3;
end
function [d, c, b, a]  = CalcCoeffs(x, y)
% подфункция для вычисления коэффициентов полиномов кубического сплайна

% находим число узлов
n = length(x);
% вычисляем вектор расстояний между узлами
h = diff(x);
% вычисляем коэффициенты трехдиагональной матрицы
alpha = 1./h(2:n-2);
beta = 2* (1./h(1:n-2) + 1./h(2:n-1));
gamma = 1./h(2:n-2);

% формируем прямоугольный массив с диагоналями, дополняя вектора 
% alpha и  gamma нулями
D3 = [[alpha; 0] beta [0; gamma]];
% формируем трехдиагональную матрицу в компактной форме
A = spdiags(D3, [-1 0 1], n-2, n-2);

% вычисляем правую часть системы
df = diff(y)./h;
delta = 3*(df(2:n-1)./h(2:n-1) + df(1:n-2)./h(1:n-2));

% решаем систему относительно тангенсов углов наклона сплайна в узлах
b = A\delta;

% расширяем полученный вектор, добавляя граничные условия
b = [0; b; 0];
% находим остальные коэффициенты полиномов
a = y(1:n-1);
c = (3*df - b(2:n) - 2*b(1:n-1))./h;
d = (b(1:n-1) + b(2:n) - 2*df)./h.^2;
% удаляем последний элемент вектора b, который содержит тангенс угла наклона в x(n),
% т.к. он не является коэффициентом последнего полинома
b(end) = [];

Проверим работу функции mycubic10, проинтерполировав кубическим сплайном (с нулевыми условиями на первую производную на границах отрезка интерполирования) следующую табличную функцию

>> x = [-4 -3.5 -2  0 2.1 2.5 5];
>> y = [0.5    0.2   -0.7  0    1.1    0.8 -1.6];

Для этого вызываем функцию mycubic10

>> [yy, C] = mycubic10(x, y, xx);

и строим по полученным значениями yy график кубического сплайна, на котором маркерами отложены табличные данные

>> plot(x, y, 'o', xx, yy)
>> xlim([x(1)-0.1 x(end)+0.1])

Получаем следующий график

Приближение табличной функции кубическим сплайном.

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

Проведем более подробную проверку, построив графики первой, второй и третьей производной кубического сплайна. Для этого вычислим производные полученного сплайна в 30 равноотстоящих точках на каждом из отрезков [xk,xk+1], на которых полиномы имеют следующий вид

Тогда

В цикле по отрезкам [xk,xk+1] вычисляем значения производных кубического сплайна и строим их графики различными цветами на каждом из отрезков, предварительно создав три пары осей в графическом окне при помощи функции subplot.

Примечание

Функция subplot(m, n, k) создает в графическом окне пару осей в соответствии с выбранным способом разбиения окна (m осей по вертикали и n по горизонтали) и указанным номером пары осей k (оси нумеруются слева направо и сверху вниз).

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

% создание графического окна белого цвета
figure('Color', 'w')

% создание трех пар осей с сеткой
subplot(3,1,1)
hold on
grid on
subplot(3,1,2)
hold on
grid on
subplot(3,1,3)
hold on
grid on

% проход по отрезкам между узлами интерполирования и построение графиков
for k = 1 : length(x) - 1 
    % извлечение коэффициентов полинома на текущем отрезке из матрицы C
    d = C(k,1);
    c = C(k,2);
    b = C(k,3);
    a = C(k,4);
    % создание 30-и равноотстоящих точек на отрезке [x(k), x(k+1)]
    xx = linspace(x(k), x(k+1), 30);
    % вычисление производных k-го полинома
    d1p = b + 2*c*(xx - x(k)) + 3*d*(xx - x(k)).^2;
    d2p = 2*c + 6*d*(xx - x(k));    
    d3p = 6*d*ones(size(xx));
    % задание случайного цвета
    randcolor = rand(1,3);    
    % построение графика первой производной 
    subplot(3,1,1)
    plot(xx, d1p, 'Color', randcolor, 'LineWidth', 2)
    % построение графика второй производной
    subplot(3,1,2)
    plot(xx, d2p, 'Color', randcolor, 'LineWidth', 2)
    % построение графика третьей производной
    subplot(3,1,3)
    plot(xx, d3p, 'Color', randcolor, 'LineWidth', 2)     
end

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

Первая, вторая и третья производные кубического сплайна.

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

>> pp = csape(x,[0 y 0],[1 1]);
>> Cpp = pp.coefs;

Получаем совпадение коэффициентов (с точностью до ошибок округления)

>> err = max(max(abs(pp.coefs - Cpp)))
err =
  4.4409e-016

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

Система Orphus

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