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

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

del 2(11) — возвращает матрицу L дискретной аппроксимации дифференциального оператора Лапласа, примененного к функции U:

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

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

L = del2(U,h) — использует шаг h для установки расстояния между точками в каждом направлении (h — скалярная величина). По умолчанию h=1;

L = de!2(U,hx,hy) — использует hx и hy для точного определения расстояния между точками. Если hx — скалярная величина, то задается расстояние между точками в направлении оси

х,

если hx — вектор, то он должен иметь размер. равный числу столбцов- матрицы U , и точно определять координаты точек по оси

х.

Аналогично если hy — скалярная величина, то задается расстояние между точками в направлении оси

у,

если hy — вектор, то он должен иметь размер. равный числу строк матрицы U, и точно определять координаты точек по оси

у:

L = del2(U,hx.hy,hz....) — если U является многомерным массивом, то расстояния задаются с помощью параметров hx, hy, hz,.... Пример:

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

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

U =

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	1

» subplot(1,2,1)

» surfl(U)

» subplot(1,2,2)

» surfl(V)

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

Рис. 16.1.

Графики функций U и V

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

diff(X) — возвращает конечные разности смежных элементов массива X. Если X — вектор, то diff(X) возвращает вектор разностей соседних элементов [Х(2)-Х(1) Х(3)-Х(2) ... X(n)-X(n-D], у которого количество элементов на единицу меньше, чем у исходного вектора X. Если X — матрица, то diff(X) возвращает матрицу разностей столбцов: [X(2:m, :)-X(l:m-l. :)];

Y = diff(X.n.dim) — возвращает конечные разности для матрицы X по строкам или по столбцам в зависимости от значения параметра dim. Если порядок п равен величине dim или превышает ее, то diff возвращает пустой массив.

Используя функцию 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)

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

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

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

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

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

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

Функции quad и quadl используют два различных алгоритма квадратуры для вычисления определенного интеграла. Функция quad выполняет интегрирование по методу низкого порядка, используя рекурсивное правило Симпсона [4, 40]. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB 6 изменен по сравнению с предшествовавшими версиями, точность по умолчанию по

сравнению с версиями 5.3х повышена в 1000 раз (с 10-

до 10-

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

quad(@fun,a.b) — возвращает численное значение определенного интеграла от заданной функции @fun на отрезке [а Ь]. Используется значительно усовершенствованный в MATLAB 6 адаптивный метод Симпсона;

quad(@fun,a.b.tol) — возвращает численное значение определенного интеграла с заданной относительной погрешностью tol. По умолчанию to1=l.e-6. Можно также использовать вектор, состоящий из двух элементов tol =[rel_tol abs_tol], чтобы точно определить комбинацию относительной и абсолютной погрешности;

quad(@fun,a.b,tol .trace) — возвращает численное значение определенного интеграла и при значении trace, не равном нулю, строит график, показывающий ход вычисления интеграла;

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

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

o=odeget(options, 'name') — извлекает значение свойства, определенного строкой 'name', из структуры параметров options; возвращает пустую матрицу, если значение данного свойства в структуре options не определено. Можно ввести только первые буквы, которые однозначно определяют имя свойства. Пустая матрица [ ] — допустимый аргумент options;

Пример:

» 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

options=odeset( 'namel' .valuel, 'name2' ,value2,...) — создает структуру параметров, в которой указанные свойства по имени ' name...' принимают следующие за ними значения. Вместо 'name...' можно ввести только первые буквы, которые однозначно определяют имя свойства (abs — Abstol, inaxit —maxiter и т. д.);

options=odeset (ol dopts, newopts) — изменяет существующую структуру параметров oldopts путем объединения ее с новой структурой newopts. Все новые параметры, не равные пустой матрице, заменяют соответствующие параметры в структуре oldopts;

options=odeset(ol dopts, 'namel' .valuel,...) — изменяет в существующей структуре параметров соответствующие значения. Пример:

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 ]

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

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

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

Он реализован указанной ниже функцией.

biсд(А. В) — возвращает решение X СЛУ А*Х=В. Матрица коэффициентов А должна быть квадратной размера

пхп,

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

п.

Функция bicg начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером

п,

состоящий из нулей. Итерации производятся или до сходимости к решению, или до появления ошибки, или до достижения максимального числа итераций (по умолчанию равно min(20,n) — либо 20, либо числу уравнений). Сходимость достигается, когда относительный остаток norm(B-A*x)/norm(B) меньше или равен погрешности метода (по умолчанию le-6). Благодаря использованию двунаправленного метода сопряженных градиентов bicg сходится за меньшее число итераций, чем lsqr (в нашем примере быстрее на одну итерацию), но требует квадратную матрицу А, отбрасывая информацию, содержащуюся в дополнительных уравнениях, в то время как 1 sqr работает и с прямоугольной матрицей;

bicgCA.B.tol) — выполняет и возвращает решение с погрешностью (порогом отбора) tol;

bicg(A,b.tol, max it) — выполняет и возвращает решение при заданном максимальном числе итераций maxit;

