Лекция 4. Решение типовых задач алгебры и анализа

Лекция 4
Лекция 4. Решение типовых задач алгебры и анализа.
1. Встроенные средства решения типовых задач алгебры и анализа..............................................	1
1.1. Задачи линейной алгебры. .............................................................................................................	1
1.2. Нахождение нулей функций............................................................................................................	2
1.3. Поиск минимума функций. ..............................................................................................................	5
Функции одной вещественной переменной............................................................................................	5
Функции нескольких переменных. ........................................................................................................	6
1.4. Управление ходом вычислением. .................................................................................................	8
1.5. Вычисление определенных интегралов. .....................................................................................	9
Вычисление определенных интегралов..................................................................................................	9
Вычисление двойных интегралов. .......................................................................................................	10

1. Встроенные средства решения типовых задач алгебры и анализа.

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

1.1. Задачи линейной алгебры.

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

В системе MATLAB для решения довольно сложной задачи линейной алгебры – нахождение корней систем линейных уравнений – предусмотрены знаки операций / и \. Чтобы решить систему линейных уравнений:

a11 x1 + a12 x2 +... + a1n xn = b1

...

an1 x1 + an2 x2 +... +ann xn = bn

которую обычно кратко записывают в виде:

Ax = B

где A – заданная квадратная матрица (размером NxN) коэффициентов при неизвестных aij , a B – заданный вектор-столбец длины N свободных членов bi , достаточно применить операцию \ и вычислить выражение x=A\B.

1

Лекция 4

Операция \, как уже было сказано в первой лекции, называется левым делением матриц и, будучи примененная к матрицам A и B в виде A\B, эквивалентна вычислению выражения inv(A)*B. Здесь под inv(A) понимается вычисление матрицы, обратной к матрице A.

Для вычисления определителя квадратной матрицы А предназначена встроенная функция det(A).

Собственные числаλi и собственные векторы ui ≠ 0 квадратной матрицы А удовлетворяют равенствам A ui = λiui . Функция eig, вызванная с входным аргументом матрицей, находит все собственные числа матрицы и записывает их в выходной аргумент – вектор:

>>A = [2 3; 3 5];

>>lam = eig(A) lam =

0.1459

6.8541

Для одновременного вычисления всех собственных векторов и чисел следует вызвать eig с двумя выходными аргументами.

>> [U, Lam] = eig(A) U =

-0.8507 0.5257 0.5257 0.8507

Lam =

0.1459 0

0 6.8541

Первый выходной аргумент и представляет собой матрицу, столбцы которой являются собственными векторами. Для доступа, например, к первому собственному вектору следует использовать индексацию при помощи двоеточия – U(:,1).

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

Проверьте, правильно ли найдены, например, второе собственное число и соответствующий ему собственный вектор. Воспользуйтесь определением:

>> A*U(:, 2) - Lam(2, 2)*U(:, 2) ans =

1.0e-015 *

0.4441 -0.8882

1.2. Нахождение нулей функций.

Функция fzero позволяет приближенно вычислить корень уравнения на некотором интервале или ближайший к заданному начальному приближению. В простейшем варианте fzero вызывается с двумя входными и одним выходным аргументом x = fzero(func_name, x0), где func_name – имя файл-функции, вычисляющей левую часть уравнения, x0 — начальное приближение к корню, х

– найденное приближенное значение корня. Решим, например, на отрезке		[-5, 5] уравнение:
	sin(x) − x2 cos(x) = 0
Перед	нахождением корней полезно построить график функции, входящей в левую часть
уравнения.	Для получения графика можно прибегнуть к plot, но	все равно понадобится

2

Лекция 4

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

function y = myf(x)

y = sin(x) - x.^2.*cos(x);

Теперь построим график функции myf, используя fplot, и нанесем сетку.

>>fplot('myf', [-5 5])

>>grid on

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

Рис. 1 График функции myf.

Уточним значение корня, расположенного вблизи х = - 5, при помощи fzero:

>> xl = fzero('myf', -5) xl =

-4.7566

Итак, приближенное значение корня равно -4.7566. При указании начального приближения к корню алгоритм fzero автоматически отделяет корень, т.е. вблизи заданного начального приближения находится отрезок, содержащий корень. Проверим ответ, вычислив значение функции myf в точке x1:

>> myf(xl) ans =

2.6645e-015

3

Лекция 4

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

Вместо начального приближения вторым параметром fzero можно задать интервал, на котором следует найти корень:

>> x2 = fzero('myf',[-3 -1]) x2 =

