Численные методы

1. Интерполяция
2. Сплайн-интерполяции
3. Решение систем линейных уравнений
4. Линейное программирование
5. Решение нелинейных уравнений
6. Решение систем нелинейных уравнений
7. Минимум функции одной переменной
8. Минимум функции нескольких переменных
9. Интегрирование

Интерполяция

В практике часто возникает необходимость аппроксимации некоторой табличной функции, заданной массивом аргументов \(x_1, x_2, n_3 \ldots, x_n\) и массивом значений функции в этих точках \(y(x_1), y(x_2), y(x_3), \ldots, y(x_n)\). Под аппроксимацией понимается подбор некоторой просто вычисляемой функции, которая приближала бы значения заданной табличной функции.

Известны значения некоторой функции \(y(x)\) в точках \(x_0, x_1, x_2, \ldots, x_n\): \[y_0, y_1, y_2, \ldots, y_n.\]

Необходимо построить интерполирующий многочлен \(L_{n}(x)\), совпадающий со значениями табличной функции в точках \(x_0, x_1, x_2, \ldots, x_n\) и приближающий функцию \(y(x)\) на интервале \([x_0,x_n]\).

Для определения коэффициентов интерполяционного полинома, проходящего через заданные точки с координатами \(x_1, f(x_1)\), \(x_2, f(x_2)\), … \(x_n, f(x_n)\) используется функция polyfit:

Первый аргумент функции массив значений \(x\), второй – массив значений функции.

>> x = [ 1  2  3  4 ];
>> y = [ 0  0.6931  1.0986  1.3863 ];
>> p = polyfit(x, y, 3)

p =
    0.0283   -0.3137    1.4362   -1.1507

Результат работы функции polyfit – коэффициенты полинома, расположенные в порядке убывания степени аргумента, при которых эти множители стоят: \[P_3(x) = 0.0283 x^3 - 0.3137 x^2 + 1.4362 x - 1.1507\]

Для вычисления значения полинома используется функция polyval. Первый аргумент функции polyval – массив коэффициентов полинома, полученный при помощи функции polyfit, второй – одно или несколько значений аргумента \(x\), для которых необходимо вычислить значение полиномиальной функции:

>> x0 = 0.1;
>> polyval(p,x0)

ans =
   -1.0100

% График табличной функции (точки)
plot(x, y, 'ro', 'LineWidth', 2);
xlabel('x');
ylabel('y');
hold on;
% График полинома
fplot(@(xa) polyval(p,xa),[x(1), x(end)],'LineWidth',2);
hold off;

Если степень многочлена, которая передается в функцию polyfit меньше, чем уменьшенное на 1 количество точек табличной функции, то коэффициенты полинома вычисляются по методу наименьших квадратов.

% Столбец значений аргумента: 50 точек равномерно расположенных в интервале от 0 до 5  
x  = linspace(0,5,50)';

Квадрат аргумента с “шумом”: к каждому вычисленному значению добавляется случайное число от 0 до 5

y  = x.^2 + rand(50,1)*5;

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

p1 = polyfit(x, y, 1);
p2 = polyfit(x, y, 2);

Точками строим график табличной функции

plot(x,y,'bo')

Удерживая (не стирая) первый график (hold on) рисуем на этом же рисунке графики аппроксимирующих многочленов при помощи функции fplot, первым аргументом которой является анонимная функция от x, вызывающая функцию polyval:

hold on;
fplot(@(x) polyval(p2,x), [0 5], 'r', 'LineWidth', 2)
fplot(@(x) polyval(p1,x), [0 5], 'g', 'LineWidth', 2)
hold off;
xlabel('x'); ylabel('y');

Сплайн-интерполяции

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

Для сплайн интерполяции используется функция spline. Первый аргумент функции – строка или столбец значений аргумента, второй – строка или столбец со значениями функции, соответствующими значениям x. Третий аргумент – значения x, для которых необходимо вычислить значения функции, используя кубическую интерполяцию

x  = [0 1 2 3];
y  = exp(x);

xq = linspace(0,3,10);
yq = spline(x,y,xq)

plot(x,y,'ro');
hold on;
plot(xq,yq,'b-');
hold off;
legend('Табличная функция','Сплайн');

Для линейно интерполяции используется функция interp1 с такой же структурой аргументов, что и функция spline. Отличие в том, что известные точки соединяются прямой линией, т.е. для определения промежуточных значений используется линейная интерполяция

x  = [0 1 2 3];
y  = exp(x);

xq = linspace(0,3,10);
yq = interp1(x,y,xq)

plot(x,y,'ro');
hold on;
plot(xq,yq,'b-');
hold off;
legend('Табличная функция','Линейная интерполяция');

Решение систем линейных уравнений

Для решения системы линейных уравнений, например \[\left\{ \begin{aligned} 2 x_1 + 3 x_2 + x_3 = 5\\ x_1 + 7 x_2 + 2 x_3 = 1\\ 3 x_1 + 2 x_2 + x_3 = 2 \end{aligned} \right.\]

необходимо представить её в матричном виде \[\boldsymbol{A} \cdot \boldsymbol x = \boldsymbol B\]

где \[\boldsymbol A = \begin{bmatrix} 2 & 3 & 1 \\ 1 & 7 & 2 \\ 3 & 2 & 1\end{bmatrix}, \quad \boldsymbol B = \begin{bmatrix} 5 \\ 1 \\ 2 \end{bmatrix}, \quad \boldsymbol x = \begin{bmatrix} x_1 \\ x_2 \\ x_3 \end{bmatrix}, \quad\]

Для определения столбца неизвестных \(\boldsymbol x\) используется оператор \

>> A = [2 3 1; 1 7 2; 3 2 1];
>> B = [5;1;2];
>> x = A\B

x =
    6.0000
    9.0000
  -34.0000

Линейное программирование

Линейное программирование — это математический численный метод для оптимизации моделей, в которых целевые функции \[f(x_1,x_2,\ldots,x_n) \rightarrow min\]

и ограничения, например \[g(x_1,x_2,\ldots,x_n) > 0\]

являются уравнениями линейной алгебры.

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

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

Прибыль (целевая функция) \[S = x_1 s_1 + x_2 s_2 + x_3 s_3 + x_4 s_4 \rightarrow max\]

Ограничения \[x_1 a_{11} + x_2 a_{12} + x_3 a_{13} + x_4 a_{14} \leq 1200\] \[x_1 a_{21} + x_2 a_{22} + x_3 a_{23} + x_4 a_{24} \leq 600\] \[x_1 a_{31} + x_2 a_{32} + x_3 a_{33} + x_4 a_{34} \leq 1500\]

В MATLAB для решения этой задачи используется функция linprog

x = linprog(f,Ane,Bne,Ae,Be,lb,ub)

При условиях равенствах \[\boldsymbol{A}_e \boldsymbol{x} = \boldsymbol{B_e}\]

При условиях неравенствах \[\boldsymbol{A}_{ne} \boldsymbol{x} \leq \boldsymbol{B_{ne}}\]

При условии \[\boldsymbol{L} \leq \boldsymbol{x} \leq \boldsymbol{U}\]

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

>> f = -[15 6 12 24];

Условия-неравенства (запасы сырья на складе) описываются матрицей коэффициентов Ane и матрицей запасов Bne

>> Ane = [4 2 1 8; 2 10 6 0; 3 0 6 1];
>> Bne = [1200; 600; 1500];

Вызываем функцию linprog

>> x = linprog(f,Ane,Bne,[],[],[0 0 0 0],[])

Optimal solution found.

x =

         0
         0
  100.0000
  137.5000

Таким образом, оптимальным планом будет производство 100 ед продукта 3 и 137 ед продукта 4.

Решение нелинейных уравнений

Для решения уравнений вида \[f(x) = 0\]

используется функция fzero, первым аргументом которой является ссылка на функцию \(f(x)\), вторым – начальное приближение для искомого значения решения уравнения или интервал внутри которого находится решение:

>> fzero( @(x) cos(x) - x, 1)

ans =
    0.7391

Решение систем нелинейных уравнений

Для решения системы нелинейных уравнений используется функция fsolve. Например, для решения системы \[\left\{ \begin{aligned} x^2 + y^2 + z^2 - 1 & = 0 \\ 2x^2 + y^2 - 4z & = 0 \\ 3x^2 - 4y + z^2 & = 0 \end{aligned} \right.\]

необходимо написать m-файл с функцией

% Файл f.m
function y = f(x)
    y = [x(1)^2+x(2)^2+x(3)^2-1;
         2*x(1)^2+x(2)^2-4*x(3);
         3*x(1)^2-4*x(2)+x(3)^2];

Для поиска решения необходимо вызвать функцию fsolve, передав ей ссылку на функцию и вектор начального приближения \([x_0, y_0, z_0]\)

>> fsolve(@y,[1;1;1])

ans =
    0.7852
    0.4966
    0.3699

Минимум функции одной переменной

Поиск минимума функции \(f(x)\) на интервале \([x_1, x_2]\)

>> fun = @(x) x^2 - 0.3*x^3 + 3 + cos(4*x)^2;
>> fplot(fun,[-1 3]);
>> fminbnd(fun,-1,3)

ans = 
    0.372

Минимум функции нескольких переменных

Найти минимум функции нескольких переменных \[\min_{\boldsymbol x} f(\boldsymbol x), \quad \boldsymbol{x} = [x_1,x_2,\ldots,x_n]\]

Пример: \[f(x)=100(x_2 − x_1^2)^2+(1−x_1)^2\]

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

>> fun = @(x1,x2) 100*(x2-x1.^2).^2+(1-x1).^2;
>> x1 = linspace(0,3,100);
>> x2 = linspace(0,3,100);
>> [X1,X2] = meshgrid(x1,x2);
>> f = fun(X1,X2);
>> contour(X1,X2,f,200);
>> fminsearch(@(x) fun(x(1),x(2)),[0, 0.5])

ans =
    1.0000    1.0000

Интегрирование

Для функции, заданной таблично, интеграл можно вычислить используя формулу трапеций, используя функцию trapz:

>> x = 1.0:0.1:2.0;
>> y = log(x);
>> trapz(x, y)

ans =
    0.3859

Точное значение интеграла равно \(\log 4 - 1 \approx 0.386294361119891\).

Для известной подинтегральной функции более точный результат дает использование функции quad, в которую нужно передать ссылку на подинтегральную функцию и пределы интегрирования \[F = \int_1^2 \ln(x) dx\]

>> format long
>> quad(@log, 1, 2)

ans =
   0.386294334336416

Третьим аргументом функции quad является желаемая точность результата (по умолчанию – \(10^{-6}\) )

>> f = @(x) sin(x).*log(x).^2;
>> quad(f, 1, 2, 1e-8)

ans =
   0.1822