bicgCA.b.tol,maxit,М) и bicg(A,b.tol .maxit,Ml,М2) — при решении используются матрица предусловий М или М=М1*М2, так что производится решение системы inv(M)*A*x=inv(M)*b относительно х. Если Ml или М2 — пустые матрицы, то они рассматривается как единичные матрицы, что эквивалентно отсутствию входных условий вообще;

bicgCA.B.tol, maxit. Ml. M2.X0) — точно задается начальное приближение Х0. Если Х0 — пустая матрица, то по умолчанию используется вектор, состоящий из нулей;

X = bi eg (А, В, tol, maxit. Ml, M2.X0) — при наличии единственного выходного параметра возвращает решение X. Если метод bicg сходится, выводится соответствующее сообщение. Если метод не сходится после максимального числа итераций или по другой причине, на экран выдается относительный остаток norm(B-A*X)/ norm(B) и номер итерации, на которой метод остановлен;

[X.flag.relres] = bicg(A,X,tol .maxit.Ml,M2.X0) — также возвращает относительную вторую норму вектора остатков relres=nQnr)(B-A*X)/norm(B). Если флаг flag равен 0, то rel res<tol;

[X, flag, rel res, iter] = bicgCA.B.tol,maxit,Ml,M2.XO) — также возвращает номер итерации, на которой был вычислен X. Значение iter всегда удовлетворяет условию 0<iter<maxit;

[X.flag.relres.iter.resvec] = bicgCA.B.tol,maxit,Ml,M2.XO) — также возвращает вектор вторых норм остатков resvec для каждой итерации начиная с res-vec(l)=norm(B-A*X0). Если флаг flag равен 0, то resvec имеет длину iter+1 и resvec(end)<tol*norm(B). Возможны значения flag, равные 0, 1, 2, 3 и 4. Эти значения предоставляют следующие данные о сходимости решения:

flag=0 - решение сходится при заданной точности tol и числе итераций не более заданного maxit;

flag=l - число итераций равно заданному maxit, но сходимость не достигнута;

f l ag=2 - матрица предусловий М плохо обусловлена;

fl ag=3 - процедура решения остановлена, поскольку две последовательные оценки решения оказались одинаковыми;

fl ag=4 - одна из величин в процессе решения вышла за пределы допустимых величин чисел (разрядной сетки компьютера).

Пример:

» 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

[X.flag] = bicg(A,X.tol ,maxit.Ml,M2,X0) — возвращает решение Х и флаг flag, описывающий сходимость метода.

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

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

X =lscov(A,B.V) —возвращаетвекторXрешения СЛУ видаА*Х=В+е, гдее — вектор шумов, который имеет ковариационную матрицу V. Реализуется метод наименьших квадратов в присутствии шумов с известной ковариацией. Прямоугольная матрица А должна быть размера

тхп,

где

т>п.

При решении минимизируется следующее выражение: (AX-b)'*inv(V)*(AX-b). Решение имеет вид X=inv(A'* inv(V)*A)*A'*inv(V)*B. Алгоритм решения, однако, построен так, что операция инвертирования матрицы V не проводится;

[X.dX] = lscov(A,B,V) — возвращает также стандартную погрешность X, помещая ее в переменную dX. Статистическая формула для стандартной погрешности вычисления коэффициентов имеет следующий вид:

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))

X =isqnonneg(A,B) — решение СЛУ АХ=В методом наименьших квадратов с неотрицательными ограничениями. А — действительная прямоугольная матрица, В — действительный вектор. Вектор X содержит неотрицательные элементы X>=Q, где

= 1, 2,...

п.

Критерий: минимизация второй нормы вектора В=АХ;

X = Isqnonneg(A.B.X0) — решение СЛУ с явно заданным неотрицательным начальным значением X для итераций;

[X,W] = Isqnonneg (...) — вместе с решением возвращает вторую норму вектора остатков в квадрате;

[X.W.W1] = Isqnonneg(..) — вместе с решением возвращает вторую норму вектора остатков в квадрате и вектор остатков W1;

[X.W.Wl.exitflag] = Isqnonneg (...) — вместе с решением возвращает вторую норму вектора остатков в квадрате, вектор остатков W1 и флаг exi tflаg, равный 1, если решение сходится после заданного по умолчанию числа итераций, и 0 — в противном случае;

