Урок №16.

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

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

Элементарные средства решения СЛУ

Решение систем линейных уравнений (СЛУ) относится к самой массовой области применения матричных методов, описанных в уроках 10-12. В этом разделе вы найдете ответы на вопросы, каким образом применяются указанные методы и какие дополнительные функции имеет система MATLAB для решения систем линейных уравнений.

Как известно, обычная СЛУ имеет вид:

а11X1, а12,X2..., а1nXn=b1

Здесь а11, а, 2,..., аппкоэффициенты, образующие матрицу А, которые могут иметь действительные или комплексные значения, x1, х2,..., хпнеизвестные, образующие вектор X, и b1, b2,..., bп.свободные члены (действительные или комплексные), образующие вектор В. Эта система может быть представлена в матричном виде как АХ=В, где А — матрица коэффициентов уравнений, X — искомый вектор неизвестных и В — вектор свободных членов. В зависимости от вида матрицы А и ее характерных особенностей MATLAB позволяет реализовать различные методы решения.

Для реализации различных алгоритмов решения СЛУ и связанных с ними матричных операций применяются следующие операторы: +,-,*,/, \, *, ' . Как отмечалось ранее, MATLAB имеет два различных типа арифметических операций — поэлементные и для массивов (векторов и матриц) в целом. Матричные арифметические операции определяются правилами линейной алгебры.

Арифметические операции сложения и вычитания над массивами выполняются поэлементно. Знак точки «.» отличает операции над элементами массивов от матричных операций. Однако, поскольку операции сложения и вычитания одинаковы для матрицы и элементов массива, знаки «.+» и «.-» не используются. Рассмотрим другие операторы и выполняемые ими операции.

 

Для случая нескалярных А и В число столбцов матрицы А должно равняться числу строк матрицы В. Скаляр может умножаться на матрицу любого размера.

Если А — матрица размера пхп, а В — вектор-столбец с п компонентами или матрица с несколькими подобными столбцами, тогда Х=А\В — решение уравнения АХ=В, которое находится хорошо известным методом исключения Гаусса.

Если А — матрица размера тхп и тхп, а В представляет собой вектор-столбец с m компонентами или матрицу с несколькими такими столбцами, тогда система оказывается недоопределенной или переопределенной и решается на основе минимизации второй нормы невязок.

Ранг k матрицы А находится на основе QR-разложения (урок 11) с выбором ведущего элемента. Полученное решение X будет иметь не больше чем k ненулевых компонентов на столбец. Если k<n, то решение, как правило, не будет совпадать с pinv(A)*B, которое имеет наименьшую норму невязок ||Х||.

При записи СЛУ в матричной форме необходимо следить за правильностью записи матрицы А и вектора В. Пример (в виде m-файла):

А-[2  1

0

1:

1    -3

2

4;

-5    0

-1

-7:

1    -6

2

6]:

В=[8  9

-5

0]:

Х1=В/А



Х2=В*А^-1



X3=B*inv(A)



Эта программа выдает результаты решения тремя способами:

X1 =

3.0000 -4.0000-1.00001.0000

Х2 =

3.0000 -4.0000-1.00001.0000

X3 =

3.0000 -4.0000-1.00001.0000

Как и следовало ожидать, результаты оказываются одинаковыми для всех трех методов. При решении систем линейных уравнений, особенно с разреженной матрицей коэффициентов, полезно применение функций colmmd (colamd), symmmd (symamd), описанных ранее в уроке 12.

Функции для решения систем линейных уравнений с ограничениями

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

mse = B'*(inv(V)-inv(V)*A*inv(A'*inv(V)*A)*A'*inv(V))*B./(m-n) 

