Вычисление определенных интегралов

В MATLAB определены команды quad и quadl для приближенного вычисления определенных интегралов

I = dx.

Команда quad (или quadl) имеет следующие модификации:

quad('f(x)', a,b);

quad('f(x)', a,b,tol).

В выражениях функции приняты следующие обозначения:

'f(x)' – подынтегральная функция или имя файл - функции, ее вычисляющей, взятые в кавычки;

a, b – пределы интегрирования;

tol – относительная погрешность, задаваемая пользователем, по умолчанию tol = 10^-6.

Команды quad (quadl) имеют и другие модификации, о которых можно узнать с помощью команды doc funfun.

Реализованный в quad алгоритм основан на квадратурной формуле Симпсона с автоматическим подбором шага интегрирования для достижения требуемой погрешности вычислений. Интегрирование достаточно гладких функций имеет смысл производить командой quadl, основанной на более точных квадратурных формулах Гаусса - Лобатто, что экономит время вычислений.

В определенном интеграле

I = dx

подынтегральная функция

f(x) =

имеет бесконечный разрыв в точке x = 2, но интеграл сходится при любых конечных значениях a и b. Известна первообразная sign(x-2)·2 от f(x), поэтому точное значение интеграла находим по формуле Ньютона - Лейбница:

I = sign(b-2)·2 - sign(a-2)·2.

В частности,

I1 = dx = 1, I2 = dx = 2, I3 = dx = 4.

Интегралы I1, I2, I3 расположены в порядке убывания степени «гладкости» подынтегральной функции f(x):

I1 – f(x) непрерывна на [1; 1,75];

I2 – точка разрыва x = 2 совпадает с правым концом промежутка интегрирования [1; 2];

I3 – точка разрыва x = 2 лежит внутри промежутка интегрирования [1; 3].

I1 называется собственным интегралом, а I2 и I3 – несобственными интегралами.

Вычислим эти интегралы. Текст файл-функции f, вычисляющей подынтегральное выражение в векторизованном коде, приведен ниже:

function y=f(x)

y=1./sqrt(abs(x-2));

Вычислим интеграл I1 функцией quad:

>> format long

>> tic,I1=quad(@f,1,1.75),toc

I1 =

1.00000005518294

elapsed_time =

0.01000000000000

Напомним, что набор команд tic и toc является средством для приблизительного замера (оценки) быстродействия выполнения оператора или M-функции (см. разд. 1.7).

Возникает вопрос, сколько верных значащих цифр в результе I1=1,00000005518294, т. е. с какой точностью вычислен интеграл? В этом случае ответ очевиден, т. к. точное значение интеграла известно и равно 1.

Приведем результаты вычислений интеграла I1 командами quad и quadl с погрешностями 10^-6, 10^-8, 10^-10, 10^-12, eps ≈ 2,22·10^-16:

quadquadl

10^-6 I1=1,00000005518294, t=0,01 с. I1=1,000000000000057, t=0,01 с.

10^-8 I1=1,00000000037531, t=0,02 с. I1=1,000000000000057, t=0,01 с.

10^-10 I1=1,00000000000121, t=0,04 с. I1=1,000000000000012, t=0,02 с.

10^-12 I1=1,00000000000000, t=0,09 с. I1=1,000000000000000, t=0,04 с.

eps I1=1,00000000000000, t=0,46 с. I1=1,000000000000000, t=0,09 с.

Системная переменная eps (см. разд. 1.2) задает минимальную допустимую погрешность вычислений.

Как видим, с ростом точности вычислений результаты приближаются к точному значению интеграла I1=1. В данном случае quadl оказалась эффективнее quad по точности и быстродействию.

При вычислениях с quad и quadl указатель @f на файл-функцию f можно заменить на 'f' или на '1./sqrt(abs(x-2))', например:

>> tic,I1=quad('1./sqrt(abs(x-2))',1,1.75,eps),toc

I1 =

elapsed_time =

1.71200000000000

При такой форме задания f(x) результат вычислений прежний, но время вычислений увеличилось с 0,46 до 1,712 с.

Вывод: первый аргумент quad или quadl лучше вычислять файл-функцией.

Приведем результаты вычислений интеграла I2:

quadquadl

10^-6 I2=2,00001060926878, t=0,06 с. I2=2,00000041416806, t=0,09 с.

10^-8 I2=2,00000006861257, t=0,13 с. I2=1,99999998793697, t=0,16 с.

10^-10 I2=1,99999998659242, t=0,31 с. I2=1,99999998412995, t=0,31 с.

10^-12 I2=Inf, t=0,81 с. I2=1,99999998563725, t=0,54 с.

eps I2=1,99940376288895, t=2,8 с. I2=2,01543021355643, t=1,76 с.

Точное значение интеграла I2=2. В данном случае сравнить quadl и quad по эффективности сложнее. Отметим, что при вычислениях с меньшей точностью получены более точные результаты (наиболее точные при 10^-10). В одном случае quad не смогла вычислить интеграл.

Вычислим интеграл I3 командами quad и quadl:

>> tic,I3=quad(@f,1,3),toc

I3 =

NaN

elapsed_time =

0.01000000000000

>> tic,I3=quadl(@f,1,3),toc

I3 =

Inf

elapsed_time =

0.01000000000000

Команды quad и quadl не смогли вычислить несобственный интеграл I3.

Вычислить несобственный интеграл с бесконечным верхним пределом

dx = 1.

Первообразная -e^-x(x+1) от подынтегральной функции известна, поэтому точное значение интеграла найдено по формуле Ньютона – Лейбница.