[X.'W.Wl.exitflag,output] = Isqnonneg(...) — дополнительно возвращает число итераций в output;

[X.W.Wl.exitflag,output,lambda] = lsqnonneg(...) — дополнительно возвращает вектор lambda, минимизирующий произведение lambda W1 на последнем шаге итераций решения, lambda (i) близко к нулю, когда X(i) положительно, lambda (i) отрицательно, когда Х(1) равно 0;

[X.W.Wl.exitflag,output,lambda] = lsqnonneg(A.В,ХО,options) — обычно используется, если при предыдущей попытке решения системы exitflag=l или если необходимо изменить заданный по умолчанию порог отбора по X - tol X. равный 10*max(size(A))*norm(A, l)*eps (произведению первой нормы матрицы, большего из измерений матрицы, машинной точности и 10). Также используется такая форма записи, как X=lsqnonneg(A,B,XO,options). Параметры options должны быть предварительно заданы при помощи функции optimset. Функция Isqnonneg использует только поля 'display' и 'tolX' структуры options. Поэтому наиболее часто используемая вместе с Isqnonneg форма записи функции — options= optimset С tol X'.tol value), где tolvalue — новое значение порога отбора по X.

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

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

D =

2.2242

2.6954

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

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

options — аргумент, создаваемый функцией odeset (еще одна функция — odeget или bvpget (только для bvp4c)— позволяет вывести параметры, установленные по умолчанию или с помощью функции odeset /bvpset);

tspan — вектор, определяющий интервал интегрирования [tO tfinal]. Для получения решений в конкретные моменты времени t0, tl,..., tfinal (расположенные в порядке уменьшения или увеличения) нужно использовать tspan = [t0 tl ... tfinal];

у0 — вектор начальных условий;

pi, р2,.„ — произвольные параметры, передаваемые в функцию F;

Т, Y — матрица решений Y, где каждая строка соответствует времени, возвращенном в векторе-столбце Т.

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

[T.Y] = solver(@F,tspan,у0) — где вместо solver подставляем имя конкретного решателя — интегрирует систему дифференциальных уравнений вида

у'=F(t,y)

на интервале tspan с начальными условиями у0. @F — дескриптор ODE-функции. Каждая строка в массиве решений Y соответствует значению времени, возвращаемому в векторе-столбце Т;

[T.Y] = solver(@F, tspan, yO. options) — дает решение, подобное описанному выше, но с параметрами, определяемыми значениями аргумента options, созданного функцией odeset. Обычно используемые параметры включают допустимое значение относительной погрешности RelTol (по умолчанию le-З) и вектор допустимых значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны 1е-6);

[T.Y] = solver(@F,tspan,yO,options,pl,p2...) — дает решение, подобное описанному выше, передавая дополнительные параметры pi, р2,... в m-файл F всякий раз, когда он вызывается. Используйте options=[], если никакие параметры не задаются;

[T.Y.TE.YE.IE] = solver(@F.tspan,yO,options) — в дополнение к описанному решению содержит свойства Events, установленные в структуре options ссылкой на функции событий. Когда эти функции событий от (t, у, равны нулю, производятся действия в зависимости от значения трех векторов value, isterminal, di rection (их величины можно установить в m-файлах функций событий). Для i-й функции событий value(i) —значение функции, isterminal (i) — прекратить интеграцию при достижении функцией нулевого значения, direction^) = 0, если все нули функции событий нужно вычислять (по умолчанию), +1 — только те нули, где функция событий увеличивается, -1 — только те нули, где функция событий уменьшается. Выходной аргумент ТЕ — вектор-столбец времен, в которые происходят события (events), строки YE являются соответствующими решениями, а индексы в векторе IE определяют, какая из i функций событий (event) равна нулю в момент времени, определенный ТЕ. Когда происходит вызов функции без выходных аргументов, по умолчанию вызывается выходная функция odeplot для построения вычисленного решения. В качестве альтернативы можно, например, установить свойство OutputFcn в значение ' odephas2' или 'odephas3' для построения двумерных или трехмерных фазовых плоскостей.

[T.X.Y] = sim(@model,tspan.-y0.options,ut.p1,р2..„) — использует модель SIMULINK, вызывая соответствующий решатель из нее. Пример:

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

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

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

NormControl — управление ошибкой в зависимости от нормы вектора решения [on | {off}]. Установите 'on', чтобы norm(e) <= max(RelTol*norm(y), AbsTol). По умолчанию все решатели используют более жесткое управление по каждой из составляющих вектора решения;

RelTol — относительный порог отбора [положительный скаляр]. По умолчанию 1е-3 (0.1% точность) во всех решателях; оценка ошибки на каждом шаге интеграции e(i) <= max(RelTol*abs(y(i)), AbsTol(i));

AbsTol — абсолютная точность [положительный скаляр или вектор {1е-6}].Скаляр вводится для всех составляющих вектора решения, а вектор указывает на компоненты вектора решения. AbsTol по умолчанию 1е-6 во всех решателях;

Refine - фактор уточнения вывода [положительное целое число] — умножает число точек вывода на этот множитель. По умолчанию всегда равен 1, кроме ODE45, где он 4. Не применяется, если tspan > 2;

OutputFcn — дескриптор функция вывода [function] — имеет значение в том случае, если решатель вызывается без явного указания функции вывода, OutputFcn по умолчанию вызывает функцию odeplot. Эту установку можно поменять именно здесь;

OutputSel — индексы отбора [вектор целых чисел] Установите компоненты, которые поступают в OutputFcn. OutputSel по умолчанию выводит все компоненты;

Stats — установите статистику стоимости вычислений [on {off}];

Jacobian — функция матрицы Якоби [function constant matrix]. Установите это свойство на дескриптор функции FJac (если FJac(t, у) возвращает dF/dy) или

на имя постоянной матрицы dF/dy;

Jpattern — график разреженности матрицы Якоби [имя разреженной матрицы]. Матрица S с S(i,j) = 1, если составляющая i F(t, у) зависит от составляющей j величины у, и 0 в противоположном случае;

