Решение систем уравнений в частных производных matlab

Средства решения систем линейных уравнений (СЛУ)

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

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

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

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

Операции с полиномами

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

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

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

Информация о пакете для решения уравнений с частными производными

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

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

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

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

а 11 X 1 , а 12 ,X 2 . а_1nX n =b₁

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

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

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

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

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

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

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

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

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

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

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

^ — возведение матрицы в степень. Х А р — это X в степени р, если р — скаляр. Если р — целое число, то степень матрицы вычисляется путем умножения X на себя р раз. Если р — целое отрицательное число, то X сначала инвертируется. Для других значений р вычисляются собственные значения и собственные векторы, так что если [V.D]=eig(X), то X*p=V*D. A p/V. Если X — скаляр и Р — матрица, то Х А Р — это скаляр X, возведенный в матричную степень Р. Если X и Р — матрицы, то Х Л Р становится некорректной операцией и система выдает сообщение об ошибке. Возможный вариант решения матричного уравнения АХ=В с применением оператора ^ можно представить как Х=В*А ^ -1.

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

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

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

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

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

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

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. Статистическая формула для стандартной погрешности вычисления коэффициентов имеет следующий вид:

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

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

Решение СЛУ с разреженными матрицами — хотя и не единственная, но несомненно одна из основных сфер применения аппарата разреженных матриц, описанного в уроке 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

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

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

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];

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

Isqr converged at iteration 5 to a solution

with relative residual

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

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

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

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

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

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

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

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

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

BICG converged at iteration 4

to a solution with relative residual

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

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

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

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

BICGSTAB converged at iteration 3.5

to a solution with relative residual

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

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

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

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

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

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

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

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

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

CGS converged at iteration 4 to a solution

with relative residual 4e-014

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

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

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

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

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

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

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

QMR converged at iteration 4 to a solution

with relative residual l.le-014

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

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

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

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

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

Довольно часто возникает задача решения нелинейного уравнения вида f(x) = О или/, (г) =/ 2 (дг). Последнее, однако, можно свести к виду f(x) =f 1 (х) — f 2 (х) = 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 пересекает ось х, а не касается ее. В зависимости от формы задания функции fzero реализуются следующие хорошо известные численные методы поиска нуля функции: деления отрезка пополам, секущей и обратной квадратичной интерполяции. Приведенный ниже пример показывает приближенное вычисление р/2 из решения уравнения cos(x)=0 с представлением косинуса дескриптором:

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

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

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

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

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 33 и являются поэтому исключением: они пока не поддерживают новую нотацию с использованием дескрипторов функций. — Примеч. ред.]

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

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

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

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

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 вычисление минимума выполняется известными методами золотого сечения или параболической интерполяции. Пример:

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

Значительно сложнее задача минимизации функций нескольких переменных 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₁,x 2 ,а) = 100*(x 2 — x 1 ) 2 + (а — x 1 ) 2 .

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

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

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

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

algorithm: ‘Nelder-Mead simplex direct search’ .

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

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

» [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

algorithm: ‘medium-scale: Quasi-Newton line search’

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

Для выполнения аппроксимации Лапласиана в 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. Пример:

На рис. 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, можно строить графики производных заданной функции. Пример этого показан ниже:

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

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

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

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

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

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.

Columns 1 through 7

2.0000 2.0000 2.0000 2.0000 -1.0000-1.50001.0000

Columns 8 through 9

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

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

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 — предельно малая величина.

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

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

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

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

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

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

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

0 1.5000 4.0000 7.5000

10.5000 6.5000 6.5000 10.5000

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

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

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

Функции quad и quadl используют два различных алгоритма квадратуры для вычисления определенного интеграла. Функция quad выполняет интегрирование по методу низкого порядка, используя рекурсивное правило Симпсона [4, 40]. Но она может быть более эффективной при негладких подынтегральных функциях или при низкой требуемой точности вычислений. Алгоритм quad в MATLAB 6 изменен по сравнению с предшествовавшими версиями, точность по умолчанию по сравнению с версиями 5.3х повышена в 1000 раз (с 10- 3 до 10- 6 ). Новая функция 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, не равном нулю, строит график, показывающий ход вычисления интеграла;

quad(@fun,a,b.tol,trace,PI,P2. ) — возвращает численное значение определенного интеграла по хот подынтегральной функции fun, использует дополнительные аргументы P1, P2. которые напрямую передаются в подынтегральную функцию: G=fun(X.Pl,P2. ). Примеры:

dblquad(@fun,inmin,inmax.outmin,outmax) — вычисляет и возвращает значение двойного интеграла для подынтегральной функции fun (Inner, outer), по умолчанию используя квадратурную функцию quad. Inner — внутренняя переменная, изменяющаяся на закрытом интервале от inmin до inmax, a outer — внешняя переменная, изменяющаяся на закрытом интервале от outmin до outmax. Первый аргумент @fun — строка, описывающая подынтегральную функцию. Этс может быть либо дескриптор функции, либо объект inline (в последнем случае символ «@» в ее записи отсутствует). Обычная запись в апострофах тепер недопустима. Эта функция должна быть функцией двух переменных вид. fout=fun(inner.outer). Функция должна брать вектор inner и скаляр outer возвращать вектор fout, который является функцией, вычисленной в outer и каждом значении inner [Функция inime(‘expr’, ‘argl’. ‘argn’) так же создает объект, но без дескриптора, ‘ехрг’ — выражеши Строки ‘argx’ —входные аргументы. При их отсутствии по умолчанию подставляется х. Если вмест ‘arg’ — скаляр, то он означает количество дополнительных переменных Р. Примеры запис; g = inline(exp); g — inline(‘t A 2’); g = inline(‘sin(2*pi*f + theta)’); g = inline(‘sin(2*pi*f + theta)’, : ‘theta’); g — inline(‘x A Pl+P2’, 2). — Примеч. ред.];

dblquad(@fun,inmin.inmax.outmin,oiitmax,tol .trace) — передает в функцию dblquad параметры tol и trace. Смотрите справку по функции quad для получения информации о параметрах to! и trace;

dblquad(@fun,inmin,inmax,outmin,outmax.tol .trace,order) — передает параметры tol и trace для функции quad или quadl в зависимости от значения строки order. Допустимые значения для параметра order — @quad , @quadl или имя любого определенного пользователем квадратурного метода с таким же вызовом и такими же возвращаемыми параметрами, как у функций quad и quad!. (Например, при проверке старых программ можно использовать @quad8 для большей совместимости с прежними версиями MATLAB). По умолчанию (без параметра order) вызывается @quad. поскольку подинтегральные функции могут быть негладкими.

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

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

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

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

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

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

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

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

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

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

14 37 65 91 63 18

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

2.0000 3.0000 5.0000 6.0000

0 0 0.1421 -0.1421-0.2132-0.1066

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

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

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

1.0000 -14.0000 -1.0000-40.0000

1.0000 -58.0000 681.0000 818.0000

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

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

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

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

В=[1.0000 6.4286 1.7145 3.2859]

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

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

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

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

1.0000 -99.0000 626.0000 -626.0000 99.0000-1.0000

-0.0034 -0.0131 -0.0447 -0.0696 -0.1897

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

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

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

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

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

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

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

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

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

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

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

-0.7594

-0.1314

0.3771

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

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

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

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

4.0000 3.0000 1.0000

1.0000 3.0000 7.0000 1.0000

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

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

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

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

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

В этом разделе обобщенное название sol ver (решатель) означает один из возможных численных методов решения ОДУ: ode45, ode23, odel!3, odelSs, 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 -умеренно жестких дифференциальных и дифференциально-алгебраических уравнений.

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

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

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..„) — использует модель SIMU LINK, вызывая соответствующий решатель из нее. Пример:

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

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

NormControl — управление ошибкой в зависимости от нормы вектора решения [on | ]. Установите ‘on’, чтобы norm(e)

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

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

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

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

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

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

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 | ]. Устанавливается на ‘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 | | strong] — установите ‘nоnе’ для уравнений M(t)*y’ = F(t, у). И слабая (‘weak’), и сильная (‘strong’) зависимости означают M(t, у), но ‘weak’ приводит к неявным алгоритмам решения, использующим аппроксимации при решении алгебраических уравнений;

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

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 | ] — указывает, нужно ли использовать формулы обратного дифференцирования (методы Gear) вместо формул численного дифференцирования, используемых в ode 15s по умолчанию;

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

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

Решение систем нелинейных уравнений в Matlab

Доброго времени суток! В этой статье мы поговорим о решении систем нелинейных алгебраических уравнений в Matlab. Вслед за решением нелинейных уравнений, переходим к их системам, рассмотрим несколько методов реализации в Matlab.

Общая информация

Итак, в прошлой статье мы рассмотрели нелинейные уравнения и теперь необходимо решить системы таких уравнений. Система представляет собой набор нелинейных уравнений (их может быть два или более), для которых иногда возможно найти решение, которое будет подходить ко всем уравнениям в системе.

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

Чтобы решить СНАУ, необходимо воспользоваться итеративными методами. Это методы, которые за определенное количество шагов получают решение с определенной точностью. Также очень важно при решении задать достаточно близкое начальное приближение, то есть такой набор переменных, которые близки к решению. Если решается система из 2 уравнений, то приближение находится с помощью построение графика двух функций.

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

Оператор Matlab для решения СНАУ

В среде Matlab существует оператор fsolve, который позволяет решить систему нелинейных уравнений. Сразу рассмотрим задачу, которую, забегая вперед, решим и другими методами для проверки.

Решить систему нелинейных уравнений с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

Нам дана система из 2 нелинейных уравнений и сначала лучше всего построить график. Воспользуемся командой ezplot в Matlab, только не забудем преобразовать уравнения к стандартному виду, где правая часть равна 0:

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

Как видно из графика, есть одно пересечение функций — то есть одно единственное решение данной системы нелинейных уравнений. И, как было сказано, по графику найдем приближение. Возьмем его как (3.0, 1.0). Теперь найдем решение с его помощью:

Создадим функцию m-файлом fun.m и поместим туда следующий код:

Заметьте, что эта функция принимает вектор приближений и возвращает вектор значений функции. То есть, вместо x здесь x(1), а вместо y — x(2). Это необходимо, потому что fsolve работает с векторами, а не с отдельными переменными.

И наконец, допишем функцию fsolve к коду построения графика таким образом:

Таким образом у нас образуется два m-файла. Первый строит график и вызывает функцию fsolve, а второй необходим для расчета самих значений функций. Если вы что-то не поняли, то в конце статьи будут исходники.

И в конце, приведем результаты:

xr (это вектор решений) =
3.3559 1.2069

fr (это значения функций при таких xr, они должны быть близки к 0) =
1.0e-09 *
0.5420 0.6829

ex (параметр сходимости, если он равен 1, то все сошлось) =
1

И, как же без графика с ответом:

Метод простых итераций в Matlab для решения СНАУ

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

Если возможно, строим график.
Из каждого уравнения выражаем неизвестную переменную след. образом: из 1 уравнения выражаем x1, из второго — x2, и т.д.
Выбираем начальное приближение X0, например (3.0 1.0)
Рассчитываем значение x1, x2. xn, которые получили на шаге 2, подставив значения из приближения X0.
Проверяем условие сходимости, (X-X0) должно быть меньше точности
Если 5 пункт не выполнился, то повторяем 4 пункт.

И перейдем к практике, тут станет все понятнее.
Решить систему нелинейных уравнений методом простых итераций с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

График мы уже строили в предыдущем пункте, поэтому переходим к преобразованию. Увидим, что x из первого уравнения выразить сложно, поэтому поменяем местами уравнения, это не повлияет на решение:

x-cos(y) = 3
cos(x-1) + y = 0.5

Далее приведем код в Matlab:

В этой части мы выразили x1 и x2 (у нас это ‘x’ и ‘y’) и задали точность.

В этой части в цикле выполняются пункты 4-6. То есть итеративно меняются значения x и y, пока отличия от предыдущего значения не станет меньше заданной точности.

k = 10
x = 3.3587
y = 1.2088

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

Метод Ньютона в Matlab для решения СНАУ

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

Если возможно, строим график.
Выбираем начальное приближение X0, например (3.0 1.0)
Рассчитываем матрицу Якоби w, это матрица частных производных каждого уравнения, считаем ее определитель для X0.
Находим вектор приращений, который рассчитывается как dx = -w -1 * f(X0)
Находим вектор решения X = X0 + dx
Проверяем условие сходимости, (X-X0) должно быть меньше точности

Далее, решим тот же пример, что и в предыдущих пунктах. Его график мы уже строили и начальное приближение останется таким же.
Решить систему нелинейных уравнений методом Ньютона с точность 10 -2 :
cos(x-1) + y = 0.5
x-cos(y) = 3

Перейдем к коду:

Сначала зададим начальное приближение. Затем необходимо просчитать матрицу Якоби, то есть частные производные по всем переменным. Воспользуемся символьным дифференцированием в Matlab, а именно командой diff с использованием символьных переменных.

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

В этой части кода выполняем первую итерацию, чтобы получить вектор решения и сравнивать его с вектором начального приближения. Отметим, чтобы посчитать значение символьной функции в Matlab, необходимо воспользоваться функцией subs. Эта функция заменяет переменную на числовое значение. Затем функция double рассчитает это числовое значение.

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

В этой части кода выполняется цикл по расчету решения с заданной точностью. Еще раз отметим, что если в первой итерации до цикла были использованы функции diff, double и subs для вычисления производной в Matlab, то в самом цикле матрица якоби была явно задана этими частными производными. Это сделано, чтобы показать возможности среды Matlab.

За 3 итерации достигнут правильный результат. Также важно сказать, что иногда такие методы могут зацикливаться и не закончить расчеты. Чтобы такого не было, мы прописали проверку на количество итераций и запретили выполнение более 100 итераций.

Заключение

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

Применение MATLAB для решения уравнений в частных производных

Для своего варианта функционала выведите дифференциальное уравнение Эйлера-Остроградского и решите его c помощью PDE Toolbox MATLAB.

Примечание: Если на рисунке указано “r = ” – то это часть окружности, если символа “r ” нет – то часть эллипса.

Начало координат − в левом нижнем углу. Граничные условия: на верхней дуге u = 10−3(y−0.15), на остальных сторонах u = 0.

Приложить файл для программы PDE Toolbox MATLAB.
Ниже указан пример готового решения задачи

Применение MATLAB для решения уравнений в частных производных

Рассмотрим графические возможности МАТЛАБ на примере вывода решения дифференциального уравнения в частных производных Эйлера-Остроградского.

Пример 1. Найти экстремаль функционала:

в прямоугольной области xÎ[0,a]; yÎ [0,b], показанной на рис.3.

Рис.3. Область решения примера 1

Граничные условия: на правой стороне x = a:

на остальных сторонах z = 0.

Выведем вначале уравнение Эйлера-Остроградского вида:

или, после сокращения на 2:

. (8)

Граничные условия по переменной y однородные, поэтому будем искать решение в виде ряда Фурье по собственным функциям Yn(y), которые равны:

Ищем решение в виде ряда:

. (9)

Это решение удовлетворяет граничным условиям на нижней и верхней сторонах: при y=0 и y=b.

Для нахождения функций Xn(x) подставим решение (9) в уравнение (8):

Левая часть данного уравнения является разложением в ряд Фурье по Yn(y). Разложим в такой же ряд правую часть этого уравнения. Собственно, она уже разложена: в этом ряду присутствует только первый член (n=1), а коэффициенты при остальных гармониках равны нулю. Известно, что два ряда Фурье тождественны друг другу тогда и только тогда, когда равны все их коэффициенты. Поэтому из полученного уравнения можно получить бесконечную систему дифференциальных уравнений для функций Xn(x):

(10)

Граничные условия для данного уравнения можно получить, раскладывая граничные условия на левой и правой сторонах в ряд Фурье по Yn(y). На левой стороне x=0 имеем z=0, и, следовательно, коэффициенты разложения данной функции нулевые:

. (11)

На правой стороне x=a имеем граничное условие:

Подставим в него решение (9):

Это возможно тогда и только тогда, когда:

(12)

Решаем систему дифференциальных уравнений (10) при граничных условиях (11) и (12). При n>1 имеем:

;

. (13)

Подставляем граничные условия:

Во втором уравнении второй множитель (гиперболический синус) не равен нулю, поэтому C2=0. Следовательно, .

Далее найдем X1(x). Дифференциальное уравнение для него – это 1-е уравнение системы (10). Для решения такого уравнения следует взять сумму общего решения соответствующего однородного уравнения вида (13) и частного решения неоднородного уравнения. В результате получим: