Дифференцирование в Octave осуществляется в технике символьных переменных. Для работы с символьными переменными в Octave подключите специальный пакет расширений octave-symbolic. Установка пакетов расширений описана в первой главе, техника работы с символьными переменными — в п. 2.7. В функциях интегрирования реализованы различные численные алгоритмы.
Дифференцирование в Octave осуществляется с помощью функции $$differentiate(f(a,x [,n]))$$, где $$a$$ — символьное выражение, $$x$$ — переменная дифференцирования, $$n$$ — порядок дифференцирования (при $$n = 1$$ параметр можно опустить). Иными словами, функция вычисляет $$n$$-ю производную выражения $$a$$ по переменной $$x$$.
Производной функции $$f (x)$$ в точке $$x_0$$ называется предел, к которому стремится отношение бесконечно малого приращения функции к соответствующему бесконечно малому приращению аргумента. Геометрический смысл этого понятия заключается в том, что если к графику функции $$f (x)$$ провести касательную в точке $$x_0$$, то её угловой коэффициент, будет равен значению производной в этой точке $$k = f’(x)$$. Следовательно, уравнение касательной к линии в заданной точке имеет вид: $$y(x) = f’(x)(x-x-0)+f(x-0)$$
(рис 8.1) График функции и её касательной
Пример 8.1. Записать уравнение касательной к функции$$f(x)=\frac{3x^2-7}{2x+1}$$, в точке $$x-o=-1$$.
Из листинга 8.1 видим, что уравнение касательной к функции в заданной точке имеет вид y$$y(x)=14x+18$$.
clear all;
x0=-1;
symbols
x = sym("x");
f =(3 - x^2-7) /(2 - x+1);
f1=differentiate(f,x)% Первая производная от заданной функции
% Уравнение касательной: k=subs(f1,x,x0), f(x0)=subs(f,x,x0)
y=subs (f1, x, x0) - (x-x0)+subs(f, x, x0)
f1 =(6.0)*x*(1.0+(2.0)*x)^(-1)-(2.0)*(-7.0+(3.0)*x^(2.0))*(1.0+(2.0)*x)^(-2)
y =18.0-9.029803704631804845E-19*I+(14.0-6.0198691364212032297E-19*I)*x
На рис. 8.1 представлены графики заданной функции и её касательной. Рисунок построен с помощью команд из листинга 8.2.
clear all; clf; cla;
symbols
x=sym ("x");
L1=ezplot(’(3*x^2-7)/(2*x+1)’);
set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’18+14*x’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’, ’Marker’, ’o’)
set(gca, ’xlim’, [-4, 2]); set(gca, ’ylim’, [-6, 10]);
set(gca, ’xtick’, [-4:2]); set(gca, ’ytick’, [-6:2:10]);
grid on; xlabel(’x’); ylabel(’y’);
title(’(3*x^2-7)/(2*x+1), 18+14*x’);
Пример 8.2. Найти a) $$f'(x)=\frac{5sin(2x)}{\sqrt{cos(2x)}$$ и б) $$f'(x)=tan\sqrt[3]{ln(x)}$$.
Решение примера показано влистинге 8.3.
>>> clear all;
>>> symbols
>>> x = sym("x");
% Пример а)
>>> f =(5- Sin(2 - x))/Sqrt(Cos(2 - x));
>>> f1=differentiate(f, x)
f1 =(5.0)*sin((2.0)*x)^2*cos((2.0)*x)^(-3/2)+(10.0)*sqrt(cos((2.0)*x))
% Пример б)
>>> f=Tan(Log(x)^(1/3));
>>> f1=differentiate(f, x)
f1 = (0.333)*(1+tan(log(x)^(0.333))^2)*x^(-1)*log(x)^(-0.666)
Если функция $$y(x)$$ задана параметрическими уравнениями $$x=\phi(t),y=\psi(t)$$, то производная вычисляется по формуле $$x=\frac{\phi'(t)}{\psi'(t)}=\frac{y'_t}{x_t}$$
Пример 8.3. Найти производную функции, заданной параметрически:
$$\left\{\begin{aligned}x(t)=3cos^3(t)\\y(t)=3sin^3(t)\end{aligned}$$Листинг 8.4 содержит решение примера.
>>> clear all;
>>> symbols
>>> t = sym("t");
>>> x=3-Cos(t)^3;
>>> y=3-Sin(t)^3;
>>> xt=differentiate(x, t);
>>> yt=differentiate(y, t);
>>> f=yt/xt
f = -sin(t) - cos(t) ^( -1.0)
Пример 8.4. Найти производные a) $$f''(x)=ln(cos(x))$$ б) $$f^{IV}(x)=tan(x$$ (Листинг 8.5).
>>> clear all;
>>> symbols
>>> x=sym("x");
% Пример а)
>>> f=Log(Cos(x));
>>> differentiate(f, x, 2)
ans = -1-cos(x)^(-2)-sin(x)^2
% Пример б)
>>> f=Tan(x);
>>> differentiate(f, x, 4)
ans = 16-(1+tan(x)^2)^2-tan(x)+8-(1+tan(x)^2)-tan(x)^3
Пример 8.5. Найти производные $$y''(x)$$
$$\left\{\begin{aligned}x(t)=t-sin(t)\\y(t)=t-cos(t)\end{aligned}$$Выражение для вычисления второй производной параметрической функции: $$y''(x)=\frac{\left(\frac{y_t'}{x_t'}\right)'}{x_t'}=\frac{y_t''x_t'-x_t''y_t'}{(x_t')^3}$$.
В листинге 8.6 представлено решение примера
>>> clear all;
>>> symbols
>>> t=sym("t");
>>> x=t-Sin(t);
>>> y=1-Cos(t);
>>> xt=differentiate(x, t);
>>> yt=differentiate(y, t);
>>> xt2=differentiate(x, t, 2);
>>> yt2=differentiate(y, t, 2);
>>> z=(yt2-xt-xt2-yt)/xt^3
z=-(1-cos(t))^(-3.0)-(cos(t)-(-1+cos(t))+sin(t)^2)
(рис 8.2) Исследование функции на возрастание и убывание
Понятие производной тесно связано с задачей исследования функции. Из курса математического анализа известно, что если производная функции $$f (x)$$ положительна на всём интервале $$[a, b]$$, то функция на нём возрастает, если всюду отрицательна, то $$f (x)$$ убывает.
Пример 8.6. Построить график функции $$f(x)=1-2x-x^2$$ и её производной. Исследовать функцию на возрастание и убывание.
Вычислим производную заданной функции и построим оба графика в одном окне :
>>> symbols
>>> x=sym("x");
>>> f=1-2-x-x^2;
>>> differentiate(f, x)
ans=-2.0-(2.0)-x
% построим график заданной функции и её производной.
clf; cla;
L1=ezplot(’1-2*x-x^2’); set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’-2-2*x’); set(L2, ’LineWidth’, 2, ’Color’, ’k’)
grid on; xlabel(’x’); ylabel(’y’); title(’1-2x-x^2, -2-2x’)
На рис. 8.2 видим, что там, где $$y=f'(x)$$ принимает положительные значения, $$f (x)$$ возрастает, соответственно, при отрицательныхзначениях $$y=f'(x)$$ функция $$f (x)$$ убывает.
Говорят, что непрерывная функция $$f(x)$$ имеет максимум в точке $$x = a$$, если в достаточной близости от этой точки производная $$f'(x)$$ положительна слева от a и отрицательна справа от a. Если наоборот, то $$f(x)$$ имеет минимум в точке $$x = a$$. Максимум и минимум объединяют названием экстремум. Если первая производная в этой точке $$f'(a)$$ либо равна нулю, либо не существует, то в этой точке может быть экстремум.
Пример 8.7. Исследовать функцию $$f(x)=\frac{x^3}{3}-2x^2+3x+1$$ на экстремум.
Найдём производную функции и отобразим её на графике:
clear all;
symbols
x=sym("x");
f=x^3/3-2-x^2+3-x+1;
% Производная от функции f(x)
y=differentiate(f, x)
y=3.0+x^(2.0)-(4.0)-x
% Изобразим функцию и её производную на графике
clf; cla;
L1=ezplot(’x^3/3-2*x^2+3*x+1’);
set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’3.0+x^(2.0)-(4.0)*x’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’)
set(gca, ’xlim’, [-2,5]); set(gca, ’ylim’, [-5,5]);
grid on; xlabel(’x’); ylabel(’y’);
title(’x^3/3-2x^2+3x+1, 3+x^2-4x’)
% Корни уравнения
>>> x1 = symfsolve(y, 1)
>>> x2 = symfsolve(y, 3)
x1 = 1
x2 = 3
На рис. 8.3 и в листинге 8.8 видно, что первая производная обращается в нуль в точках $$x = 1$$ и $$x = 3$$. При переходе через точку $$x = 1\ f'(x)$$ меняет знак с плюса на минус, следовательно, это точка максимума функции $$f (x)$$, а в точке $$x = 3$$ знак первой производной меняется с минуса на плюс, то есть это точка минимума.
График функции называется выпуклым на промежутке $$[a, b]$$, если он расположен выше касательной, проведённой в любой точке этого интервала. Если же график функции лежит ниже касательной, то он называется вогнутым. Функция будет выпуклой на интервале $$[a, b]$$, если вторая производная $$f''(x)$$ на нём положительна. И наоборот, если вторая производная отрицательна, то функция вогнута. Если же вторая производная равна нулю в некоторой точке $$a$$, а слева и справа от неё имеет значения разных знаков, то точка $$a$$ — точка перегиба.
(рис 8.3) Исследование функции на экстремум
Пример 8.8. Определить точки перегиба функции $$f(x)=\frac{3x-2}{x^2+1}$$
Найдём вторую производную заданной функции. Построим графики функции и её второй производной. Определим точки в которых вторая производная обращается в ноль (листинг 8.9).
clear all;
symbols
x=sym("x");
f=(3-x-2)/(x^2+1);
y=differentiate(f, x, 2)
y=-(2.0)*(1.0+x^(2.0))^(-2)*(-2.0+(3.0)*x)-(12.0)*(1.0+x^(2.0))^(-2)
*x+(8.0)*(1.0+x^(2.0))^(-3)*x^2*(-2.0+(3.0)*x)
clf; cla;
L1=ezplot(’(3*x-2)/(x^2+1)’); set(L1, ’LineWidth’, 4, ’Color’, ’k’)
hold on
L2=ezplot(’-(2.0)*(1.0+x^(2.0))^(-2)*(-2.0+(3.0)*x)-(12.0)*(1.0+x^(2.0))^(-2)
*x+(8.0)*(1.0+x^(2.0))^(-3)*x^2*(-2.0+(3.0)*x)’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’)
set(gca, ’xlim’, [-5, 5]); set(gca, ’ylim’, [-5, 7]);
grid on; xlabel(’x’); ylabel(’y’); title(’ ’)
>>> x1 = symfsolve(y, -1)
>>> x2 = symfsolve(y, 0)
>>> x3 = symfsolve(y,2)
x1 = -1.1411
x2 = 0.19855
x3 = 2.9425
(рис 8.4) Определение точки перегиба функции
Иллюстрации приведены на рис. 8.4. Исследование второй производной функции $$f''(x)$$ показывает, что она определена на всей числовой оси и обращается в нуль в трёх точках $$x_1=-1.1411,x_2=0.19855,x_3=2.9425$$, причём при переходе через них она меняет знак. Следовательно, на интервале $$(-\infty,x_1)$$ функция $$f (x)$$ вогнутая, так как $$f'' (x) < 0$$, на $$(x_1,x_2)$$ — выпуклая $$( f''(x) > 0)$$, на $$(x_2,x_3)$$ — вогнутая $$( f''(x) < 0)$$ и на $$(x_3,+\infty)$$ опять выпуклая, потому что $$f''(x) > 0$$.
Пусть дана функция $$f (x)$$, известно, что она непрерывна на интервале $$[a, b]$$ и уже определена её первообразная $$F (x)$$, тогда определённый интеграл от этой функции можно вычислить в пределах от $$a$$ до $$b$$ по формуле Ньютона–Лейбница:
$$\int\limits_a^b f(x)dx=F(b)-F(a),$ где $F'(x)=f(x)$$Пример 8.9. Вычислить определённый интеграл
$$I=\int\limits_2^5\sqrt{2x-1}dx$$К сожалению в Octave не предусмотрены средства символьного интегрирования, поэтому обратимся к таблице интегралов и найдём, что
$$I=\int\sqrt{2x-1}dx=\frac{1}{3}\sqrt[3]{(2x-1)^2}+C$$Теперь вычислим интеграл по формуле Ньютона–Лейбница:
clear all; % Функция, определяющая подынтегральное выражение % x — переменная интегрирования, C — постоянная интегрирования. function y=F(x,C) y=1/3-(2-x-1)^(3/2)+C; end; >>> a =2; b=5; % Вычисление интеграла по формуле Ньютона-Лейбница >>> I = F(b, 0)-F(a, 0) I = 7.2679
На практике часто встречаются интегралы с первообразной, которая не может быть выражена через элементарные функции или является слишком сложной, что затрудняет, или делает невозможным, вычисления по формуле Ньютон-Лейбница. Кроме того, нередко подынтегральная функция задаётся таблицей или графиком и тогда понятие первообразной вообще теряет смысл. В этом случае большое значение имеют численные методы интегрирования, основная задача которых заключается в вычислении значения определённого интеграла на основании значений подынтегральной функции.
Численное вычисление определённого интеграла называют механической квадратурой. Формулы, соответствующие тому или иному численному методу приближённого интегрирования, называют квадратурными. Подобное название связано с геометрическим смыслом определённого интеграла: значение определённого интеграла
$$I=\int\limits_a^b f(x)dx,$f(x)\not=0$$равно площади криволинейной трапеции с основаниями $$[a, b]$$ и $$f (x)$$.
Вообще говоря, классические учебники по численной математике предлагают немало методов интегрирования, но здесь мы рассмотрим только те методы, которые имеют непосредственное отношение к функциям Octave.
(рис 8.5) Геометрическая интерпретация метода трапеций
Изложим геометрическую интерпретацию интегрирования по методу трапеций. Для этого участок интегрирования $$[a, b]$$ разобьём точками на n равных частей (рис. 8.5), причём $$x_0=a, x_n=b$$
Тогда длина каждой части будет равна $$h=\frac{b-a}{n}$$, а значение абсциссы каждой из точек разбиения можно вычислить по формуле $$x_i=x_0+ih,i=1,2,...,n-1$$. Теперь из каждой точки $$x_i$$ проведём перпендикуляр до пересечения с кривой $$f (x)$$, а затем заменим каждую из полученных криволинейных трапеций прямолинейной. Приближённое значение интеграла будем рассматривать как сумму площадей прямолинейных трапеций, причём площадь отдельной трапеции составляет $$S_i=\frac{y_{i-1}+y_i}{2}h$$, следовательно, площадь искомой фигуры вычисляют по формуле:
$$S=\int\limits_a^b f(x)dx=\sum \limits_{i=1}^nS_i=\frac{h}{2}\sum \limits_{i=1}^n(y_{i-1}+y_i)=h\left(\frac{y_0+y_n}{2}+\sum \limits_{i=1}^ny_i\right)$$Таким образом, получена квадратурная формула трапеций для численного интегрирования:
$$I=\int\limits_a^b f(x)dx=h\left(\frac{f(a)+f(b)}{2}+\sum \limits_{i=1}^{n-1}f(x_i)\right).$$Функции $$trapz$$ и $$cumtrapz$$ реализуют численное интегрирование по методу трапеций в Octave.
Площадь фигуры под графиком функции $$y(x)$$, в котором все точки заданы векторами $$x$$ и $$y$$, вычисляет команда $$trapz (x, y)$$.
| x | -1.5708 | -1.0708 | -0.5708 | -0.0708 | 0.4292 | 0.9292 | 1,4292 |
| y | 0 | 0.47943 | 0.84147 | 0.99749 | 0.90930 | 0.59847 | 0.14112 |
Если вызвать функцию $$trapz (y)$$ с одним аргументом, то будет вычислена площадь фигуры под графиком функции $$y(x)$$, в котором все точки заданы векторами $$x$$ и $$y$$, причём по умолчанию элементы вектора $$x$$ принимают значения номеров элементов вектора $$y$$.
Пример 8.10. Вычислить интеграл от функции y(x) = cos(x). Значения функции представлены в табл. 8.1.
Решение примера представлено в листинге 8.11.
>>> clear all; >>> x=[-1.5708 -1.0708 -0.5708 -0.0708 0.4292 0.9292 1.4292]; >>> y=[0 0.47943 0.84147 0.99749 0.90930 0.59847 0.14112]; >>> I=trapz(x, y) I = 1.9484
Пример 8.11. Вычислить интеграл $$I=\int\limits_2^5\sqrt{2x-1}dx$$
Листинг 8.12 содержит несколько вариантов решения данного примера. В первом случае интервал интегрирования делится на отрезки с шагом 1, во втором 0.5, в третьем 0.1 и в четвёртом 0.05. Не трудно заметить, что чем больше точек разбиения, тем точнее значение искомого интеграла. Решение можно сравнить с результатом полученным в задаче 8.9, где этот же интеграл был найден по формулам Ньютона-Лейбница (листинг 8.10).
clear all; % Вариант 1. h=1 >>> x = 2:5; y=sqrt(2-x-1); I1=trapz(x, y) % Вариант 2. h=0.5 >>> x = 2:0.5:5; y=sqrt(2-x-1); I2=trapz(x, y) % Вариант 3. h=0.1 >>> x = 2:0.1:5; y=sqrt(2-x-1); I3=trapz(x, y) % Вариант 4. h=0.05 >>> x = 2:0.05:5; y=sqrt(2-x-1); I4=trapz(x, y) % Результаты интегрирования I1 = 7.2478 I2 = 7.2629 I3 = 7.2677 I4 = 7.2679
В листинге 8.13 приведён пример использования функции $$trapz$$ с одним аргументом. Как видим, в первом случае значение интеграла, вычисленного при помощи этой функции, не точно и совпадает со значением, полученным функцией $$trapz (x, y)$$ на интервале [2, 5] с шагом 1 (листинг 8.12, первый вариант). То есть мы нашли сумму площадей трёх прямолинейных трапеций с основанием $$h = 1$$ и боковыми сторонами, заданными вектором y. Во втором случае, при попытке увеличить точность интегрирования, значение интеграла существенно увеличивается. Дело в том что, уменьшив шаг разбиения интервала интегрирования до 0.05, мы увеличили количество элементов векторов $$x$$ и $$y$$ и применение функции $$trapz (y)$$ приведёт к вычислению суммы площадей шестидесяти трапеций с основанием $$h = 1$$ и боковыми сторонами, заданными вектором $$y$$. Таким образом, в первом и втором примерах листинга 8.13 вычисляются площади совершенно разных фигур.
% Пример 1. >>> x = 2:5; y=sqrt(2-x-1); I=trapz(y) I = 7.2478 % Пример 2. >>> x = 2:0.05:5; y=sqrt(2-x-1); I=trapz(y) I = 145.36
Функция $$cumtrapz$$ выполняет так называемое "интегрирование с накоплением" по методу трапеций. Это означает, что она, так же как и $$trapz$$, вычисляет площадь фигуры под графиком функции $$y(x)$$, но результатом её работы является вектор, состоящий из промежуточных вычислений. То есть, если общая площадь $$S$$ криволинейной трапеции сформирована из суммы площадей $$\sum_{i=1}^nS_i$$ прямолинейных трапеций, то элементы вектора представляют собой следующую последовательность $$S_1=0,S_2=S_1+S_2,S_3=S_1+S_2+S_3,...,S_n=S_1+S_2+S_3+\cdot\cdot\cdot+S_n$$.
Таким образом, последний элемент вектора будет равен искомой площади фигуры $$S$$. Функцию интегрирования с накоплением можно вызывать в форматах $$cumtrapz(x, y)$$ и $$cumtrapz(y)$$, где $$x$$ и $$y$$ векторы, определяющие функцию $$y(x)$$.
Пример 8.12. Вычислить интеграл $$I=\int\limits_{0}^{\frac{\pi}{2}}\frac{1}{5+sin(x)}dx$$
Листинг 8.14 демонстрирует применение функции интегрирования с накоплением $$cumtrapz$$ к поставленной задаче. Там же приведена интерпретация работы этой функции с помощью команды $$trapz$$.
>>> x = 0:0.1:pi/2; y=(5+sin(x)).^(-1); % 1. Интегрирование с накоплением >>> I1=cumtrapz(x, y) I1 = Columns 1 through 8: 0.0 0.0198 0.03923 0.05829 0.07701 0.09541 0.11352 0.13136 Columns 9 through 16: 0.1489 0.1663 0.18356 0.2006 0.2175 0.23434 0.25108 0.2677 % 2. Обычное интегрирование >>> I2=trapz(x, y) % Значение I2 совпадает с последним значением вектора I1 I2 = 0.26777 % 3. Интегрирование на левой части интервала от 0 до 1 >>> x = 0:0.1:1; y=(5+sin(x)).^(-1); >>> I3=trapz(x, y) % Значение I3 совпадает с 11-м значением вектора I1 I3 = 0.18356
Изложим идею интегрирования по методу Симпсона. Пусть $$n = 2m$$ — чётное число, а $$y_i=f(x_i)(i=0,1,...,n)$$ — значения функции $$y = f (x)$$ для равноотстоящих точек $$a=x_0,x_1,x_2,...,x_n =b$$ с шагом $$h=\frac{b-a}{n}=\frac{b-a}{2m}$$. На паре соседних участков (рис. 8.6) кривая $$y = f (x)$$ заменяется параболой $$y = L(x)$$, коэффициенты которой подобраны так, что она проходит через точки $$Y_0,Y_1,Y_2$$
Площадь криволинейной трапеции, ограниченной сверху параболой, составит: $$S_i=\frac{h}{3}(y_{i-1}+4y_i+y_{i+1})$$
Суммируя площади всех криволинейных трапеций, получим:
(рис 8.6) Геометрическая интерпретация интегрирования по методу Симпсона
$$S=\int\limits_{a}^bf(x)dx=\frac{h}{3}(y_0+4y_1+2y_2+4y_3+\cdot\cdot\cdot+2y_{2m-2}+4y_{2m-1}+y_{2m})=\frac{h}{3}\left(y_0+y_{2m}+\sum \limits_{i=1}^{2m-1}py_i\right)$$
,
Следовательно, формула Симпсона для численного интегрирования имеет вид:
$$I=\int\limits_{a}^bf(x)dx=\frac{h}{3}\left(f(a)+f(b)+\sum \limits_{i=1}^{2m-1}py_i\right)$$Методы трапеций и Симпсона являются частными случаями квадратурных формул Ньютона–Котеса, которые, вообще говоря, имеют вид,
$$\int\limits_{a}^bydx=(b-a)\sum \limits_{i=0}^nH_iy_i$$где $$H_i$$ — это некоторые константы называемые постоянными Ньютона–Котеса.
Если для квадратурных формул Ньютона-Котеса принять $$n = 1$$, то получим метод трапеций, а при $$n = 2$$ — метод Симпсона. Поэтому эти методы называют квадратурными методами низших порядков. Для $$n > 2$$ получают квадратурные формулы Ньютона–Котеса высших порядков.
В Octave реализован вычислительный алгоритм метода Симпсона с автоматическим выбором шага. Автоматический выбор шага интегрирования заключается в том, что интервал интегрирования разбивают на n отрезков и вычисляют значение интеграла, если полученное значение не удовлетворяет заданной точности вычислений, то n увеличивают вдвое и вновь вычисляют значение интеграла, так повторяют до тех пор пока не будет достигнута заданная точность. Итак, вычисление интеграла по методу Симпсона обеспечивает функция $$[F, K] = quadv(name, a, b[, tol, trace])$$, где:
$$name$$ — имя функции, задающей подынтегральное выражение;
$$a, b$$ — пределы интегрирования;
$$tol$$ — точность вычислений;
$$trace$$ — параметр позволяющий получить информацию о ходе вычислений в виде таблицы, в столбцах которой представлены: значение количества вычислений, начальная точка текущего промежутка интегрирования, его длина и значение интеграла;
$$F$$ — значение интеграла;
$$K$$ — количество итераций.
Пример 8.13. Вычислить интеграл $$\int\limits_{0}^1\sqrt{4-x^2}dx$$
Решение примера с применением функции $$quadv$$ приведено влистинге 8.15.
% Подынтегральная функция function y=G(x) y=(4-x^2).^(1/2); end; >>> format long % Вычисление интеграла по методу Симпсона % Точность установлена по умолчанию 1.0e-06 >>> [F1, K1]=quadv(’G’, 0, 1) % Результат — значение интеграла и количество итераций F1 = 1.91322288999134 K1 = 17 % Точность установлена пользователем 1.0e-07 >>> [F2, K2]=quadv(’G’, 0, 1, 1.0 e-07) % Результат — значение интеграла и количество итераций F2 = 1.91322295090669 K2 = 33 >>> format short % Вызов функций с заданной степенью точности % Вывод дополнительной информации о вычислениях >>> quadv (’G’, 0, 1, 1.0 e-05,5) 5.00000 0.00000 1.00000 1.91321 7.00000 0.00000 0.50000 0.98948 9.00000 0.50000 0.50000 0.92374 11.0000 0.50000 0.25000 0.47458 13.0000 0.75000 0.25000 0.44916 ans = 1.9132
Запишем в общем виде квадратурную формулу для функции заданной на промежутке [-1; 1]
$$\int\limits_{-1}^1ydx=\sum \limits_{i=0}^nA_if(x_i)$$Попытаемся найти коэффициенты $$A_i$$ и узловые точки $$t_i$$, таким образом, чтобы квадратурная формула была точной для всех полиномов $$f(t)=1,t,t^2,...,t^{2n-1}$$
В этом случае построение квадратурной формулы приводит к определению $$A_i$$ и $$t_i$$ из нелинейной системы $$2n$$ уравнений:
$$\left\{ \begin{array}{ll} \displaystyle\sum _{i=1}^{n}A_{i}=2,\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}=0,\\ \hdotsfor{2}\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}^{2n-2}=\frac{2}{2n-1},\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}^{2n-1}=0. \end{array} \right.$$Решение нелинейной системы задача не тривиальная, но её можно обойти, если знать, что значениями $$t_i$$ квадратурной формулы служат корни многочлена Лежандра
$$P_n(t)=\frac{1}{2^n\cdot n!}\cdot\frac{d^n}{dy^n}((t^2-1)^2),(n=0,1,2,...)$$Как известно, корни полинома Лежандра существуют при любом $$n$$, различны и принадлежат интервалу [-1; 1].
Итак, квадратурной формулой Гаусса называют выражение
$$\int\limits_{a}^bf(x)dx=\frac{b-a}{2}\sum \limits_{i=1}^nA_if \left(\frac{a+b}{2}+\frac{b-a}{2}t_i\right)$$где $$t_i$$ — корни полинома Лежандра, а $$A_i$$ определяется интегрирова-нием базисных многочленов Лежандра $$P_i(t)$$ степени $$n - 1$$:
$$A_i=\int\limits_{-1}^1\frac{(t-t_1)...(t-t_{i-1})(t-t_{i+1})...(t-tn)}{(t-t_1)...(t-t_{i-1})(t-t_{i+1})...(t_i-tn)}dt$$В Octave интегрирование по квадратуре Гаусса выполняет функция $$[F, kod, K, err] = quad(name, a, b, tol, sing)$$, где: $$name$$ — имя функции, задающей подынтегральное выражение; $$a, b$$ — пределы интегрирования; $$tol$$ — точность вычислений; $$sing$$ — вектор значений, близких к тем, в которых подынтегральная функция терпит разрыв; $$F$$ — значение интеграла; $$kod$$ — код ошибки в решении (0 — решение завершено успешно); $$K$$ — количество итераций; $$err$$ — погрешность вычислений.
Пример 8.14. Вычислить интеграл $$\int\limits_{0}^1t^2\sqrt{\left(3+sin \left(\frac{1}{t}\right)\right)}dt$$
Обратите внимание, что в нижней границе интегрирования подынтегральная функция терпит разрыв. Решение примера с применением функции $$quad$$ приведено в листинге 8.16.
clear all; function y=f(x) y=(x.^2).-sqrt(3+sin(1./x)); end; >>> format long >>> [F, kod, K, err]=quad(’f’, 0, 1) F = 0.654343719149802 kod = 0 K = 1323 err = 1.37557012147481 e-08 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-05) F = 0.654343738854992 kod = 0 K = 315 err = 3.82733563379833 e -06 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-20) F = 0.654343718970708 kod = 0 K = 1491 err = 9.39557628735834 e-09 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-20,0.1) F = 0.654343710193938 kod = 0 K = 840 err = 5.80259740257105 e-09 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e- 20,0.001) F = 0.654343718720156 kod = 0 K = 1596 err = 8.35248716562893 e-09
Функции
$$F = quadl(f, a, b[, tol, trace])$$ и $$[F, err] = quadgk(f, a, b[, tol, trace])$$, где: $$name$$ — имя функции, задающей подынтегральное выражение; $$a, b$$ — пределы интегрирования; $$tol$$ — точность вычислений; $$trace$$ — таблица промежуточных вычислений; $$F$$ — значение интеграла; $$err$$ — погрешность вычислений; также выполняют интегрирование по квадратуре Гаусса. В этих функциях специальным образом подбирается шаг. В первом случае по методу Гаусса-Лобатто, во втором Гаусса-Конрада.
Пример 8.15. Вычислить интеграл $$\int\limits_{-\frac{\pi}{3}}^{\frac{\pi}{3}}tg^4(x)dx$$
Решение примера с применением функций $$quadl$$ и $$quadgk$$ приведено в листинге 8.17.
function y=f(x) y = tan(x).^4; end; >>> format long >>> [F]= quadl (’f’,-pi/3, pi/3, 1.0 e-05) F = 2.09439512983937 >>> [F, err]=quadgk(’f’, -pi/3, pi/3, 1.0 e-05) F = 2.09439510239319 err = 1.02555919485880 e-12
Дифференцирование в Octave осуществляется в технике символьных переменных. Для работы с символьными переменными в Octave подключите специальный пакет расширений octave-symbolic. Установка пакетов расширений описана в первой главе, техника работы с символьными переменными — в п. 2.7. В функциях интегрирования реализованы различные численные алгоритмы.
Дифференцирование в Octave осуществляется с помощью функции $$differentiate(f(a,x [,n]))$$, где $$a$$ — символьное выражение, $$x$$ — переменная дифференцирования, $$n$$ — порядок дифференцирования (при $$n = 1$$ параметр можно опустить). Иными словами, функция вычисляет $$n$$-ю производную выражения $$a$$ по переменной $$x$$.
Производной функции $$f (x)$$ в точке $$x_0$$ называется предел, к которому стремится отношение бесконечно малого приращения функции к соответствующему бесконечно малому приращению аргумента. Геометрический смысл этого понятия заключается в том, что если к графику функции $$f (x)$$ провести касательную в точке $$x_0$$, то её угловой коэффициент, будет равен значению производной в этой точке $$k = f’(x)$$. Следовательно, уравнение касательной к линии в заданной точке имеет вид: $$y(x) = f’(x)(x-x-0)+f(x-0)$$
(рис 8.1) График функции и её касательной
Пример 8.1. Записать уравнение касательной к функции$$f(x)=\frac{3x^2-7}{2x+1}$$, в точке $$x-o=-1$$.
Из листинга 8.1 видим, что уравнение касательной к функции в заданной точке имеет вид y$$y(x)=14x+18$$.
clear all;
x0=-1;
symbols
x = sym("x");
f =(3 - x^2-7) /(2 - x+1);
f1=differentiate(f,x)% Первая производная от заданной функции
% Уравнение касательной: k=subs(f1,x,x0), f(x0)=subs(f,x,x0)
y=subs (f1, x, x0) - (x-x0)+subs(f, x, x0)
f1 =(6.0)*x*(1.0+(2.0)*x)^(-1)-(2.0)*(-7.0+(3.0)*x^(2.0))*(1.0+(2.0)*x)^(-2)
y =18.0-9.029803704631804845E-19*I+(14.0-6.0198691364212032297E-19*I)*x
На рис. 8.1 представлены графики заданной функции и её касательной. Рисунок построен с помощью команд из листинга 8.2.
clear all; clf; cla;
symbols
x=sym ("x");
L1=ezplot(’(3*x^2-7)/(2*x+1)’);
set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’18+14*x’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’, ’Marker’, ’o’)
set(gca, ’xlim’, [-4, 2]); set(gca, ’ylim’, [-6, 10]);
set(gca, ’xtick’, [-4:2]); set(gca, ’ytick’, [-6:2:10]);
grid on; xlabel(’x’); ylabel(’y’);
title(’(3*x^2-7)/(2*x+1), 18+14*x’);
Пример 8.2. Найти a) $$f'(x)=\frac{5sin(2x)}{\sqrt{cos(2x)}$$ и б) $$f'(x)=tan\sqrt[3]{ln(x)}$$.
Решение примера показано влистинге 8.3.
>>> clear all;
>>> symbols
>>> x = sym("x");
% Пример а)
>>> f =(5- Sin(2 - x))/Sqrt(Cos(2 - x));
>>> f1=differentiate(f, x)
f1 =(5.0)*sin((2.0)*x)^2*cos((2.0)*x)^(-3/2)+(10.0)*sqrt(cos((2.0)*x))
% Пример б)
>>> f=Tan(Log(x)^(1/3));
>>> f1=differentiate(f, x)
f1 = (0.333)*(1+tan(log(x)^(0.333))^2)*x^(-1)*log(x)^(-0.666)
Если функция $$y(x)$$ задана параметрическими уравнениями $$x=\phi(t),y=\psi(t)$$, то производная вычисляется по формуле $$x=\frac{\phi'(t)}{\psi'(t)}=\frac{y'_t}{x_t}$$
Пример 8.3. Найти производную функции, заданной параметрически:
$$\left\{\begin{aligned}x(t)=3cos^3(t)\\y(t)=3sin^3(t)\end{aligned}$$Листинг 8.4 содержит решение примера.
>>> clear all;
>>> symbols
>>> t = sym("t");
>>> x=3-Cos(t)^3;
>>> y=3-Sin(t)^3;
>>> xt=differentiate(x, t);
>>> yt=differentiate(y, t);
>>> f=yt/xt
f = -sin(t) - cos(t) ^( -1.0)
Пример 8.4. Найти производные a) $$f''(x)=ln(cos(x))$$ б) $$f^{IV}(x)=tan(x$$ (Листинг 8.5).
>>> clear all;
>>> symbols
>>> x=sym("x");
% Пример а)
>>> f=Log(Cos(x));
>>> differentiate(f, x, 2)
ans = -1-cos(x)^(-2)-sin(x)^2
% Пример б)
>>> f=Tan(x);
>>> differentiate(f, x, 4)
ans = 16-(1+tan(x)^2)^2-tan(x)+8-(1+tan(x)^2)-tan(x)^3
Пример 8.5. Найти производные $$y''(x)$$
$$\left\{\begin{aligned}x(t)=t-sin(t)\\y(t)=t-cos(t)\end{aligned}$$Выражение для вычисления второй производной параметрической функции: $$y''(x)=\frac{\left(\frac{y_t'}{x_t'}\right)'}{x_t'}=\frac{y_t''x_t'-x_t''y_t'}{(x_t')^3}$$.
В листинге 8.6 представлено решение примера
>>> clear all;
>>> symbols
>>> t=sym("t");
>>> x=t-Sin(t);
>>> y=1-Cos(t);
>>> xt=differentiate(x, t);
>>> yt=differentiate(y, t);
>>> xt2=differentiate(x, t, 2);
>>> yt2=differentiate(y, t, 2);
>>> z=(yt2-xt-xt2-yt)/xt^3
z=-(1-cos(t))^(-3.0)-(cos(t)-(-1+cos(t))+sin(t)^2)
(рис 8.2) Исследование функции на возрастание и убывание
Понятие производной тесно связано с задачей исследования функции. Из курса математического анализа известно, что если производная функции $$f (x)$$ положительна на всём интервале $$[a, b]$$, то функция на нём возрастает, если всюду отрицательна, то $$f (x)$$ убывает.
Пример 8.6. Построить график функции $$f(x)=1-2x-x^2$$ и её производной. Исследовать функцию на возрастание и убывание.
Вычислим производную заданной функции и построим оба графика в одном окне :
>>> symbols
>>> x=sym("x");
>>> f=1-2-x-x^2;
>>> differentiate(f, x)
ans=-2.0-(2.0)-x
% построим график заданной функции и её производной.
clf; cla;
L1=ezplot(’1-2*x-x^2’); set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’-2-2*x’); set(L2, ’LineWidth’, 2, ’Color’, ’k’)
grid on; xlabel(’x’); ylabel(’y’); title(’1-2x-x^2, -2-2x’)
На рис. 8.2 видим, что там, где $$y=f'(x)$$ принимает положительные значения, $$f (x)$$ возрастает, соответственно, при отрицательныхзначениях $$y=f'(x)$$ функция $$f (x)$$ убывает.
Говорят, что непрерывная функция $$f(x)$$ имеет максимум в точке $$x = a$$, если в достаточной близости от этой точки производная $$f'(x)$$ положительна слева от a и отрицательна справа от a. Если наоборот, то $$f(x)$$ имеет минимум в точке $$x = a$$. Максимум и минимум объединяют названием экстремум. Если первая производная в этой точке $$f'(a)$$ либо равна нулю, либо не существует, то в этой точке может быть экстремум.
Пример 8.7. Исследовать функцию $$f(x)=\frac{x^3}{3}-2x^2+3x+1$$ на экстремум.
Найдём производную функции и отобразим её на графике:
clear all;
symbols
x=sym("x");
f=x^3/3-2-x^2+3-x+1;
% Производная от функции f(x)
y=differentiate(f, x)
y=3.0+x^(2.0)-(4.0)-x
% Изобразим функцию и её производную на графике
clf; cla;
L1=ezplot(’x^3/3-2*x^2+3*x+1’);
set(L1, ’LineWidth’, 3, ’Color’, ’k’)
hold on
L2=ezplot(’3.0+x^(2.0)-(4.0)*x’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’)
set(gca, ’xlim’, [-2,5]); set(gca, ’ylim’, [-5,5]);
grid on; xlabel(’x’); ylabel(’y’);
title(’x^3/3-2x^2+3x+1, 3+x^2-4x’)
% Корни уравнения
>>> x1 = symfsolve(y, 1)
>>> x2 = symfsolve(y, 3)
x1 = 1
x2 = 3
На рис. 8.3 и в листинге 8.8 видно, что первая производная обращается в нуль в точках $$x = 1$$ и $$x = 3$$. При переходе через точку $$x = 1\ f'(x)$$ меняет знак с плюса на минус, следовательно, это точка максимума функции $$f (x)$$, а в точке $$x = 3$$ знак первой производной меняется с минуса на плюс, то есть это точка минимума.
График функции называется выпуклым на промежутке $$[a, b]$$, если он расположен выше касательной, проведённой в любой точке этого интервала. Если же график функции лежит ниже касательной, то он называется вогнутым. Функция будет выпуклой на интервале $$[a, b]$$, если вторая производная $$f''(x)$$ на нём положительна. И наоборот, если вторая производная отрицательна, то функция вогнута. Если же вторая производная равна нулю в некоторой точке $$a$$, а слева и справа от неё имеет значения разных знаков, то точка $$a$$ — точка перегиба.
(рис 8.3) Исследование функции на экстремум
Пример 8.8. Определить точки перегиба функции $$f(x)=\frac{3x-2}{x^2+1}$$
Найдём вторую производную заданной функции. Построим графики функции и её второй производной. Определим точки в которых вторая производная обращается в ноль (листинг 8.9).
clear all;
symbols
x=sym("x");
f=(3-x-2)/(x^2+1);
y=differentiate(f, x, 2)
y=-(2.0)*(1.0+x^(2.0))^(-2)*(-2.0+(3.0)*x)-(12.0)*(1.0+x^(2.0))^(-2)
*x+(8.0)*(1.0+x^(2.0))^(-3)*x^2*(-2.0+(3.0)*x)
clf; cla;
L1=ezplot(’(3*x-2)/(x^2+1)’); set(L1, ’LineWidth’, 4, ’Color’, ’k’)
hold on
L2=ezplot(’-(2.0)*(1.0+x^(2.0))^(-2)*(-2.0+(3.0)*x)-(12.0)*(1.0+x^(2.0))^(-2)
*x+(8.0)*(1.0+x^(2.0))^(-3)*x^2*(-2.0+(3.0)*x)’);
set(L2, ’LineWidth’, 2, ’Color’, ’k’)
set(gca, ’xlim’, [-5, 5]); set(gca, ’ylim’, [-5, 7]);
grid on; xlabel(’x’); ylabel(’y’); title(’ ’)
>>> x1 = symfsolve(y, -1)
>>> x2 = symfsolve(y, 0)
>>> x3 = symfsolve(y,2)
x1 = -1.1411
x2 = 0.19855
x3 = 2.9425
(рис 8.4) Определение точки перегиба функции
Иллюстрации приведены на рис. 8.4. Исследование второй производной функции $$f''(x)$$ показывает, что она определена на всей числовой оси и обращается в нуль в трёх точках $$x_1=-1.1411,x_2=0.19855,x_3=2.9425$$, причём при переходе через них она меняет знак. Следовательно, на интервале $$(-\infty,x_1)$$ функция $$f (x)$$ вогнутая, так как $$f'' (x) < 0$$, на $$(x_1,x_2)$$ — выпуклая $$( f''(x) > 0)$$, на $$(x_2,x_3)$$ — вогнутая $$( f''(x) < 0)$$ и на $$(x_3,+\infty)$$ опять выпуклая, потому что $$f''(x) > 0$$.
Пусть дана функция $$f (x)$$, известно, что она непрерывна на интервале $$[a, b]$$ и уже определена её первообразная $$F (x)$$, тогда определённый интеграл от этой функции можно вычислить в пределах от $$a$$ до $$b$$ по формуле Ньютона–Лейбница:
$$\int\limits_a^b f(x)dx=F(b)-F(a),$ где $F'(x)=f(x)$$Пример 8.9. Вычислить определённый интеграл
$$I=\int\limits_2^5\sqrt{2x-1}dx$$К сожалению в Octave не предусмотрены средства символьного интегрирования, поэтому обратимся к таблице интегралов и найдём, что
$$I=\int\sqrt{2x-1}dx=\frac{1}{3}\sqrt[3]{(2x-1)^2}+C$$Теперь вычислим интеграл по формуле Ньютона–Лейбница:
clear all; % Функция, определяющая подынтегральное выражение % x — переменная интегрирования, C — постоянная интегрирования. function y=F(x,C) y=1/3-(2-x-1)^(3/2)+C; end; >>> a =2; b=5; % Вычисление интеграла по формуле Ньютона-Лейбница >>> I = F(b, 0)-F(a, 0) I = 7.2679
На практике часто встречаются интегралы с первообразной, которая не может быть выражена через элементарные функции или является слишком сложной, что затрудняет, или делает невозможным, вычисления по формуле Ньютон-Лейбница. Кроме того, нередко подынтегральная функция задаётся таблицей или графиком и тогда понятие первообразной вообще теряет смысл. В этом случае большое значение имеют численные методы интегрирования, основная задача которых заключается в вычислении значения определённого интеграла на основании значений подынтегральной функции.
Численное вычисление определённого интеграла называют механической квадратурой. Формулы, соответствующие тому или иному численному методу приближённого интегрирования, называют квадратурными. Подобное название связано с геометрическим смыслом определённого интеграла: значение определённого интеграла
$$I=\int\limits_a^b f(x)dx,$f(x)\not=0$$равно площади криволинейной трапеции с основаниями $$[a, b]$$ и $$f (x)$$.
Вообще говоря, классические учебники по численной математике предлагают немало методов интегрирования, но здесь мы рассмотрим только те методы, которые имеют непосредственное отношение к функциям Octave.
(рис 8.5) Геометрическая интерпретация метода трапеций
Изложим геометрическую интерпретацию интегрирования по методу трапеций. Для этого участок интегрирования $$[a, b]$$ разобьём точками на n равных частей (рис. 8.5), причём $$x_0=a, x_n=b$$
Тогда длина каждой части будет равна $$h=\frac{b-a}{n}$$, а значение абсциссы каждой из точек разбиения можно вычислить по формуле $$x_i=x_0+ih,i=1,2,...,n-1$$. Теперь из каждой точки $$x_i$$ проведём перпендикуляр до пересечения с кривой $$f (x)$$, а затем заменим каждую из полученных криволинейных трапеций прямолинейной. Приближённое значение интеграла будем рассматривать как сумму площадей прямолинейных трапеций, причём площадь отдельной трапеции составляет $$S_i=\frac{y_{i-1}+y_i}{2}h$$, следовательно, площадь искомой фигуры вычисляют по формуле:
$$S=\int\limits_a^b f(x)dx=\sum \limits_{i=1}^nS_i=\frac{h}{2}\sum \limits_{i=1}^n(y_{i-1}+y_i)=h\left(\frac{y_0+y_n}{2}+\sum \limits_{i=1}^ny_i\right)$$Таким образом, получена квадратурная формула трапеций для численного интегрирования:
$$I=\int\limits_a^b f(x)dx=h\left(\frac{f(a)+f(b)}{2}+\sum \limits_{i=1}^{n-1}f(x_i)\right).$$Функции $$trapz$$ и $$cumtrapz$$ реализуют численное интегрирование по методу трапеций в Octave.
Площадь фигуры под графиком функции $$y(x)$$, в котором все точки заданы векторами $$x$$ и $$y$$, вычисляет команда $$trapz (x, y)$$.
| x | -1.5708 | -1.0708 | -0.5708 | -0.0708 | 0.4292 | 0.9292 | 1,4292 |
| y | 0 | 0.47943 | 0.84147 | 0.99749 | 0.90930 | 0.59847 | 0.14112 |
Если вызвать функцию $$trapz (y)$$ с одним аргументом, то будет вычислена площадь фигуры под графиком функции $$y(x)$$, в котором все точки заданы векторами $$x$$ и $$y$$, причём по умолчанию элементы вектора $$x$$ принимают значения номеров элементов вектора $$y$$.
Пример 8.10. Вычислить интеграл от функции y(x) = cos(x). Значения функции представлены в табл. 8.1.
Решение примера представлено в листинге 8.11.
>>> clear all; >>> x=[-1.5708 -1.0708 -0.5708 -0.0708 0.4292 0.9292 1.4292]; >>> y=[0 0.47943 0.84147 0.99749 0.90930 0.59847 0.14112]; >>> I=trapz(x, y) I = 1.9484
Пример 8.11. Вычислить интеграл $$I=\int\limits_2^5\sqrt{2x-1}dx$$
Листинг 8.12 содержит несколько вариантов решения данного примера. В первом случае интервал интегрирования делится на отрезки с шагом 1, во втором 0.5, в третьем 0.1 и в четвёртом 0.05. Не трудно заметить, что чем больше точек разбиения, тем точнее значение искомого интеграла. Решение можно сравнить с результатом полученным в задаче 8.9, где этот же интеграл был найден по формулам Ньютона-Лейбница (листинг 8.10).
clear all; % Вариант 1. h=1 >>> x = 2:5; y=sqrt(2-x-1); I1=trapz(x, y) % Вариант 2. h=0.5 >>> x = 2:0.5:5; y=sqrt(2-x-1); I2=trapz(x, y) % Вариант 3. h=0.1 >>> x = 2:0.1:5; y=sqrt(2-x-1); I3=trapz(x, y) % Вариант 4. h=0.05 >>> x = 2:0.05:5; y=sqrt(2-x-1); I4=trapz(x, y) % Результаты интегрирования I1 = 7.2478 I2 = 7.2629 I3 = 7.2677 I4 = 7.2679
В листинге 8.13 приведён пример использования функции $$trapz$$ с одним аргументом. Как видим, в первом случае значение интеграла, вычисленного при помощи этой функции, не точно и совпадает со значением, полученным функцией $$trapz (x, y)$$ на интервале [2, 5] с шагом 1 (листинг 8.12, первый вариант). То есть мы нашли сумму площадей трёх прямолинейных трапеций с основанием $$h = 1$$ и боковыми сторонами, заданными вектором y. Во втором случае, при попытке увеличить точность интегрирования, значение интеграла существенно увеличивается. Дело в том что, уменьшив шаг разбиения интервала интегрирования до 0.05, мы увеличили количество элементов векторов $$x$$ и $$y$$ и применение функции $$trapz (y)$$ приведёт к вычислению суммы площадей шестидесяти трапеций с основанием $$h = 1$$ и боковыми сторонами, заданными вектором $$y$$. Таким образом, в первом и втором примерах листинга 8.13 вычисляются площади совершенно разных фигур.
% Пример 1. >>> x = 2:5; y=sqrt(2-x-1); I=trapz(y) I = 7.2478 % Пример 2. >>> x = 2:0.05:5; y=sqrt(2-x-1); I=trapz(y) I = 145.36
Функция $$cumtrapz$$ выполняет так называемое "интегрирование с накоплением" по методу трапеций. Это означает, что она, так же как и $$trapz$$, вычисляет площадь фигуры под графиком функции $$y(x)$$, но результатом её работы является вектор, состоящий из промежуточных вычислений. То есть, если общая площадь $$S$$ криволинейной трапеции сформирована из суммы площадей $$\sum_{i=1}^nS_i$$ прямолинейных трапеций, то элементы вектора представляют собой следующую последовательность $$S_1=0,S_2=S_1+S_2,S_3=S_1+S_2+S_3,...,S_n=S_1+S_2+S_3+\cdot\cdot\cdot+S_n$$.
Таким образом, последний элемент вектора будет равен искомой площади фигуры $$S$$. Функцию интегрирования с накоплением можно вызывать в форматах $$cumtrapz(x, y)$$ и $$cumtrapz(y)$$, где $$x$$ и $$y$$ векторы, определяющие функцию $$y(x)$$.
Пример 8.12. Вычислить интеграл $$I=\int\limits_{0}^{\frac{\pi}{2}}\frac{1}{5+sin(x)}dx$$
Листинг 8.14 демонстрирует применение функции интегрирования с накоплением $$cumtrapz$$ к поставленной задаче. Там же приведена интерпретация работы этой функции с помощью команды $$trapz$$.
>>> x = 0:0.1:pi/2; y=(5+sin(x)).^(-1); % 1. Интегрирование с накоплением >>> I1=cumtrapz(x, y) I1 = Columns 1 through 8: 0.0 0.0198 0.03923 0.05829 0.07701 0.09541 0.11352 0.13136 Columns 9 through 16: 0.1489 0.1663 0.18356 0.2006 0.2175 0.23434 0.25108 0.2677 % 2. Обычное интегрирование >>> I2=trapz(x, y) % Значение I2 совпадает с последним значением вектора I1 I2 = 0.26777 % 3. Интегрирование на левой части интервала от 0 до 1 >>> x = 0:0.1:1; y=(5+sin(x)).^(-1); >>> I3=trapz(x, y) % Значение I3 совпадает с 11-м значением вектора I1 I3 = 0.18356
Изложим идею интегрирования по методу Симпсона. Пусть $$n = 2m$$ — чётное число, а $$y_i=f(x_i)(i=0,1,...,n)$$ — значения функции $$y = f (x)$$ для равноотстоящих точек $$a=x_0,x_1,x_2,...,x_n =b$$ с шагом $$h=\frac{b-a}{n}=\frac{b-a}{2m}$$. На паре соседних участков (рис. 8.6) кривая $$y = f (x)$$ заменяется параболой $$y = L(x)$$, коэффициенты которой подобраны так, что она проходит через точки $$Y_0,Y_1,Y_2$$
Площадь криволинейной трапеции, ограниченной сверху параболой, составит: $$S_i=\frac{h}{3}(y_{i-1}+4y_i+y_{i+1})$$
Суммируя площади всех криволинейных трапеций, получим:
(рис 8.6) Геометрическая интерпретация интегрирования по методу Симпсона
$$S=\int\limits_{a}^bf(x)dx=\frac{h}{3}(y_0+4y_1+2y_2+4y_3+\cdot\cdot\cdot+2y_{2m-2}+4y_{2m-1}+y_{2m})=\frac{h}{3}\left(y_0+y_{2m}+\sum \limits_{i=1}^{2m-1}py_i\right)$$
,
Следовательно, формула Симпсона для численного интегрирования имеет вид:
$$I=\int\limits_{a}^bf(x)dx=\frac{h}{3}\left(f(a)+f(b)+\sum \limits_{i=1}^{2m-1}py_i\right)$$Методы трапеций и Симпсона являются частными случаями квадратурных формул Ньютона–Котеса, которые, вообще говоря, имеют вид,
$$\int\limits_{a}^bydx=(b-a)\sum \limits_{i=0}^nH_iy_i$$где $$H_i$$ — это некоторые константы называемые постоянными Ньютона–Котеса.
Если для квадратурных формул Ньютона-Котеса принять $$n = 1$$, то получим метод трапеций, а при $$n = 2$$ — метод Симпсона. Поэтому эти методы называют квадратурными методами низших порядков. Для $$n > 2$$ получают квадратурные формулы Ньютона–Котеса высших порядков.
В Octave реализован вычислительный алгоритм метода Симпсона с автоматическим выбором шага. Автоматический выбор шага интегрирования заключается в том, что интервал интегрирования разбивают на n отрезков и вычисляют значение интеграла, если полученное значение не удовлетворяет заданной точности вычислений, то n увеличивают вдвое и вновь вычисляют значение интеграла, так повторяют до тех пор пока не будет достигнута заданная точность. Итак, вычисление интеграла по методу Симпсона обеспечивает функция $$[F, K] = quadv(name, a, b[, tol, trace])$$, где:
$$name$$ — имя функции, задающей подынтегральное выражение;
$$a, b$$ — пределы интегрирования;
$$tol$$ — точность вычислений;
$$trace$$ — параметр позволяющий получить информацию о ходе вычислений в виде таблицы, в столбцах которой представлены: значение количества вычислений, начальная точка текущего промежутка интегрирования, его длина и значение интеграла;
$$F$$ — значение интеграла;
$$K$$ — количество итераций.
Пример 8.13. Вычислить интеграл $$\int\limits_{0}^1\sqrt{4-x^2}dx$$
Решение примера с применением функции $$quadv$$ приведено влистинге 8.15.
% Подынтегральная функция function y=G(x) y=(4-x^2).^(1/2); end; >>> format long % Вычисление интеграла по методу Симпсона % Точность установлена по умолчанию 1.0e-06 >>> [F1, K1]=quadv(’G’, 0, 1) % Результат — значение интеграла и количество итераций F1 = 1.91322288999134 K1 = 17 % Точность установлена пользователем 1.0e-07 >>> [F2, K2]=quadv(’G’, 0, 1, 1.0 e-07) % Результат — значение интеграла и количество итераций F2 = 1.91322295090669 K2 = 33 >>> format short % Вызов функций с заданной степенью точности % Вывод дополнительной информации о вычислениях >>> quadv (’G’, 0, 1, 1.0 e-05,5) 5.00000 0.00000 1.00000 1.91321 7.00000 0.00000 0.50000 0.98948 9.00000 0.50000 0.50000 0.92374 11.0000 0.50000 0.25000 0.47458 13.0000 0.75000 0.25000 0.44916 ans = 1.9132
Запишем в общем виде квадратурную формулу для функции заданной на промежутке [-1; 1]
$$\int\limits_{-1}^1ydx=\sum \limits_{i=0}^nA_if(x_i)$$Попытаемся найти коэффициенты $$A_i$$ и узловые точки $$t_i$$, таким образом, чтобы квадратурная формула была точной для всех полиномов $$f(t)=1,t,t^2,...,t^{2n-1}$$
В этом случае построение квадратурной формулы приводит к определению $$A_i$$ и $$t_i$$ из нелинейной системы $$2n$$ уравнений:
$$\left\{ \begin{array}{ll} \displaystyle\sum _{i=1}^{n}A_{i}=2,\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}=0,\\ \hdotsfor{2}\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}^{2n-2}=\frac{2}{2n-1},\\ \displaystyle\sum _{i=1}^{n}A_{i}\cdot t_{i}^{2n-1}=0. \end{array} \right.$$Решение нелинейной системы задача не тривиальная, но её можно обойти, если знать, что значениями $$t_i$$ квадратурной формулы служат корни многочлена Лежандра
$$P_n(t)=\frac{1}{2^n\cdot n!}\cdot\frac{d^n}{dy^n}((t^2-1)^2),(n=0,1,2,...)$$Как известно, корни полинома Лежандра существуют при любом $$n$$, различны и принадлежат интервалу [-1; 1].
Итак, квадратурной формулой Гаусса называют выражение
$$\int\limits_{a}^bf(x)dx=\frac{b-a}{2}\sum \limits_{i=1}^nA_if \left(\frac{a+b}{2}+\frac{b-a}{2}t_i\right)$$где $$t_i$$ — корни полинома Лежандра, а $$A_i$$ определяется интегрирова-нием базисных многочленов Лежандра $$P_i(t)$$ степени $$n - 1$$:
$$A_i=\int\limits_{-1}^1\frac{(t-t_1)...(t-t_{i-1})(t-t_{i+1})...(t-tn)}{(t-t_1)...(t-t_{i-1})(t-t_{i+1})...(t_i-tn)}dt$$В Octave интегрирование по квадратуре Гаусса выполняет функция $$[F, kod, K, err] = quad(name, a, b, tol, sing)$$, где: $$name$$ — имя функции, задающей подынтегральное выражение; $$a, b$$ — пределы интегрирования; $$tol$$ — точность вычислений; $$sing$$ — вектор значений, близких к тем, в которых подынтегральная функция терпит разрыв; $$F$$ — значение интеграла; $$kod$$ — код ошибки в решении (0 — решение завершено успешно); $$K$$ — количество итераций; $$err$$ — погрешность вычислений.
Пример 8.14. Вычислить интеграл $$\int\limits_{0}^1t^2\sqrt{\left(3+sin \left(\frac{1}{t}\right)\right)}dt$$
Обратите внимание, что в нижней границе интегрирования подынтегральная функция терпит разрыв. Решение примера с применением функции $$quad$$ приведено в листинге 8.16.
clear all; function y=f(x) y=(x.^2).-sqrt(3+sin(1./x)); end; >>> format long >>> [F, kod, K, err]=quad(’f’, 0, 1) F = 0.654343719149802 kod = 0 K = 1323 err = 1.37557012147481 e-08 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-05) F = 0.654343738854992 kod = 0 K = 315 err = 3.82733563379833 e -06 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-20) F = 0.654343718970708 kod = 0 K = 1491 err = 9.39557628735834 e-09 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e-20,0.1) F = 0.654343710193938 kod = 0 K = 840 err = 5.80259740257105 e-09 >>> [F, kod, K, err]=quad(’f’, 0, 1, 1.0 e- 20,0.001) F = 0.654343718720156 kod = 0 K = 1596 err = 8.35248716562893 e-09
Функции
$$F = quadl(f, a, b[, tol, trace])$$ и $$[F, err] = quadgk(f, a, b[, tol, trace])$$, где: $$name$$ — имя функции, задающей подынтегральное выражение; $$a, b$$ — пределы интегрирования; $$tol$$ — точность вычислений; $$trace$$ — таблица промежуточных вычислений; $$F$$ — значение интеграла; $$err$$ — погрешность вычислений; также выполняют интегрирование по квадратуре Гаусса. В этих функциях специальным образом подбирается шаг. В первом случае по методу Гаусса-Лобатто, во втором Гаусса-Конрада.
Пример 8.15. Вычислить интеграл $$\int\limits_{-\frac{\pi}{3}}^{\frac{\pi}{3}}tg^4(x)dx$$
Решение примера с применением функций $$quadl$$ и $$quadgk$$ приведено в листинге 8.17.
function y=f(x) y = tan(x).^4; end; >>> format long >>> [F]= quadl (’f’,-pi/3, pi/3, 1.0 e-05) F = 2.09439512983937 >>> [F, err]=quadgk(’f’, -pi/3, pi/3, 1.0 e-05) F = 2.09439510239319 err = 1.02555919485880 e-12
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.