Vectorized — векторизованная ODE-функция [on | {off}]. Устанавливается на 'on', если ODE-функция F F(t,[yl y2...]) возвращает вектор [F(t, yl) F(t, y2) ...];

Events — [function] — введите дескрипторы функций событий, содержащих собственно функцию в векторе value, и векторы isterminal и direction (см выше);

Mass — матрица массы [constant matrix function]. Для задач М*у' = f(t, у) установите имя постоянной матрицы. Для задач с переменной М введите дескриптор функции, описывающей матрицу массы;

MstateDependence — зависимость матрицы массы от у [none | {weak} | strong] — установите 'nоnе' для уравнений M(t)*y' = F(t, у). И слабая ('weak'), и сильная ('strong') зависимости означают M(t, у), но 'weak' приводит к неявным алгоритмам решения, использующим аппроксимации при решении алгебраических уравнений;

MassSingular — матрица массы М сингулярная [yes no | {maybe}] (да/нет/может быть);

MvPattern — разреженность (dMv/dy), график разреженности (см функцию spy) — введите имя разреженной матрицы S с S(i,j) = 1 для любого k, где (i, k) элемент матрицы массы M(t, у) зависит от проекции] переменной у, и 0 в противном случае;

Initial Step — предлагаемый начальный размер шага, по умолчанию каждый решатель определяет его автоматически по своему алгоритму;

Initial SI ope — вектор начального уклона ур0 ур0 = F(t0,y0)/M(t0, y0);

MaxStep — максимальный шаг, по умолчанию во всех решателях равен одной десятой интервала tspan;

BDF (Backward Differentiation Formulas) [on | {off}] — указывает, нужно ли использовать формулы обратного дифференцирования (методы Gear) вместо формул численного дифференцирования, используемых в ode 15s по умолчанию;

MaxOrder - Максимальный порядок ODE15S [1 | 2 | 3 4 | {5}].

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

Параметры	Ode45	Ode23	Ode11s	Ode15s	ode23s
RelTol,AbsTol	+	+	+	+	+
OutputFcaOutputSel, Refine, Stats	+	+	+	+	+
Events	+	+	+	+	+
MaxStep, InitlalStep	+	+	+	+	+
Jconstant, Jacobl an,
Jpattern, Vectorized	-	-	-	+	+
Mass	-	-	-	+	+
MassConstant	-	-	-	+	-
MaxOrder, BOF	-		-	+	-

<

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

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

y'

1=

y

2 ;

y'

2=

100*(1-y

1

)^2

*

y

2

-y

1

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

y

1

,(0) = 0;

y

2