Вычислим интеграл командами quad и quadl:

>> quad('x.*exp(-x)',0,inf)

ans =

NaN

>> quadl('x.*exp(-x)',0,inf)

ans =

NaN

Команды quad и quadl не смогли вычислить интеграл.

В рассмотренных примерах точные значения интегралов были заранее известны. Команды quad и quadl хорошо справились с собственным интегралом I1 и хуже с несобственным интегралом I2. Остальные несобственные интегралы они не смогли вычислить.

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

Вычислить определенный интеграл

I = dx.

Аналитического выражения для первообразной не существует, поэтому точное значение интеграла неизвестно.

Текст файл-функции f1, вычисляющей подынтегральное выражение, приведен ниже:

function f=f1(x)

f=exp(abs(sin(x)));

Вычислим интеграл командой quad:

>> format long

>> tic,I=quad(@f1,0,2*pi),toc

I =

12.41751613152453

elapsed_time =

0.03000000000000

Число верных значащих цифр результа I=12,41751613152453 неизвестно. Приведем результаты вычислений интеграла командами quad и quadl:

quadquadl

10^-6 I=12,41751613152453, t=0,03 с. I=12,41751607142046, t=0,04 с.

10^-8 I=12,41751607216332, t=0,07 с. I=12,41751607142223, t=0,05 с.

10^-10 I=12,41751607142415, t=0,15 с. I=12,41751607142222, t=0,08 с.

10^-12 I=12,41751607142224, t=0,39 с. I=12,41751607142222, t=0,13 с.

eps I=12,41751607142222, t=2,13 с. I=12,41751607142222, t=0,52 с.

Длина мантиссы формата вывода long ограничена 16 значащими цифрами. Поэтому I=12,41751607142222 – точное наблюдаемое значение интеграла для этого формата. Правильность всех цифр подтверждается символьным интегрированием командой int в разделе 7.10. В данном случае quadl оказалась эффективнее quad по точности и быстродействию.

В MATLAB определена команда dblquad для приближенного вычисления двойных интегралов.

Вычислить двойной интеграл

I = e^sin(x²+y²)dxdy.

Файл-функция f2, вычисляющая подынтегральное выражение должна содержать два входных аргумента x и y. Ее текст приведен ниже:

function f=f2(x,y)

f=exp(sin(x.^2+y.^2));

Команда dblquad имеет семь входных аргументов, причем два последних необязательны. При ее вызове необходимо учесть, что первыми задаются пределы внутреннего интеграла по х, а вторыми – внешнего по y:

>> format long

>> tic,I=dblquad(@f2,0,pi,0,pi),toc

I =

13.44629073627405

elapsed_time =

3.08500000000000

Число верных значащих цифр результата I=13,44629073627405 неизвестно. Дополнительным шестым аргументом можно задать погрешность вычислений. При вычислении двойного интеграла dblquad использует quad. Если седьмым аргументом указать 'quadl', то для достаточно гладких функций интеграл будет вычислен быстрее. Приведем результаты вычислений интеграла командой dblquad без 'quadl' и с 'quadl':

без 'quadl' с 'quadl'

10^-6 I=13,44629073627405, t=3,08 с. I=13,44629794896821, t=5,21 с.

10^-8 I=13,44629793673247, t=20,5 с. I=13,44629793729471, t=15,8 с.

10^-10 I=13,44629793749324, t=111 с. I=13,44629793733824, t=80,3 с.

10^-12 I=13,44629793749208, t=719 с. I=13,44629793733828, t=206 с.

eps I=13,44629793733828, t=22972 с. I=13,44629793733828, t=2800 с.

С ростом точности вычислений результаты стремятся к точному значению интеграла I=13,44629793733828. Подтвердить павильность всех цифр интегрированием командой int не удается из - за большого времени вычислений (см. разд. 7.10). В данном случае dblquad с аргументом 'quadl' намного эффективнее int по быстродействию.

В MATLAB определена также команда triplequad для приближенного вычисления тройных интегралов. Вычисления с triplequad аналогичны вычислениям с dblquad.

Вычислить тройной интеграл

I = sin(x²+y²+z²)dxdydz.

Файл-функция f3, вычисляющая подынтегральное выражение должна содержать три входных аргумента x, y и z. Ее текст приведен ниже:

function f=f3(x,y,z)

f=sin(x.^2+y.^2+z.^2);

Вычислим интеграл с погрешностью «по умолчанию»:

>> format long

>> tic,I=triplequad('f3',0,pi,0,pi,0,pi,[],'quadl'),toc

I =

0.28050086212989

elapsed_time =

1.448080000000000e+002

Число верных значащих цифр результата I=0,28050086212989 неизвестно.

Приведем результаты вычислений интеграла командой triplequad с разными погрешностями вычислений:

10^-6 I=0,28050086212989, t=145 с.

10^-7 I=0,28050099966310, t=580 с.

10^-8 I=0,28050100157222, t=1448 с.

10^-9 I=0,28050099359960, t=4165 с.

Время вычислений быстро возрастает с повышением их точности. Результат I=0,2805009936125356921 с 20 верными цифрами будет получен в разделе 7.10 примерно за 2 с. В данном случае int оказалась эффективнее triplequad по точности и быстродействию.

Вывод: точность вычислений собственных интегралов командами quad, quadl, dblquad и triplequad (с погрешностью «по умолчанию») достаточна для практических расчетов; вычисление кратных интегралов с более высокой точностью требует значительных временных затрат.

Date: 2016-08-30; view: 202; Нарушение авторских прав