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

 

Поиск локального минимума функции одной переменной на некотором отрезке осуществляется командой fminbnd. Одна из модификаций fminbnd имеет вид fminbnd('file', x1,x2). Здесь

file – имя файл-функции, вычисляющей значение функции,

x1 и x2 – границы отрезка изоляции локального минимума.

Первый входной аргумент можно задать как указатель на файл-функцию @»file. О других модификациях команды fminbnd можно узнать с помощью команды doc fminbnd.

Найти локальные минимумы функции e^-xcos2πx на отрезке [0;2]. Создадим файл-функцию gr, вычисляющую значение функции e^-xcos2πx при заданном значении аргумента x:

function y=gr(x)

y=exp(-x)*cos(2*pi*x);

Перед нахождением локальных минимумов построим график исследуемой функции командой fplot (рис. 6.3):

>> fplot(@gr,[0,2])

 

 

Рис. 6.3

На рис. 6.3 видно, что исследуемая функция имеет два локальных минимума. Вычислим значение х, при котором достигается второй локальный минимум:

>> x2=fminbnd(@gr,1.4,1.6)

x2 =

1.4749

Итак, второй локальный минимум достигается при х ≈ 1,4779.

Для одновременного вычисления значения функции в точке минимума следует вызвать fminbnd с двумя аргументами:

>> [x2,f]=fminbnd(@gr,1.4,1.6)

x2 =

1.4749

f =

-0.2260

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

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

 

Для поиска минимума функции нескольких переменных применяется команда fminsearch.

Использование команды продемонстрируем на примере нахождения локального минимума функции f(x,y)=sinsin, зависящей от двух переменных x и y. Сначала получим представление о поведении функции, построив ее линии уровня при помощи следующих команд (рис. 6.4):

>> [X,Y]=meshgrid(0:.01:2);

>> Z=sin(pi/2*X.^2).*sin(pi/2*Y.^2);

>> [CMatr,h]=contour(X,Y,Z,[-.9,-.5,-.1,.1,.5,.9]);

>> clabel(CMatr,h)

На графике, приведенном на рис. 6.4, видно расположение локальных минимумов и максимумов. Один из локальных минимумов имеет стартовое приближение (1;1,7) в системе координат xOy.

 

 

Рис. 6.4

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

function f=gr2(arg)

x=arg(1);

y=arg(2);

f=sin(pi/2*x.^2).*sin(pi/2*y.^2);

Теперь для уточнения координат точки минимума вызовем fminsearch с двумя входными аргументами – указателем на файл-функцию gr2 и стартовым приближением (1;1,7). В результате выходной аргумент – вектор M, возвращает координаты искомой точки минимума:

>> M=fminsearch(@gr2,[1,1.7])

M =

1.0000 1.7320

Для получения значения функции в точке минимума следует вызвать fminsearch с двумя выходными аргументами:

>> [M,f]=fminsearch(@gr2,[1,1.7])

M =

1.0000 1.7320

f =

-1.0000

Команды fzero, fsolve, fminbnd и fminsearch позволяют задать дополнительный параметр options, контролирующий вычислительный процесс. Значение options следует предварительно сформировать при помощи команды optimset в соответствии с характером требуемого контроля.

В разделе 6.3 найдено приближенное значение локального минимума функции e^-xcos2πx. Представим результат в формате вывода long:

>> format long

>> [x2,f]=fminbnd(@gr,1.4,1.6)

x2 =

1.47489614357390

f =

-0.22596214557056

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

>> options=optimset('TolX',eps);

>> [x2,f]=fminbnd(@gr,1.4,1.6,options)

x2 =

1.47488038336724

f =

-0.22596214670612

Здесь eps – системная переменная MATLAB (см. разд. 1.2):

>> disp(eps)

2.220446049250313e-016

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

В следующем разделе проводится сравнение результатов вычислений с различной заданной точностью.

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

 

В 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 (с погрешностью «по умолчанию») достаточна для практических расчетов; вычисление кратных интегралов с более высокой точностью требует значительных временных затрат.