(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(@vdp100.[0 30].[2 0]);

» plot(T.Y)

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

Последние команды позволяют с помощью мыши нанести на графики решений y

1

= y(1) и у

2

= y(2) помечающие их надписи.

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

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

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

, а

..., а

Здесь а

, а,

,...,

пп

—

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

, х

,..., х

—

неизвестные, образующие вектор X, и b

, b

,..., b

—

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

Для реализации различных алгоритмов решения СЛУ и связанных с ними матричных операций применяются следующие операторы:

+,-,*,/, \, *, ' .

Как отмечалось ранее, MATLAB имеет два различных типа арифметических операций — поэлементные и для массивов (векторов и матриц) в целом. Матричные арифметические операции определяются

правилами линейной алгебры.

Арифметические операции сложения и вычитания над массивами выполняются поэлементно. Знак точки «.» отличает операции над элементами массивов от

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

* — матричное умножение;

С = А*В — линейное алгебраическое произведение матриц А и В:

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

/ —

правое

деление. Выражение Х=В/ А дает решение ряда систем линейных уравнений АХ=В, где А — матрица размера

тхп

и В — матрица размера

nxk;

\ —

левое

деление. Выражение Х=В\А дает решение ряда систем линейных уравнений ХА=В, где А — матрица размера

тхп

и В — матрица размера

nxk.

Если А — квадратная матрица, то А\В — примерно то же самое, что и inv(A)*B, в остальных случаях возможны варианты, отмеченные ниже.

Если А — матрица размера

пхп,

а В — вектор-столбец с

п

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

Если А — матрица размера

тхп

и

тхп,

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

Ранг

k

матрицы А находится на основе QR-разложения (урок 11) с выбором ведущего элемента. Полученное решение X будет иметь не больше чем

k

ненулевых компонентов на столбец. Если

k<n,

то решение, как правило, не будет совпадать с pinv(A)*B, которое имеет наименьшую норму невязок ||Х||.

^

—

возведение матрицы в степень. Х

А

р — это X в степени р, если р — скаляр. Если р — целое число, то степень матрицы вычисляется путем умножения X на себя р раз. Если р — целое отрицательное число, то X сначала инвертируется. Для других значений р вычисляются собственные значения и собственные векторы, так что если [V.D]=eig(X), то X*p=V*D.

A

p/V. Если X — скаляр и Р — матрица, то Х

А

Р — это скаляр X, возведенный в матричную степень Р. Если X и Р — матрицы, то Х

Л

Р становится некорректной операцией и система выдает сообщение об ошибке. Возможный вариант решения матричного уравнения АХ=В с применением оператора

^

можно представить как Х=В*А

^

-1.

' — транспонирование матрицы, то есть замена строк столбцами и наоборот. Например, А' — транспонированная матрица А. Для комплексных матриц транспонирование дополняется комплексным сопряжением. Транспонирование при решении СЛУ полезно, если в матрице А переставлены местами столбцы и строки.

При записи СЛУ в матричной форме необходимо следить за правильностью записи матрицы А и вектора В. Пример (в виде 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.

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

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

cgs(A.B) — возвращает решение X СЛУ А*Х=В. А — квадратная матрица. Функция cgs начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размера

п,

состоящий из нулей. Итерации производятся либо до сходимости метода, либо до появления ошибки, либо до достижения максимального числа итераций. Сходимость метода достигается, когда относительный остаток norm(B-A*X)/norm(B) меньше или равен погрешности метода (по умолчанию le-6). Функция cgs(...) имеет и ряд других форм записи, аналогичных описанным для функции bieg(...). Пример:

» 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

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

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

qmr (А, В) — возвращает решение X СЛУ А*Х=Ь. Матрица А должна быть квадратной. Функция qmr начинает итерации от начальной оценки, представляющей собой вектор длиной

п,

состоящий из нулей. Итерации производятся либо до сходимости метода, либо до появления ошибки, либо до достижения максимального числа итераций. Максимальное число итераций — минимум из п и 20. Функция qmr(...) имеет и ряд других форм записи, аналогичных описанным для функции bi eg (...). Пример:

» qmr(A.B)

QMR converged at iteration 4 to a solution

with relative residual 1.1e-014

ans =

1.0000

2.0000

3 .0000

4.0000

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

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

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

» 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

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

Итерационный метод сопряженных градиентов реализован функцией 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-алгоритм итерационного метода сопряженных градиентов и также не требует, чтобы ее входной аргумент — квадратная симметрическая матрица — была положительно определенной. Эта функция тоже не может решить наш пример, так как наша матрица А — квадратная, но не симметрическая.

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

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

трапеций

и методом

трапеций с накоплением.

trapz(Y) — возвращает определенный интеграл, используя интегрирование методом трапеций с единичным шагом между отсчетами. Если Y — вектор, то trapz(Y) возвращает интеграл элементов вектора Y, если Y — матрица, то trapz(Y) возвращает вектор-строку, содержащую интегралы каждого столбца этой матрицы;

trapz(X.Y) — возвращает интеграл от функции Y по переменной X, используя метод трапеций (пределы интегрирования в этом случае задаются начальным и конечным элементами вектора X);

trapz(...,dim) — возвращает интеграл по строкам или по столбцам для входной матрицы в зависимости от значения переменной dim. Примеры:

»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) — возвращает численное значение определенного интеграла для функции, заданной ординатами в векторе или матрице Y с шагом интегрирования, равным единице (интегрирование методом трапеций с накоплением). В случае когда шаг отличен от единицы, но постоянен, вычисленный интеграл достаточно умножить на величину шага. Для векторов эта функция возвращает вектор, содержащий результат интегрирования с накоплением элементов вектора Y. Для матриц — возвращает матрицу того же размера, что и Y, содержащую результаты интегрирования с накоплением для каждого столбца матрицы Y;

cumtrapz(X, Y) — выполняет интегрирование с накоплением от Y по переменной X, используя метод трапеций. X и Y должны быть векторами одной и той же длины или X должен быть вектором-столбцом, a Y — матрицей;

cumtrapz(...,dim) — выполняет интегрирование с накоплением элементов по размерности, точно определенной скаляром dim. Длина вектора X должна быть равна size(Y.dim). Примеры:

» 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

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

Значительно сложнее задача минимизации функций нескольких переменных

f(х

,...).

При этом значения переменных представляются вектором

х,

причем начальные значения задаются вектором

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

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

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

fminsearch(@fun,xO) — возвращает вектор х, который является локальным минимумом функции fun(x) вблизи хО.хО может быть скаляром, вектором (отрезком) при минимизации функции одной переменной или матрицей (для функции нескольких переменных);

fminsearch(@fun,xO,options) — аналогична описанной выше функции, но использует вектор параметров options точно так же, как функция fminbnd;

fminsearch(@fun,xO,options.P1.P2,...) — сходна с описанной выше функцией, но передает в минимизируемую функцию нескольких переменных fun(x.P1,P2....) ее дополнительные аргументы Р1. Р2,.... Если требуется использовать параметры вычислений по умолчанию, то вместо options перед Р1, Р2 необходимо ввести [ ].;

[x.fval] = fminsearchC...) — дополнительно возвращает значение целевой функции fval в точке минимума;

[x.fval .exitflag] = fminsearchC...) —дополнительно возвращает параметр exitflag, положительный, если процесс итераций сходится с использованием options. tol X, отрицательный, если итерационный процесс не сходится к полученному решению х, и 0, если превышено максимальное число итераций options. maxi ten;

[х. fval .exitflag.output] - fminsearch(...) возвращает структуру (запись) output,

output.algorithm — использованный алгоритм;

output. funcCount — число оценок целевой функции;

output.Iterations — число проведенных итераций.

Классическим примером применения функции fminsearch является поиск минимума тестовой функции Розенброка, точка минимума которой находится в «овраге» с «плоским дном»: rb(x

1

,x

2

,а) = 100*(x

2

- x

1

)

2

+ (а - x

1

)

2

.

Минимальное значение этой функции равно нулю и достигается в точке [

а

а

2

]. В качестве примера уточним значения

x

1

и

х

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',1.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=

1

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

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

Еще одна важная задача численных методов — поиск

минимума

функции

f(x)

в некотором интервале изменения

х —

от

до

. Если нужно найти

максимум

такой функции, то достаточно поставить знак «минус» перед функцией. Для решения этой задачи используется следующая функция:

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

fminbnd(@fun,xl,x2) — возвращает значение х, которое является локальным минимумом функции fun(x) на интервале xl<x<x2;

fminbnd(@fun,xl,x2.options) — сходна с описанной выше формой функции, но использует параметры to!X, maxfuneval, maxiter, display из вектора options, предварительно установленные при помощи команды optimset (смотрите описание lsqnonneg);

fminbnd(@fun,xl.x2,options.P1.P2...) — сходна с описанной выше, но передает в целевую функцию дополнительные аргументы: Р1, Р2..... Если требуется использовать параметры вычислений по умолчанию, то вместо options перед P1, Р2 необходимо ввести [ ] (пустой массив);

[x.fval] = fminbnd(...) — дополнительно возвращает значение целевой функции fval в точке минимума;

[x.fval .exitflag] = fminbndL.) —дополнительно возвращает параметр exitflag, равный 1, если функция сошлась с использованием options.tolX, и 0, если достигнуто максимальное число итераций options.maxiter.

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

которая будет минимизирована; options — вектор параметров вычислений. В зависимости от формы задания функции fminbnd вычисление минимума выполняется известными методами золотого сечения или параболической интерполяции. Пример:

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

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

х =

3.1416

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

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

y'= F(t,

у), My' = F(t, у)

или

M(t)y'

Y(t, у),

где t — независимая переменная (скаляр), которая обычно представляет время;

— вектор зависимых переменных; F — функция от

у,

возвращающая вектор-столбец такой же длины как и

у;

М и М(£) — матрицы, которые не должны быть вырожденными. М может быть и константой.

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

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

Подготовим 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

Тогда решение системы с помощью решателя 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, можно наблюдать результат решения в виде графиков. В решении вы найдете структуру узлов начальной сетки решения, которая поясняется ниже.

solinit — это структура узлов начальной сетки решения (в любой шкале), но такая, что solinit.x=a, solinit.x - b. И функция у, и функция у' должны быть непрерывны на участке [а b]. Дотадка для начальной итерации so1imt=bvpi-mt(x,yinit.params) в bvp4c отличается тем, что ваше начальное представление о функции у yinit можно вводить не только в виде вектора, но и как символьную функцию.

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

Пакет 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 — колебания плоскости (гиперболическая проблема) и т. д. К сожалению, они представлены слишком объемными файлами, что не позволяет воспроизвести их в книге учебного характера. Однако вам их несомненно стоит просмотреть самостоятельно.

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

Полиномы

(у нас их принято называть также

степенными многочленами)

— широко известный объект математических вычислений и обработки данных. Обычно полином записывается в виде

р(х)

х^n + x

п-1

x^n -1+ ...

x^2 + а

^х + а

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

р(х) = a

x^n + а

x^n-1 + ... + а

х + а

п+1

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

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

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

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

[r,p,k] = res I due (b, a) — возвращает вычеты, полюса и многочлен целой части отношения двух полиномов b(s) и a(s) в виде

[b.a] = residue(r.p.k) — выполняет обратную свертку суммы простых дробей (см. более подробное описание в справочной системе) в пару полиномов с коэффициентами в векторах b и а.

Пример:

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

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

Решатели ОДУ

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

решателями

ОДУ.

Примечание

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

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

ode45 — одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка. Это классический метод, рекомендуемый для начальной пробы решения. Во многих случаях он дает хорошие результаты;

ode23 — одношаговые явные методы Рунге-Кутта 2-го и 4-го порядка. При умеренной жесткости системы ОДУ и низких требованиях к точности этот мето;. может дать выигрыш в скорости решения;

ode113 — многошаговый метод Адамса-Башворта-Мултона переменного порядка Это адаптивный метод, который может обеспечить высокую точность решения

ode23tb — неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем

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

ode15s — многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения;

ode23s — одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности решения жесткой системы дифференциальных уравнений;

ode23t — метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих колебательные системы с почти гармоническим выходным сигналом;

bvp4c служит для проблемы граничных значений систем дифференциальных уравнений вида y

'==

f(t,y), F(y(a), y(b), p)

=

0 (краевая задача);

pdepe нужен для решения систем параболических и эллиптических дифференциальных уравнений в частных производных, введен в ядро системы для поддержки новых графических функций Open GL, пакет расширения Partial Differential Equations Toolbox содержит более мощные средства.

Все решатели могут решать системы уравнений явного вида

у'

= 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 -умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

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

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

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

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

с граничными условиями y(t

end

, p)

b, где t

end

, t

— начальные и конечные точки интервалов. Параметр

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

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

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

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

[Х.е] = polyeig(AO,Al,...Ap) — решает задачу собственных значений для матричного полинома степени р вида:

где степень полинома р — целое неотрицательное число, а А

,..., А

—

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

п.

Выходная матрица 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

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

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

Использование всех этих функций, кроме Isqr, как правило, ограничено системами уравнений, в которых А — нормальная квадратная матрица, т. е. А* А -АА*, где А*— сопряженная (эрмитова) матрица А. — Примеч. ред.

]

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

lsqr(A, В)—возвращает точное решение X СЛУ А*Х=В, если матрица последовательная, в противном случае — возвращает решение, полученное итерационным методом наименьших квадратов. Матрица коэффициентов А должна быть прямоугольной размера тхя, а вектор-столбец правых частей уравнений В должен иметь размер т. Условие m>=n может быть и необязательным. Функция 1 sqr начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером

п,

состоящий из нулей. Итерации производятся или до сходимости к решению, или до появления ошибки, или до достижения максимального числа итераций (по умолчанию равного min(20, m, n) — либо 20, либо числу уравнений, либо числу неизвестных). Сходимость достигается, когда отношение вторых норм векторов norm(B-Ax)/norm(B) меньше или равно погрешности метода tol (по умолчанию 1е-б);

lsqr(A.B,tol) — возвращает решение с заданной погрешностью (порогом отбора) tol;

lsqr(A,b.tol .maxlt) — возвращает решение при заданном максимальном числе итераций maxit вместо, возможно, чересчур малого числа, заданного по умолчанию;

lsqr(A,b.tol .maxit,M) и lsqr(A,b,tol .maxit.Ml.M2) — при решении используются матрица предусловий М или М=М1*М2, так что производится решение системы inv(M)*A*x=inv(M)*b относительно х. Если Ml или М2 — пустые матрицы, то они рассматривается как единичные матрицы, что эквивалентно отсутствию входных условий вообще;

lsqr(A.B,tol .maxit.Ml.M2.X0) — точно задается начальное приближение Х0. Если Х0 — пустая матрица, то по умолчанию используется вектор, состоящий из нулей;

X = lsqr(A,B.tol .maxit,Ml.M2.X0) — при наличии единственного выходного параметра возвращает решение X. Если метод 1 sqr сходится, выводится соответствующее сообщение. Если метод не сходится после максимального числа итераций или по другой причине, на экран выдается относительный остаток попп(В-А*Х)/ norm(B) и номер итерации, на которой метод остановлен;

[X.flag] = lsqr(A.X.tol.maxit.Ml.M2.X0) — возвращает решение X и флаг flag. описывающий сходимость метода;

[X.flag.relres] = lsqr(A,X,tol.maxit,Ml.M2.X0) — также возвращает относительную вторую норму вектора остатков rel res=norm(B-A*X)/norm(B). Если флаг flag равен 0, то relres<tol;

[X.flag.relres.iter] = bicg(A,B.tol,maxit,Ml,M2.X0) — также возвращает номер итерации, на которой был вычислен X. Значение iter всегда удовлетворяет условию 0<iter<maxit;

[X.flag.relres,iter,resvec]= lsqr(A.B.tol.maxit.Ml.M2.X0) — также возвращает вектор вторых норм остатков resvec для каждой итерации начиная с res-vec(l)=norm(B=A*X0). Если флаг flag равен 0, то resvec имеет длину iter+1 и resvec(end)<tol*norm(B). Возможны значения flag, равные 0, 1, 2, 3 и 4. Значения flag предоставляют следующие данные о сходимости решения:

flag=0 - решение сходится при заданной точности tol и числе итераций не более заданного maxit;

flag=l - число итераций равно заданному maxit, но сходимость не достигнута;

flag=2 - матрица предусловий М плохо обусловлена;

flag=3 - процедура решения остановлена, поскольку две последовательные оценки решения оказались одинаковыми;

fl ag=4 - одна из величин в процессе решения вышла за пределы допустимых величин чисел (разрядной сетки компьютера).

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

Пример:

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

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

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

13

.

Пример:

» lsqr(A,B.1e-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

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

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

w = conv(u.v) — возвращает свертку векторов и и v. Алгебраически свертка — то же самое, что и произведение полиномов, чьи коэффициенты — элементы векторов и и v. Если длина вектора и равна

т,

а длина вектора v —

п,

то вектор w имеет длину

т+п-1,

а его

k-й

элемент вычисляется по следующей формуле

Пример:

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

r =

14 37 65 91 63 18

[q,r] = deconv(v.u) —возвращает результат деления полинома v на полином и. Вектор q представляет собой частное от деления, а г — остаток от деления, так что выполняется соотношение v=conv(u,q)+r.

Пример:

» 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

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

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

bicgstab:

bicgstab(A.B) — возвращает решение X СЛУ А*Х=В. А — квадратная матрица. Функция bi cgstab начинает итерации от начальной оценки, по умолчанию представляющей собой вектор размером

п,

состоящий из нулей. Итерации производятся либо до сходимости метода, либо до появления ошибки, либо до достижения максимального числа итераций. Сходимость метода достигается, когда относительный остаток norm(B-A*X)/norm(B) меньше или равен погрешности метода (по умолчанию le-б). Функция bicgstab(...) имеет и ряд других форм записи, аналогичных описанным для функции bi eg (...). Пример:

» 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

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

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

FX = gradient(F) — возвращает градиент функции одной переменной, заданной вектором ее значений F. FX соответствует конечным разностям в направлении

х,

[FX.FY] = gradient(F) — возвращает градиент функции F(X,Y) двух переменных, заданной матрицей F, в виде массивов FX и FY. Массив FX соответствует конечным разностям в направлении

(столбцов). Массив FY соответствует конечным разностям в направлении

(строк);

[FX.FY.FZ,...] = gradient(F) — возвращает ряд компонентов градиента функции нескольких переменных, заданной в виде многомерного массива F;

[...] = gradient(F.h) — использует шаг h для установки расстояния между точками в каждом направлении (h — скалярная величина). По умолчанию h=l;

[...] = gradient(F.hi,h2,...) — если F является многомерным массивом, то расстояния задаются с помощью параметров h1, h2, h3,....

Пример:

» 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 часто используется для построения графиков полей градиентов.

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

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

Примечание

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

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

» feval(fe.1.0)

ans =

2.7183

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

Довольно часто возникает задача решения нелинейного уравнения вида

f(x) =

О или/, (г) =/

(дг). Последнее, однако, можно свести к виду

f(x)

(х) -

(х) =

0. Таким образом, данная задача сводится к нахождению значений аргумента

функции

f(x)

одной переменной, при котором значение функции равно нулю. Соответствующая функция MATLAB, решающая данную задачу, приведена ниже:

fzero(@fun,x) — возвращает уточненное значение х, при котором достигается нуль функции fun, представленной в символьном виде, при начальном значении аргумента х. Возвращенное значение близко к точке, где функция меняет знак, или равно NaN, если такая точка не найдена;

fzero(@fun,[xl x2]) — возвращает значение х, при котором fun(x)=0 с заданием интервала поиска с помощью вектора x=[xl х2], такого, что знак fun(x(D) отличается от знака fun(x(2)). Если это не так, выдается сообщение об ошибке. Вызов функции fzero с интервалом гарантирует, что fzero возвратит значение, близкое к точке, где fun изменяет знак;

fzero(@fun,x.tol) — возвращает результат с заданной погрешностью tol;

fzero(@fun,x.tol .trace) — выдает на экран информацию о каждой итерации;

fzero(@fun,х.tol .trace,Р1.Р2,...) — предусматривает дополнительные аргументы, передаваемые в функцию fun(x.Pl,P2,...). При задании пустой матрицы для tol или trace используются значения по умолчанию. Пример:

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

Для функции fzero ноль рассматривается как точка, где график функции fun

пересекает

ось

х,

а не

касается

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

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

x =

1.5708

В более сложных случаях настоятельно рекомендуется строить график функции

f(x)

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

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

function f=funl(x)

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

» х=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

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

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

poly(A) — для квадратной матрицы А размера

пхп

возвращает вектор-строку размером n+1, элементы которой являются коэффициентами характеристического полинома det(A-sI), где I — единичная матрица, as — оператор Лапласа. Коэффициенты упорядочены по убыванию степеней. Если вектор состоит из

п+1

компонентов, то ему соответствует полином вида c

s^n+...+c

s+c

n+1

;

poly (г) — для вектора г возвращает вектор-строку р с элементами, представляющими собой коэффициенты полинома, корнями которого являются элементы вектора г. Функция roots(p) является обратной, ее результаты, умноженные на целое число, дают poly (r ).

А =

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

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

roots (с) — возвращает вектор-столбец, чьи элементы являются корнями полинома с.

Вектор-строка с содержит коэффициенты полинома, упорядоченные по убыванию степеней. Если с имеет 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

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

polyval (p,x) — возвращает значения полинома р, вычисленные в точках, заданных в массиве х. Полином р — вектор, элементы которого являются коэффициентами полинома в порядке уменьшения степеней, х может быть матрицей или вектором. В любом случае функция polyval вычисляет значения полинома р для каждого элемента х;

[у.delta] = polyval (p. x.S) или [у,delta] = polyval (p.x.S.mu)— использует структуру S, возвращенную функцией polyfit, и данные о среднем значении (mu(l)) и стандартном отклонении (mu(2)) генеральной совокупности для оценки пр-грешности аппроксимации (y+delta).

Пример:

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

d =

35 675

polyvalm(p.X) — вычисляет значения полинома для матрицы. Это эквивалентно подстановке матрицы X в полином р. Полином р — вектор, чьи элементы являются коэффициентами полинома в порядке уменьшения степеней, а X — квадратная матрица.

Пример:


» 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

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

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

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

polyder(p) — возвращает производную полинома р;

polyder(a.b) — возвращает производную от произведения полиномов а и b;

[q,d] = polyder(b.a) — возвращает числитель q и знаменатель d производной от отношения полиномов b/а.

Примеры:

» 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