-1.8539

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

В качестве исследуемой функции может выступать и встроенная математическая функция, например:

>> x=fzero('cos', 1) x =

1.5708

Как мы уже говорили, формат вызова функции fzero, имеет вид: x = fzero(func_name, x0)

Необходимо отметить, что имя функции func_name, задающее исходное выражения, для fzero (это также справедливо для большинства других функций MATLAB, таких как fplot, fminbnd и др.) можно задавать как:

•как имя m-файла (в апострофах и без расширения m) – 'MyFunc1', 'cos'

•как указатель на функцию – @MyFunc1, @cos

•или как формула с одним неизвестным x, которую необходимо заключить в апострофах –

'cos(x)+x^3-x', 'cos(x)'

Вычисление всех корней полинома.

Полином в MATLAB задается вектором его коэффициентов, например, для определения полинома p = x7 +3.2x5 −5.2x4 +0.5x2 + x −3 следует использовать команду:

>> p = [1, 0, 3.2, -5.2, 0, 0.5, 1, -3];

Число элементов вектора, т. е. число коэффициентов полинома, всегда на единицу больше его степени, нулевые коэффициенты должны содержаться в векторе.

Функция polyval предназначена для вычисления значения полинома от некоторого аргумента:

>> polyval(p,1) ans =

-2.5000

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

Нахождение всех корней полиномов осуществляется при помощи функции roots, в качестве аргумента которой указывается вектор с коэффициентами полинома. Функция roots возвращает вектор корней полинома, в том числе и комплексных:

4

Лекция 4

>> r = roots(p)

r=

-0.5668 + 2.0698i -0.5668 - 2.0698i -0.6305 + 0.5534i -0.6305 - 0.5534i 1.2149

0.5898 + 0.6435i

0.5898 - 0.6435i

Известно, что число всех корней полинома совпадает с его степенью. Убедимся в правильности работы roots, вычислив значение полинома от вектора его корней:

>> polyval(p,r) ans =

1.0e-012 *

0.8700 + 0.2254i

0.8700 - 0.2254i

0.0031 - 0.0007i

0.0031 + 0.0007i -0.0329

0.0018 + 0.0031i

0.0018 - 0.0031i

В верхней строке результата содержится общий множитель 1.0e-012, на который следует умножить каждую компоненту получившегося вектора.

1.3. Поиск минимума функций.

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

Функции одной вещественной переменной.

Поиск локального минимума функции одной переменной на некотором отрезке осуществляется при помощи fminbnd, использование которой схоже с fzero. Найдем локальные минимумы функции:

ex 2 +sin 3πx

на отрезке [-1.5,1.5]. Требуется предварительно создать соответствующую файл-функцию, назвав ее,

кпримеру ftest, или ввести inline. Приведем пример использования функции inline:

>>ftest=inline('exp(x.^2)+sin(3*pi*x)');

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

>> fplot(ftest, [-1.5 1.5])

Обратите внимание, что для inline-функций апострофы применять не надо. Кроме того, использованием символа @ перед inline-функций также приведет к ошибке. Из графика, приведенного на рис.2, видно, что исследуемая функция имеет четыре локальных минимума.

5

Лекция 4

Рис. 2 График функции ftest.

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

>> x2 = fminbnd(ftest, -0.5, 0) x2 =

-0.1629

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

>>xx = fminbnd(ftest, -1.5, 1.5)

xx=

0.4861

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

>> [x2,f] = fminbnd(ftest, -0.5, 0) x2 =

-0.1629

f =

0.0275

Функции нескольких переменных.

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

Пример. Зададим начальное приближение равное -0.5 для уже упомянутой функции ftest и найдем точку локального минимума:

6

Лекция 4

>> x2 = fminsearch(ftest,-0.5) x2 =

-0.1629 - 0 . 1 6 2 9

Рассмотрим теперь минимизацию функций нескольких переменных на примере функции двух переменных:

f (x, y) = sin πx sin πy

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

[X, Y] = meshgrid(0:0.02:2); Z = sin(pi*X).*sin(pi*Y); surfc(X,Y,Z)

На получившемся графике, приведенном на рис.3, видно расположение локальных минимумов и максимумов.

Рис. 3 График функции f (x, y) = sin πx sin πy .

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

function f = ftest2(v) x = v(1);

y = v(2);

f = sin(pi*x).*sin(pi*y);

7

Лекция 4

Теперь для нахождения локального минимума вызовите fminsearch с двумя входными аргументами – именем файл-функции и начальным приближением и выходным аргументом – вектором с координатами искомой точки минимума:

>> M = fminsearch('ftest2', [1.4, 0.6]) M =

1.5000 0.5000

Решение найдено с точностью 10-4, как по значениям х и у, так и по значению функции. Для получения не только вектора с координатами точки минимума, но и значения функции следует вызвать fminsearch с двумя выходными аргументами:

>> [M,F] = fminsearch('ftest2', [1.4, 0.6]) M =

1.5000 0.5000

F = -1.0000

1.4. Управление ходом вычислением.

Функции fzero, fminbnd и fminsearch допускают определение дополнительных параметров для управления вычислительным процессом и контроля за ним. Формат вызова:

fzero(fname, x0, options) fminbnd(fname, x1, x2, options) fminsearch(fname, x0, options)

Параметры задаются в управляющей структуре, которую мы будем называть options, как в справочной системе MATLAB, хотя имя может быть произвольным. Перед вызовом вычислительных функций следует предварительно сформировать переменную options в соответствии с характером требуемого контроля, воспользовавшись функцией optimset. Переменная options на самом деле является структурой. Это новый тип данных – до сих пор вы имели дело только с числовыми массивами. Это аналог типа sruct в C++, и record в Pascal. В рамках данного курса умение работать со структурами не требуется, можно просто следовать приведенным правилам заполнения, просмотра и использования управляющей структуры options.

Приступим к формированию структуры options на примере минимизации функции одной переменной при помощи fminbnd. Для изменения точности следует сформировать структуру options, указав требуемую точность (например 10-9) по аргументу при помощи параметра TolX:

>>format long

>>options = optimset('TolX', 1.0e-09);

>>x2 = fminbnd(ftest, -0.5, 0, options)

x2 = -0.16289942841268

В общем случае входные аргументы optimset задаются попарно:

options = optimset('Свойство1', Знач1, 'Свойство2', Знач2, ...);

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

Свойство	Значение	Результат			Примечание
Display	'off'	Информация о вычислительном
		процессе не выводится
	'iter'	Выводится	информация	о
		каждом шаге	вычислительного
		процесса

8

Лекция 4

	'final'	Выводится информация только о
		завершении	вычислительного
		процесса
	'notify'	Выводится	предупреждение,
		если процесс	не сходится
		(используется по умолчанию)
MaxFunEvals	Целое число	Максимальное	количество	Только для
		вызовов исследуемой функции		fminbnd
				fminsearch
MaxIter	Целое число	Максимальное	количество	Только для
		итераций	вычислительного	fminbnd
		процесса		fminsearch
TolFun	Положительное	Точность по	функции для	Только для
	вещественное	останова вычислений		fminbnd
	число			fminsearch
TolX	Положительное	Точность по аргументу функции
	вещественное	для останова вычислений (по
	число	умолчанию 10-6)

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

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

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

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

Начнем с примера нахождения значения определенного интеграла

∫1 e−x sin(x)dx

−1

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

function f = fint(x) f = exp(-x).*sin(x);

Теперь для вычисления интеграла вызовите quad, задав первым аргументом ссылку на функцию fint, а вторым и третьим – нижний и верхний пределы интегрирования. В качестве выходного аргумента можно указать имя переменной, в которую следует записать найденное значение:

>>format long

>>I = quad(@fint,-1,1)

I = -0.66349146785310

По умолчанию функция quad вычисляет приближенное значение интеграла с точностью 10-6. Для изменения точности вычислений следует задать дополнительный четвертый аргумент:

9

Лекция 4

>> I = quad(@fint,-1,1,1.0e-07)

I=

-0.66349366574399

Реализованный в quad рекурсивный алгоритм основан на квадратурной формуле Симпсона с автоматическим подбором шага интегрирования для достижения требуемой относительной погрешности.

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

В MATLAB определена функция dblquad для приближенного вычисления двойных интегралов. Как и в случае вычисления определенных интегралов, следует написать файл-функцию для вычисления подынтегрального выражения. Найдем значение интеграла:

∫01∫−ππ (ex sin y +e−x cos y)dxdy

Функция fint2 должна содержать два входных аргумента х и у, ее текст:

function f = fint2(x, y)

f = exp(x).*sin(y) + exp(-x).*cos(y);

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

>> dblquad(@fint2, -pi, pi, 0, 1) ans =

30.0537

Дополнительным шестым параметром можно задать точность вычисления интеграла.

>>format long

>>dblquad(@fint2, -pi, pi, 0, 1, 1.0e-012) ans =

30.05371586580528

10