dX = sqrt(diag(inv(A'*inv(V)*A)*mse))

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

» С=[4 3:12 5:3 12];b=[1.45,41];D=lsqnonneg(C.b') 

D =

2.2242

2.6954

Решение СЛУ с разреженными матрицами

Решение СЛУ с разреженными матрицами — хотя и не единственная, но несомненно одна из основных сфер применения аппарата разреженных матриц, описанного в уроке 12. Ниже приведены функции, относящиеся к этой области применения разреженных матриц. Большинство описанных методов относится к итерационным, т. е. к тем, решение которых получается в ходе ряда шагов — итераций, постепенно ведущих к получению результата с заданной погрешностью или с максимальным правдоподобием [33]. Описанные ниже функции MATLAB могут и должны применяться и при решении обычных СЛУ - без разреженных матриц.[Использование всех этих функций, кроме Isqr, как правило, ограничено системами уравнений, в которых А — нормальная квадратная матрица, т. е. А* А -АА*, где А*— сопряженная (эрмитова) матрица А. — Примеч. ред.]

Точное решение, метод наименьших квадратов и сопряженных градиентов

Если значение flag больше нуля, то возвращается не последнее решение, а то решение, которое имеет минимальное значение отношения вторых норм векторов norm(B-A*x)/norm(B).

Пример:

» А=[0 012; 1300; 0101; 1010];

» В=[11; 7; 6; 4];

Введенные в этом примере матрица А и вектор В будут использованы и в других примерах данного раздела. В примере процесс итераций сходится на пятом шаге с относительным остатком (отношением вторых норм векторов невязки и свободных членов) 1.9 10-13.

Пример:

» lsqr(A,B.le-6.5)

Isqr converged at iteration 5 to a solution

with relative residual 

1.9e-013 

ans =

1.0000 

2.0000 

3.0000

4.0000

Двунаправленный метод сопряженных градиентов

Решение СЛУ с разреженной матрицей возможно также известным двунаправленным методом сопряженных градиентов. Он реализован указанной ниже функцией.

Пример:

 » bicg(A.B)

BICG converged at iteration 4 

to a solution with relative residual 

2.3e-015

ans= 

1.0000 

2.0000 

3.0000 

4.0000

Устойчивый двунаправленный метод

Еще один двунаправленный метод, называемый- устойчивым, реализует функция bicgstab:

» bicgstab(A.B)

BICGSTAB converged at iteration 3.5 

to a solution with relative residual 

6.5e-012

ans = 

1.0000 

2.0000 

3.0000 

4.0000

Метод сопряженных градиентов

Итерационный метод сопряженных градиентов реализован функцией peg: О рсд(А.В) — возвращает решение X СЛУ А*Х=В. Матрица А должна быть квадратной, симметрической [В нашем примере матрица А — несимметрическая, т. е. A(i,j)—*A(j,i). — Примеч. ред. ]  и положительно определенной [Матрица называется положительно определенной, если все ее собственные значения (характеристические числа) действительные и положительные. — Примеч. ред.]. Функция pcg начинает итерации от начальной оценки, представляющей собой вектор размером п, состоящий из нулей. Итерации производятся либо до сходимости решения, либо до появления ошибки, либо до достижения максимального числа итераций. Сходимость достигается, если относительный остаток norm(b-A*X)/norm(B) меньше или равен погрешности метода (по умолчанию 1е-6). Максимальное число итераций — минимум из п и 20. Функция pcg(...) имеет и ряд других форм записи, описанных для функции bieg(...). Пример:

» pcg(A.B)

Warning: PCG stopped after the maximum 4 iterations 

without converging to the desired tolerance le-006

The iterate returned (number 4) has relative residual 0.46 

> In C:\MATI_AB\toolbox\matlab\sparfun\pcg.m at line 347 

ans =

1.7006

1.2870

-2.0535

8.2912

В данном случае решение к успеху не привело, поскольку матрица А —несимметрическая. Новая функция rrrinres не требует, чтобы матрица А была положительно определенной. Достаточно, чтобы она была квадратной и симметрической. В отличие от peg минимизируется не относительная невязка, а абсолютная. Но и эта функция не может решить наш пример:

» minres(A.В.1е-6.1000000)

minres stopped at iteration 1000000 without converging 

to the desired tolerance le-006 

because the maximum number of iterations was reached. 

The iterate returned (number 1000000) has relative residual 0.011 

ans =

1.9669

3.7757

3.0789

1.9367

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

Квадратичный метод сопряженных градиентов

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

» cgs(A.B)

CGS converged at iteration 4 to a solution

with relative residual 4e-014 

ans =

1.0000

2.0000

3.0000

4.0000

Метод минимизации обобщенной невязки

Итерационный метод минимизации обобщенной невязки также реализован в системе MATLAB. Для этого используется функция gmres:

» gmres(A.B)

GMRES(4) converged at Iteration 1(4) to a solution with relative residual le-016

ans =

1.0000

2.0000

3.0000

4.0000

Квазиминимизация невязки — функция qmr

Метод решения СЛУ с квазиминимизацией невязки реализует функция qmr:

» qmr(A.B)

QMR converged at iteration 4 to a solution

with relative residual l.le-014 

ans =

1.0000

2.0000

3 .0000

4.0000

Вычисление нулей функции одной переменной

Ряд функций системы MATLAB предназначен для работы с функциями. По аналогии с дескрипторами графических объектов могут использоваться объекты класса дескрипторов функций, задаваемых с помощью символа @, например: » fe=@exp.

Подфункциями понимаются как встроенные функции, например sin(x) или ехр(х),так и функции пользователя, например f(x), задаваемые как т-файлы-функции.

Численные значения таких функций, заданных дескрипторами, вычисляются с помощью функции feval:

» feval(fe.1.0)

ans =

2.7183

Для совместимости с прежними версиями можно записывать функции в символьном виде в апострофах, использование функции eval для их вычисления может быть более наглядно, не нужно создавать m-файл, но в учебном курсе мы будем стараться использовать новую нотацию, с использованием дескрипторов функций и feval, так как при этом программирование становится «более объектно-ориентированным», повышается скорость, точность и надежность численных методов. Поэтому, хотя везде в нижеследующем тексте вместо @fun можно подставить и символьное значение функции в апострофах, мы будем использовать нотацию @fun в дидактических целях. Все же иногда в интерактивном режиме можно использовать старую запись, чтобы не создавать m-файл функции.

Довольно часто возникает задача решения нелинейного уравнения вида f(x) = О или/, (г) =/2(дг). Последнее, однако, можно свести к виду f(x) =f1(х) - f2(х) = 0. Таким образом, данная задача сводится к нахождению значений аргумента х функции f(x) одной переменной, при котором значение функции равно нулю. Соответствующая функция MATLAB, решающая данную задачу, приведена ниже:

fzero(fun.x.[ ].[ ].Р1).

Для функции fzero ноль рассматривается как точка, где график функции fun пересекает ось х, а не касается ее. В зависимости от формы задания функции fzero реализуются следующие хорошо известные численные методы поиска нуля функции: деления отрезка пополам, секущей и обратной квадратичной интерполяции. Приведенный ниже пример показывает приближенное вычисление р/2 из решения уравнения cos(x)=0 с представлением косинуса дескриптором:

» х= fzero(@cos.[l 3]) 

x =

1.5708

В более сложных случаях настоятельно рекомендуется строить график функции f(x) для приближенного определения корней и интервалов, в пределах которых они находятся. Ниже дан пример такого рода (следующий листинг представляет собой содержимое m-файла fun1.m):

%Функция, корни которой ищутся 

function f=funl(x) 

f=0.25*x+sin(x)-l;

» х-0:0.1:10;

» plot(x.funl(x));grid on;

Из рисунка нетрудно заметить, что значения корней заключены в интервалах [0.5 1], [2 3] и [5 6]. Найдем их, используя функцию fzero:

» xl=fzero(@funl.[0.5 1]) 

xl =

0.8905

» x2=fzero(@funl.[2 3]) 

x2 =

2.8500

» x3=fzero(@funl,[5.6]) 

x3 =

5.8128

» x3=fzero(@funl,5.0.001)

 x3 =

5.8111

Обратите внимание на то, что корень хЗ найден двумя способами и что его значения в третьем знаке после десятичной точки отличаются в пределах заданной погрешности tol =0.001. К сожалению, сразу найти все корни функция fzero не в состоянии. Решим эту же систему при помощи функции f sol ve из пакета Optimization Toolbox, которая решает систему нелинейных уравнений вида f(x)=0 методом наименьших квадратов, ищет не только точки пересечения, но и точки касания, f solve имеет почти те же параметры (дополнительный параметр — задание якобиана) и почти ту же запись, что и функция lsqnonneg, подробно рассмотренная ранее. Пример:

»fsolve(@funl,0:10 ) 

ans =

Columns 1 through 7

0.8905 0.8905 2.8500 2.8500 2.8500 5.8128 5.8128

Columns 8 through 11

5.8128 2.8500 2.8500 10.7429

Для решения систем нелинейных уравнений следует также использовать функцию solve из пакета Symbolic Math Toolbox. Эта функция способна выдавать результат в символьной форме, а если такого нет, то она позволяет получить решение в численном виде. Пример:

» solve('0.25*x + sin(x) -1) [Пакеты расширения Symbolic Math ToolBox и Extended Symbolic Math Toolbox MATLAB 6.0 используют ядро Maple V Release 5 [30-35] и являются поэтому исключением: они пока не поддерживают новую нотацию с использованием дескрипторов функций. — Примеч. ред.]

ans = 

.89048708074438001001103173059554

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

Еще одна важная задача численных методов — поиск минимума функции f(x) в некотором интервале изменения х — от х1 до х2. Если нужно найти максимум такой функции, то достаточно поставить знак «минус» перед функцией. Для решения этой задачи используется следующая функция:

[X.fval.exitflag,output] = fminbnd(@fun.xl,x2.options, pl,p2,...)

В этих представлениях используются следующие обозначения: xl. х2 — интервал, на котором ищется минимум функции; Р1.Р2... — дополнительные, помимо х, передаваемые в функцию аргументы; fun — строка, содержащая название функции,  которая будет минимизирована; options — вектор параметров вычислений. В зависимости от формы задания функции fminbnd вычисление минимума выполняется известными методами золотого сечения или параболической интерполяции. Пример:

» options=optimset('tolX',1.е-10):

[x]=fminbnd(@cos.3,4,options) 

х =

3.1416

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

Значительно сложнее задача минимизации функций нескольких переменных f(х1,...). При этом значения переменных представляются вектором х, причем начальные значения задаются вектором х0 Для минимизации функций ряда переменных MATLAB обычно использует разновидности симплекс-метода Нелдера-Мида.

Этот метод является одним из лучших прямых методов минимизации функций ряда переменных, не требующим вычисления градиента или производных функции. Он сводится к построению симплекса в n-мерном пространстве, заданного n+1 вершиной. В двумерном пространстве симплекс является треугольником, а в трехмерном — пирамидой. На каждом шаге итераций выбирается новая точка решения внутри или вблизи симплекса. Она сравнивается с одной из вершин симплекса. Ближайшая к этой точке вершина симплекса обычно заменяется этой точкой. Таким образом, симплекс перестраивается и обычно позволяет найти новое, более точное положение точки решения. Решение повторяется, пока размеры симплекса по всем переменным не станут меньше заданной погрешности решения.

Реализующая симплекс-методы Нелдера-Мида функция записывается в виде:

Классическим примером применения функции fminsearch является поиск минимума тестовой функции Розенброка, точка минимума которой находится в «овраге» с «плоским дном»: rb(x1,x2,а) = 100*(x2 - x1)2 + (а - x1 )2.

Минимальное значение этой функции равно нулю и достигается в точке [а а2]. В качестве примера уточним значения x1 и х2 в точке [-1.2 1]. Зададим функцию (в файле rb.m):

% Тестовая функция Розенброка 

function f=rb(x.a) 

if nargin<2 a=l: end 

f=100*(x(2)-x(i^2)^2+(a-x(l)^2:

Теперь решим поставленную задачу:

»options=optimset( 'tolX',l.e-6):

[xmin. opt, rosexflag, rosout]=fminsearch(@rb.[-1.2 1],options) 

xmin =

1.0000 1.0000 

opt =

4.1940e-014 

rosexflag =

1 rosout =

iterations: 101

funcCount: 189

algorithm: 'Nelder-Mead simplex direct search' .

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

Для минимизации функций нескольких переменных можно использовать также функцию MATLAB fminunc и функцию Isqnonlin из пакета Optimization Toolbox. Первая из них позволяет использовать предварительно заданные командой optimset порог сходимости для значения целевой функции, вектор градиентов opt ions, grad-obj, матрицу Гесса, функцию умножения матрицы Гесса или график разреженности матрицы Гесса целевой функции. Isqnonl in реализует метод наименьших квадратов и, как правило, дает наименьшее число итераций при минимизации. Ограничимся приведением примеров их применения для минимизации функции Розенброка:

» options=optimset('tolX',le-6.'To!Fun'.le-6);

» [xmin. opt. exflag. out, grad, hessian ]=fminunc(@rb,[-1.2 1].options)

Warning: Gradient must be provided for trust-region method;

using line-search method instead. 

> In C:\MATLABR12\toolbox\optim\fminunc.m at line 211 

Optimization terminated successfully:

Current search direction is a descent direction, and magnitude of directional derivative in search direction less than 2*options.TolFun 

xmin =

1.0000 1.0000

opt =

1.9116e-011 

exflag=

out =

iterations: 26

funcCount: 162

stepsize: 1.2992

firstorderopt: 5.0020e-004

algorithm: 'medium-scale: Quasi-Newton line search' 

grad=

l.Oe-003 *

-0.5002

-0.1888 

hessian =

820.4028 -409.5496

-409.5496 204.7720

firstorderopt - мера оптимальности для первой нормы градиента целевой функции в найденной точке минимума; 

»options=optimset('tolX' Gе-6. 'maxFunEvals' .162): 

» [xmin. opt]=lsqnonlin(@rb,[-1.2 1].[0 le-6].[0 le-6],options) 

Warning: Large-scale method requires at least as many equations as variables:

switching to line-search method instead. Upper and lower bounds will be ignored.

> In C:\MATLABR12\toolbox\optim\private\lsqncommon.m at line 155 

In C:\MATLABR12\toolbox\optim\lsqnonlin.m at line 121 

Maximum number of function evaluations exceeded Increase 

OPTIONS.maxFunEvals 

xmin =

0.6120 0.3715 

opt =

0.1446

Аппроксимация производных

Аппроксимация Лапласиана

Для выполнения аппроксимации Лапласиана в MATLAB используется следующая функция:

Матрица L имеет тот же размер, что и матрица U, и каждый ее элемент равен разности элемента массива U и среднего значения четырех его соседних элементов (для узлов сетки во внутренней области). Для вычислений используется пятиточечная формула аппроксимации Лапласиана.

Другие формы этой функции также возвращают матрицу L, но допускают дополнительные установки:

» [х.у]= meshgrid(-5:5.-4:4); 

» U =x.*x+y.*y

 U =

41

32

25

20

17

16

17

20

25

32

41

34

25

18

13

10

9

10

13

18

25

34

29

20

13

8

5

4

5

8

13

20

29

26

17

10

5

2

1

2

5

10

17

26

25

16

9

4

1

0

1

4

9

16

25

26

17

10

5

2

1

2

5

10

17

26

29

20

13

8

5

4

5

8

13 .

20

29

34

41

25

 32

18 

25

13 

20

10

 17

9

 16

10 

17

13

20

18 

25

25 

32

34

41

» V=del2(U)

V =

1

1

1

1   1

1   1   1

1

1

1

1

1

1

1   1

1   1  1

1

1

1

1

1

1

1   1

1  1   1

1

1

1

1

1

1

1   1

1   1   1

1

1

1

1

1

1

1   1

1   1   1

1

1

1

1

1

1

1   1

1   1   1

1

1

1

1

1

1

1   1

1   1  1

1

1

1

1

1

1

1   1

1   1   1 

1


1

1

1

1

1   1

1   1  1

1

1

1

» subplot(l,2,l)

 » surfl(U) 

» subplot(l,2,2) 

» surfl(V)

На рис. 16.1 представлены графики поверхностей U и V.

Рис. 16.1. Графики функций U и V

Апроксимация производных конечными разностями

Используя функцию diff, можно строить графики производных заданной функции. Пример этого показан ниже:

» Х=0:0. 05:10: 

» S=sin(X); 

» D=diff(S): 

» plot(D/0.05)

Для получения приближенных численных значений производной от функции sin(.r) вектор конечно-разностных значений D поделен на шаг точек по х. Как и следовало ожидать, полученный график очень близок к графику функции косинуса (рис. 16.2). Обратите внимание, что по оси х отложены номера элемента* вектора X, а не истинные значения х.

Пакет расширения Symbolic Math Toolbox позволяет выполнять дифференциро вание функций в аналитическом виде, т. е. точно. Это, однако, не всегда нужно

Рис. 16.2. Приближенный график производной от функции sin(x)

Вычисление градиента функции

Вычисление конечно-разностным методом градиента функций реализуется следующей функцией:

Пример:

» F=[l 35795678] 

F =

135795678

 » FX = gradient(F) 

FX =

Columns 1 through 7

2.0000 2.0000 2.0000 2.0000 -1.0000-1.50001.0000

Columns 8 through 9

1.0000 1.0000

» F=[l 2 3 6 7 8:1 4 5 7 9 3;5 9 5 2 5 7] 

F =

123678

145793

595257

 » [FX.FY] = gradient(F) 

FX =

1.0000 1.0000 2.0000 2.0000 1.0000 1.0000

.3.0000 2.0000 1.5000 2.0000 -2.0000 -6.0000

4.0000 0 -3.50000 2.5000 2.0000 

FY =

0 2.0000 2.0000 1.0000 2.0000-5.0000

2.0000 3.5000 1.0000 -2.0000-1.0000-0.5000

4.0000 5.0000 0 -5.0000-4.0000 4.0000

Функция gradient часто используется для построения графиков полей градиентов.

Численное интегрирование

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


одним из многочисленных численных методов, для вычисления неопределенного интеграла по рекурсивному алгоритму усредняются значения b=a+dn, где d — предельно малая величина.

Метод трапеций

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

»Y=[1.2.3.4] 

Y =

1 2 3 4 

» trapz(y) 

ans =

7.5000

» X=0:pi/70:pi/2; 

» Y=cos(X); 

» Z = trapz(Y) 

Z =

22.2780

» cumtrapz(y) 

ans=

0 1.5000 4.0000 7.5000 

» Y=magic(4) 

Y =

162 3 13

5 11 10 8

97 6 12

4 14 15 1

» Z= cumtrapz(Y.l) 

Z =

0 0 0 0

10.5000 6.5000 6.5000 10.5000

17.5000 15.500014.500020.5000

24.0000 26.000025.000027.0000 

Численное интегрирование методом квадратур

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

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

Функции quad и quadl используют два различных алгоритма квадратуры для вычисления определенного интеграла. Функция quad выполняет интегрирование по методу низкого порядка, используя рекурсивное правило Симпсона [4, 40]. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB 6 изменен по сравнению с предшествовавшими версиями, точность по умолчанию по сравнению с версиями 5.3х повышена в 1000 раз (с 10-3 до 10-6). Новая функция quadl (квадратура Лобатто) использует адаптивное правило квадратуры Гаусса— Лобатто очень высокого порядка. Устаревшая функция quads выполняла интегрирование, используя квадратурные формулы Ньютона—Котеса 8-го порядка [40]. Достижимая точность интегрирования гладких функций в MATLAB 6 поэтому также значительно выше, чем в предшествующих версиях.

» quad('(exp(x)-l)'.0.1.1e-5) 

ans =

0.7183 

» q = quad(@exp.0.2.1e-4)

q = 

6.3891 

» q = quad(@sin.0,pi,le-3)

q = 2.0000

Пример: пусть m-файл integl.m описывает функцию 2*y*sin(x)+x/2*cos(y), тогда вычислить двойной интеграл от той функции можно следующим образом:

» result = dblquad(@integl,pi,2*pi,0.2*pi) 

result =

78.9574

Работа с полиномами

Полиномы (у нас их принято называть также степенными многочленами) — широко известный объект математических вычислений и обработки данных. Обычно полином записывается в виде

р(х) = апх^n + xп-1x^n -1+ ... + а2x^2 + а1^х + а0,

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

р(х) = a1x^n + а2x^n-1 + ... + апх + ап+1

Последняя запись обычно принята в MATLAB для n, по умолчанию положительных.

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

Умножение и деление полиномов

Ниже приведены функции, осуществляющие умножение и деление полиномов, или, что то же самое, свертку двух входных векторов, в которых находятся коэффициенты полиномов, и операцию, обратную свертке.

Пример:

» f=[2.3.5.6];d=[7,8,3]:r=conv(f,d)

r =

14 37 65 91 63 18

Пример:

» t=[14,37.65.91,63,18]:r=[7.8.3];[w.e]=deconv(t.r) 

w =

2.0000 3.0000 5.0000 6.0000 

е =

1.0е-013

0 0 0.1421 -0.1421-0.2132-0.1066

Вычисление полиномов

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

А =

2 3 6

3 8 6

1 7 4 

» d=poly(A) 

d =

1.0000 -14.0000 -1.0000-40.0000 

» А=[3,6.8:12.23.5:11.12.32] 

А =

3 6 8

1223 5

1112 32

 » poly(A) 

ans =

1.0000 -58.0000 681.0000 818.0000

Приведенная ниже функция вычисляет корни (в том числе комплексные) для полинома вида

Вектор-строка с содержит коэффициенты полинома, упорядоченные по убыванию степеней. Если с имеет n+1 компонентов, то полином, представленный этим вектором, имеет вид . Пример:

» x=[7.45.12.23];d=roots(x) 

d =

-6.2382

-0.0952+0.7195i

-0.0952 -0.7195i

А=[-6.2382 -0.0952+0.71951 -0.0952 -0.71951]: 

B=Poly (А)

В=[1.0000 6.4286 1.7145 3.2859] 

В*7 

ans =

7.0000 45.000212.001523.0013

С погрешностью округления получили тот же вектор.

Пример:

» р=[3,0.4.3]; d=polyval(p,[2,6]) 

d =

35 675

Пример:

» D=pascal(5)

D =




1 1

1

1

1

1 2

3

4

5

1 3

6

10

15

1 4

10

20

35

1 5

15

35

70

f=poly(d)

f =

1.0000 -99.0000 626.0000 -626.0000 99.0000-1.0000 

» polyvalm(f.D) 

ans =

l.0e-006*

-0.0003 -0.0011-0.0038-0.0059-0.0162

-0.0012 -0.0048-0.0163-0.0253-0.0692

-0.0034 -0.0131 -0.0447 -0.0696 -0.1897

-0.0076 -0.0288-0.0983-0.1529-0.4169 

-0.0145-0.0551-0.1883-0.2929-0.7984

Данный пример иллюстрирует также погрешности численных методов, поскольку точное решение дает нулевую матрицу.

Вычисление производной полинома

Ниже приведена функция, возвращающая производную полинома р:

Примеры:

» a=[3.5.8]:b=[5,3.8]:dp=polyder(a) 

dp =

6 5

» dt=polyder(a,b) 

dt =

60102 158 64 » [q.p]=polyder(b.a)

q =

1632 -16

p =

9 30 73 80 64

Решение полиномиальных матричных уравнений

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

где степень полинома р — целое неотрицательное число, а А0, А1,..., Аpвходные матрицы порядка п. Выходная матрица X размера nхnр содержит собственные векторы в столбцах. Вектор е размером пр содержит собственные значения.

Пример:

» А=[1:4:5:8:9:12:13:16] 

А =

1 2

3 4


5 6

7 8


9 10

11 12


1314

15 16


» В=[4:7

;2:5;10:13;23:26]

3 -



4 5

6 7


2 3

4 5


1011

12 13


2324

25 26


» [F.a]=

polyeig(A.B)


F =



0.4373

0.0689 

-0.5426

-0.7594

-0.3372

-0.4969 

0.6675

-0.1314

-0.6375

0.7870 

0.2927

-0.1314

0.5374

-0.3591

- 0.4176

 0.3771

a =



4.4048



0.4425



-0.3229



-1.0000




Разложение на простые дроби

Для отношения полиномов b и а используются следующие функции: 

Пример:

» b=[4.3.1]:a=[1.3.7.1]:[r.p,k]=residue(b,a)

r=

1.9484 + 0.80641 

1.9484 - 0.80641 

0.1033

Р =

-1.4239 + 2.13051

-1.4239 - 2.13051

-0.1523 

k =

[ ]

» [bl,al]=residue(r,p.k) 

b1=

4.0000 3.0000 1.0000 

a1 =

1.0000 3.0000 7.0000 1.0000

Решение обыкновенных дифференциальных уравнений

Анализ поведения многих систем и устройств в динамике, а также решение многих задач в теории колебаний и в поведении упругих: оболочек обычно базируются на решении систем обыкновенных дифференциальных уравнений (ОДУ). Их, как правило, представляют в виде системы из дифференциальных уравнений первого порядка в форме Коши:

с граничными условиями y(t0 tend, p)=b, где tend, t0 — начальные и конечные точки интервалов. Параметр t не обязательно означает время, хотя чаще всего решение дифференциальных уравнений ищется во временной области. Вектор b задает начальные и конечные условия.

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

Решатели ОДУ

Для решения систем ОДУ в MATLAB реализованы различные методы. Их реализации названы решателями ОДУ.

В этом разделе обобщенное название sol ver (решатель) означает один из возможных  численных методов решения ОДУ: ode45, ode23, odel!3, odelSs, ode23s, ode23t , ode23tb, bvp4c или pdepe.

Решатели реализуют следующие методы решения систем дифференциальных уравнений, причем для решения жестких систем уравнений рекомендуется использовать только специальные решатели ode 15s , ode23s, ode23t. ode23tb:

Несмотря на сравнительно низкую точность, этот метод может оказаться более эффективным, чем ode15s;

Все решатели могут решать системы уравнений явного вида у' = F(£, y). Решатели ode15s и ode23t способны найти корни дифференциально-алгебраических уравнений M(t)y' = F(t, у},, где М называется матрицей массы. Решатели ode!5s, ode23s, ode23t и ode23tb могут решать уравнения неявного вида M(t,y) у' = F(t, у). И наконец, все решатели, за исключением ode23s, который требует постоянства матрицы массы, и bvp4c, могут находить корни матричного уравнения вида M(t, у) у' - F(t, у). ode23tb, ode23s служат для решения жестких дифференциальных уравнений . ode15s -жестких дифференциальных и дифференциально-алгебраических уравнений, ode23t -умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

Использование решателей систем ОДУ

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

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

[T.X.Y] - sim(@model....).

Параметры интегрирования (options) могут быть определены и в m-файле, и в командной строке с помощью команды odeset. Если параметр определен в обоих местах, определение в командной строке имеет приоритет.

Решатели используют в списке параметров различные параметры:

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

Параметры

Ode45

Ode23

OdellS

OdelSs

ode23s

RelTol,AbsTol

+

+

+

+

+

OutputFcaOutputSel, Refine, Stats

+

+

+

+

+

Events

+

+

+

+

+

MaxStep, InitlalStep

+

+

+

+

+

Jconstant, Jacobl an,






Jpattern, Vectorized

- ,

-

-

+

+

Mass

- •'

-

-

+

+

MassConstant

-

-

-

+

-

MaxOrder, BOF

-


-

+

-

Решатель bvp4c имеет очень небольшое число параметров, но можно вводить не только матрицу Якоби интегрируемой функции, но и матрицу Якоби, содержащую частные производные функции граничных условий по границам интервала и по неизвестным параметрам.

Покажем применение решателя ОДУ на ставшем классическом примере — решении уравнения Ван-дер-Поля, записанного в виде системы из двух дифференциальных уравнений:

y'1= y2 ;

y'2=100*(1-y1)^2 *y2-y1

при начальных условиях

y1,(0) = 0; 

y2(0) = 1.

Перед решением нужно записать систему дифференциальных уравнений в виде ode-функции. Для этого в главном меню выберем File > New > M-File и введем

function dydt = vdp100(t.y)

dydt = zeros(2.1); % a column vector

dydt(l) = y(2);

dydtC2) = 100*(1 -у(1^)2)*у(2) -y(1);

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

» [T,Y]=odel5s(@vdplOO.[0 30].[2 0]); 

» plot(T.Y)

» hold on:gtext('yl').gtext('y2')

Последние команды позволяют с помощью мыши нанести на графики решений y1 = y(1) и у2 = y(2) помечающие их надписи.

Рис. 16.5. Пример решения системы дифференциальных уравнений численным методом

Описание системы ОДУ

Можно использовать m-файл типа odefunction (или m-file типа odefile для совместимости с прежними версиями, но последний случай мы рассматривать не будем, чтобы определить систему дифференциальных уравнений в одной из явных (первая формула) или неявных форм:

y'= F(t, у), My' = F(t, у) или M(t)y' = Y(t, у),

где t — независимая переменная (скаляр), которая обычно представляет время; у — вектор зависимых переменных; F — функция от t и у, возвращающая вектор-столбец такой же длины как и у; М и М(£) — матрицы, которые не должны быть вырожденными. М может быть и константой.

Рассмотрим пример решения уравнения вида

 

Оно сводится к следующей системе уравнений:


Подготовим m-файл ode-функции vdp.m:

function [outl.out2.out3] = vdp(t.y.flag)

if nargin < 3 | isempty(flag)

outl = [2.*y(2).*(l-y(2).^2)-y(1); y(1)];

else

switch(flag)

case 'inlt' % Return tspan. y0 and options

out1 = [0 20];

out2 = [2; 0];

out3 = [ ];

otherwise

error([' Unknown request ''' flag '''.']);

end

end

Тогда решение системы с помощью решателя ode23 реализуется следующими командами:

» [T.Y] = ode23(@vdp.[0 20]. [2 0]);

Еще проще работать с решателями нового поколения. Рассмотрим систему уравнений: y'+abs(y)=0; y(0)=0; у(4)=-2.

Для решения в пределах отрезка [0; 4] с помощью bvp4c достаточно привести эту систему к виду: y'=-abs(y), y(0)=0; у(4)+2=0. Затем -создаем две ode-функции: twoode и twobc в разных m-файлах:

function dydx = twoode(x,у) 

dydx = [ у(2)

-abs(yd))];

function res = twobc(ya.yb) res = [ ya(l)

yb(l) + 2];

Теперь наберите в командной строке type twobvp и посмотрите само решение уравнения, которое содержится в уже имеющемся в системе файле twobvp. А исполнив команду twodvp, можно наблюдать результат решения в виде графиков. В решении вы найдете структуру узлов начальной сетки решения, которая поясняется ниже.

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

Дескрипторная поддержка параметров решателя

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

Пример:

» options = odesetCRelTol' ,[le-6 le-7].'AbsTol' ,6е-3); 

» odeget(options.'Rel') 

ans =

l.0e-006*

1.0000 0.1000 

» odegetCoptions.'Abs') 

ans =

0.0060

oldopts

F 1 [ ] 4 'S' 'S' [ ] [ ] [ ]

newopts

Т 3 F [] 'S'  [] [] [] []

odeset(oldopts.newopts)

Т 3 F . 4 ' ' 's' [ ] [ ] [ ]

Функция odeset без параметров возвращает все имена свойств и их допустимые значения.

Пример:

» odeset

AbsTol: [ positive scalar or vector {le-6} ]

RelTol: [ positive scalar {le-3} ]

NormControl: [ on | {off} ]

OutputFcn: [ function ]

OutputSel: [ vector of integers ]

Refine: [ positive integer ]

Stats: [ on | {off} ]

InitialStep: [ positive scalar ]

MaxStep: [ positive scalar ]

BDF: [ on | {off} ]

MaxOrder: [ 1 | 2 | 3 | 4 | {5} ]

Jacobian: [ matrix | function ]

JPattern: [ sparse matrix ]

Vectorized: [ on | {off} ]

Mass: [ matrix | function ] MStateOependence: [ none | weak | strong ]

MvPattern: [ sparse matrix ]

MassSingular: [ yes | no | {maybe} ]

InitialSlope: [ vector ] Events: [ function ]

Пакет Partial Differential Equations Toolbox

Специалистов в области численных методов решения дифференциальных уравнений в частных производных несомненно заинтересует обширный пакет Partial Differential Equations Toolbox (PDETB). Хотя этот пакет является самостоятельным приложением и в ядро MATLAB не входит, мы приведем краткое описание некоторых его возможностей с парой примеров. Вы можете вызвать пакет с его графическим интерфейсом командой pdetool.

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

[p,e,t]=initmesh('Ishapeg'); 

u=assempde('Ishapeb'.p.e.t.1.0,1): 

pdeplottp.e.t.'xydata',u.'zdata'.u.'mesh'.'off'):

В состав пакета входит ряд полезных демонстрационных примеров с именами от pdedemol до pdedemo8. Их можно запустить как из командной строки (путем указания имени), так из окна демонстрационных примеров Demos. Рекомендуется, однако, сделать это из командной строки, так как примеры ориентированы на управление в командном режиме — в основном оно сводится к нажатию клавиши Enter для прерывания пауз, введенных для просмотра промежуточных результатов вычислений.

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

-div( l/sqrt(l+grad|u|^2) * grad(u) ) = 0

при граничном условии u=х*2. Ниже представлен текст файла pdedemo3.m с убранными англоязычными комментариями:

%PDEDEM03 Решение проблемы минимизации параболической поверхности

echo on

clc

g='circleg': % Единичная окружность

b='circleb2': % х^2 - граничное условие

c='l./sqrt(l+ux.*2+uy.*2)':

а=0; 

f=0:

rto1=le-3; % Погрешность вычислений решателем

pause % Пауза в вычислениях

clc

% Генерация поверхности

[p,e,t]=initmesh(g); [p,e.t]-refinemesh(g.p.e.t);

% Решение нелинейной проблемы

u=pdenonlin(b.p.e,t.c.a.f,'tol'.rtol);

pdesurf(p.t.u);

pause % Пауза в вычислениях

echo off

Весьма интересны и поучительны примеры с анимацией: pdedemo2 — появление и распространение волн, pdedemoS — вздутие поверхности (пузырек газа), pdedemo6 — колебания плоскости (гиперболическая проблема) и т. д. К сожалению, они представлены слишком объемными файлами, что не позволяет воспроизвести их в книге учебного характера. Однако вам их несомненно стоит просмотреть самостоятельно.

Что нового мы узнали?

В этом уроке мы научились: