I Error: Missing MATLAB operator

Это не удивительно, т. к. MATLAB ждет от вас 5 + prog (или 5,prog) для вычисления арифметического выражения с переменной prog (или добавле­ния 5 в качестве первого элемента к вектор-строке prog). Следовательно, правильным было бы имя prog5.m (или хотя бы p5rog.m), но только начи­нающееся с буквы.

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

Файл-функции, так же как и файл-программы, создаются в редакторе М-файлов.

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

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

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

function f = myfun(x)

f = exp(-x)*sqrt((x^2 + 1)/(x^4 + 0.1));

 

Слово function в первой строке определяет, что данный файл содержит файл-функцию. Первая строка является заголовком функции, в которой раз­мещаются имя функции и списки входных и выходных аргументов. Входные аргументы -записываются в круглых скобках после имени функции. В нашем примере есть только одни входной аргумент— х. Выходной аргумент f ука­зывается слева от знака равенства в заголовке функции При выборе имени файл-функции следует позаботиться об отсутствии конфликтов с занятыми именами в МАТLAВ. Аналогичным вопрос мы обсуждали выше: как сохра­нить файл-программу в файле с уникальным именем. Тот же самый подход, основанный на обращении к функции exist. вы можете применить для за­дания имени файл-функции.

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

Теперь сохраните файл в рабочем каталоге, Обратите внимание, что выбор пунктов Save или Save as... меню File приводит к появлению диалогового

окна сохранения файла, в поле File name которого уже содержится название myfun. Сохраните файл-функцию в файле с предложенным именем. Теперь

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

» у = myfun(1.3)

У =

0.2600

При создании файл-функцнн myfun мы подавили вывод значения f в командное окно, завершив оператор присваивания точкой с запятой. Если этого не сделать, то оно выведется при обращении у = myfun(1.3). Как правило, лучше избегать вывода в командное окно результатов промежу­точных вычислений внутри файл-функции.

Имя файл-функции не обязательно должно совпадать с именем файла, од­нако обращение к ней происходит по имени файла. Например, если в файле f22.m содержится функция с заголовком g = init (z), то ее следует вызы­вать так:

» f = f22{-0.9)

и вовсе не

» f = init(-0.9)

Функции MATLAB обладают еще одним полезным качеством — возможно­стью получения информации о них при помощи команды help, например, help fplot. Собственные файл-функции так же можно наделить этим свой­ством, используя строки комментариев. Все строки комментариев после за­головка и до тела функции или пустой строки выводятся в командное окно командой help. Изучите содержимое файла fplot.m с файл-функцией fplot, который расположен в подкаталоге \toolbox\matlab\specgraph\ основного каталога MATLAB.

Первая строка комментариев после заголовка функ­ции называется Hl-line и используется при поиске командой lookfor. Эта команда ищет указанное слово в строках Hl-line всех файл-функций в ката­логах, указанных в путях поиска, в том числе и в текущем каталоге.

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

Подфункции. Использование подфункций и вложенных функций основано на выделении части алгоритма в самостоятельную функцию, текст которой содержится в том же файле, что и основная функция. Рассмотрим модельный пример. Предположим, что в файл-функции simple, хранящейся в файле simple.m, часто приходится вычислять некоторое выражение. Конечно, можно про­стым копированием строк добавить соответствующие операторы:

function simple;

x = 1.1;

У = 2.1;

fl = х^3 - 2*y^3 + 3*{x^2 + y^2) - x*y + 9

x = 3.1;

У = 4.2;

f2 = x^3 - 2*у^3 + 3*(x^2 + y^2) - x*y + 9

x = -2.8;

У = 0.7;

f3 = х^3 - 2*у^3 + 3*(х^2 + у^2) - х*у + 9

 

Проще и нагляднее определить вычисляемое выражение в подфункции f с двумя входными и одним выходным аргументом и разместить ее в том же М-файле simple.m.:

function simple;

% Основная функция

fl = f(l.1, 2.1)

f2 = f(3.1, 4.2)

f3 = f( -2.8, 0.7)

function z = f(x, y)

% Подфункция

z = х^3 - 2*y^3 + 3*(x^2 + у^2) - x*y + 9;

 

Первая функция simple является основной функцией в simple.m, именно ее операторы выполняются, если пользователь вызывает simple, например, из командной строки. Каждое обращение к подфункции f в основном функции приводит к переходу к размещенным в подфункции операторам и после­дующему возврату в основную функцию.

Файл-функция может содержать одну или несколько подфункций со своими входными и выходными параметрами, но основная функция может быть только одна. Заголовок новом подфункции одновременно является призна­ком конца предыдущей. Основная функция обменивается информацией к подфункциями только при помощи входных II выходных параметров. Пе­ременные, определенные в подфункциях и в основной функции, являются локальными, они доступны в пределах своей функции.

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

Вложенные функции

Другой разновидностью функций, доступных в одном М-файле, являются вложенные функции. Если подфункция является внешней по отношению к основной функции, то вложенная функция является внутренней. В силу это­го обстоятельства переменные из рабочей среды основной функции доступ­ны и во вложенной функции. Простейшая структура функции main, содер­жащей вложенную fnested:

 

function main;

ALP = 5.3;

BET = 9.1;

fl = fnested(1.1, 2.1)

function z = fnested(x, y}

z = х^3 - 2*y^3 + 3*(x^2 + y^2) - x*y + 9 + ALP'BET;

End

End

 

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

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

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

Приватные функции. Простейший прием сделать некоторые функции видимыми или нет, состоит в размещении М-файлов в рабочих каталогах. Определенная структура ка­талога с пользовательскими файл-функциями позволяет задать некоторые вспомогательные функции, которые используются только файл-функциями, содержащимися в М-файлах данного каталога, а для файл-функцнй из дру­гих каталогов являются недоступными. Для этого следует создать подката­лог с именем private и разместить в нем вспомогательные файл-функции. Рисунок1 поясняет доступ к приватным функциям.

 

Рисунок1

 

задание и порядок выполнения работы

1. Ознакомиться с основными теоретическими сведениями.

2. Выполнить индивидуальное задание на лабораторную работу в соответствии со своим вариантом.

3. Провести тестовый прогон.

4. Сделать выводы по работе.

СОДЕРЖАНИЕ ОТЧЕтА

1. Цель работы

2. Индивидуальное задание на лабораторную работу

3. Краткие теоретические сведения.

4. Блок-схема.

5. Тестовый прогон.

6. Выводы по работе.

7. Листинг.

ВОПРОСЫ ДЛЯ ЗАЩИТЫ ЛАБОРАТОРНОЙ РАБОТЫ

1. Типы файлов в МАТLAB.

2. Действия над файлами.

3. Типы М-файлов и их сравнение.

 

 

лабораторная работа №3

решение обыконвенных дифференциальных
уравнений в пакете МАТLАВ

Цель работы: научиться решать систему обычных дифференциальных уравнений в пакете МАТLАВ.

 

Теоретические сведения

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

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

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

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

Обыкновенными дифференциальными уравнениями называются такие уравнения, которые содержат одну или несколько производных от искомой функции y = y(x). Их можно записать в виде

F(x, y,y’,…,y⁽ⁿ⁾) = 0, (1)

Где x - независимая переменная.

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

F(x, y, y’) = 0, F(x, y, y’, y”) = 0.

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

y’ = f(x, y), (2)

y” = f(x, y, y’).

Такая форма записи называется уравнением, разрешенным относительно старшей производной.

Линейным дифференциальным уравнением называется уравнение, линейное относительно искомой функции и ее производных. Например, y' – x²y = sin x – линейное уравнение первого порядка.

Решением дифференциального уравнения (1) называется всякая n раз дифференцируемая функция y = f(x), которая после ее подстановки в уравнение превращает его в тождество.

Общее решение обыкновенного дифференциального уравнения n-го порядка содержит n

произвольных постоянных C₁, C₂, …, C_n:

y = f(x, C₁, C₂,…,C_n), (3)

где (3) является решением уравнения (1) при любых значениях C₁, C₂, …, C_n, а любое решение уравнения (1) можно представить в виде (3) при некоторых C₁, C₂, …, C_n. Частное решение дифференциального уравнения получается из общего, если произвольным постоянным придать определенные значения.

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

y = f(x, C). (4)

Если постоянная принимает определенное значение C=C₀, то получается частное решение

y = f(x, C₀).

Дадим геометрическую интерпретацию дифференциального уравнения первого порядка (2). Поскольку производная y' характеризует наклон касательной к графику решения

y = y(x) (интегральной кривой) в данной точке, то при y' = k = const из (2) получим

f(x, y) = k – уравнение линии постоянного наклона, называемой изоклиной. Меняя k, получаем, семейство изоклин.

Приведем геометрическую интерпретацию общего решения (4). Это решение описывает бесконечное семейство интегральных кривых с параметром C, а частному решению соответствует одна кривая из этого семейства. При некоторых дополнительных предположениях через каждую точку (x₀,y₀) проходит одна и только одна интегральная кривая. Это утверждение следует из следующей теоремы.

Теорема Коши. Если правая часть f(x, y) уравнения (2) и ее частная производная f'y(x, y) определены и непрерывны в некоторой области G изменения переменных x, y, то для всякой внутренней точки (x₀, y₀) этой области данное уравнение имеет единственное решение, принимающее заданное значение y = y₀ при x = x₀.

Для уравнений высших порядков геометрическая интерпретация более сложна. Через каждую точку в области решения уравнения при n > 1 проходит не одна интегральная кривая. Поэтому если для выделения некоторого частного решения уравнения первого порядка достаточно задать координаты (x₀, y₀) произвольной точки на данной интегральной кривой, то для уравнений высших порядков этого недостаточно. Здесь правило следующее: для выделения частного решения из общего нужно задавать столько дополнительных условий, сколько произвольных постоянных в общем решении, т.е. каков порядок уравнения. Следовательно, для уравнения второго порядка нужно задать два дополнительных условия, благодаря которым можно найти значения двух произвольных постоянных.

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

Если эти условия задаются в одной точке, то такая задача называется задачей Коши. Дополнительные условия в задаче Коши называются начальными условиями, а точка

x = x₀, в которой они задаются, - начальной точкой.

Для уравнения первого порядка дополнительное условие одно, поэтому в этом случае может быть сформулирована только задача Коши: для заданных x₀, y₀ найти такое решение y = y(x) уравнения (2), что y(x₀) = y₀. Таким образом, теорема Коши дает достаточные условия существования и единственности решения задачи Коши.

Если же для уравнения порядка n > 1 дополнительные условия задаются в более чем одной точке, т.е. при разных значениях независимой переменной, то такая задача называется краевой. Сами дополнительные условия называются при этом граничными (или краевыми) условиями. На практике обычно граничные условия задаются в двух точках x = a и x = b, являющихся границами отрезка, на котором рассматривается дифференциальное уравнение.

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

 

dx/dt = x²cost, t > 0, x(0) = 1;

y” = y’/x + x², x > 1, y(1) = 2, y’(1) = 0.

Краевые задачи:

 

y” + 2y’ – y = sinx, 0 <= x <= 1, y(0) = 1, y(1) = 0;

y”’ = x + yy’, 1 <= x <=3, y(1) = 0, y’(1) = 1; y’(3) = 2.

 

Решатели ОДУ в MATLAB

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

 

 

Для решения задачи Коши в MATLAB существует семь солверов:

Таблица №1

ode45	Одношаговые явные методы Рунге-Кутта 4-го и 5-го порядка. Это классический метод, рекомендуемый для начальной пробы решения. Во многих случаях он дает хорошие результаты
ode23	Одношаговые явные методы Рунге-Кутта 2-го и 4-го порядка. При умеренной жесткости системы ОДУ и низких требованиях к точности этот мето;. может дать выигрыш в скорости решения
ode113	Многошаговый метод Адамса-Башворта-Мултона переменного порядка Это адаптивный метод, который может обеспечить высокую точность решения
ode23tb	Неявный метод Рунге-Кутта в начале решения и метод, использующий формулы обратного дифференцирования 2-го порядка в последующем
ode15s	Многошаговый метод переменного порядка (от 1 до 5, по умолчанию 5), использующий формулы численного дифференцирования. Это адаптивный метод, его стоит применять, если решатель ode45 не обеспечивает решения
ode23s	Одношаговый метод, использующий модифицированную формулу Розенброка 2-го порядка. Может обеспечить высокую скорость вычислений при низкой точности решения жесткой системы дифференциальных уравнений
ode23t	Метод трапеций с интерполяцией. Этот метод дает хорошие результаты при решении задач, описывающих колебательные системы с почти гармоническим выходным сигналом

 

Методика их использова­ния одинакова, включая способы задания входных и выходных аргументов:

Таблица №2

[T.Y]=solver(@F,tspan,у0)	Где вместо solver подставляем имя конкретного решателя — интегрирует систему дифференциальных уравнений вида у'=F(t,y) на интервале tspan с начальными условиями у0. @F — дескриптор ODE-функции. Каждая строка в массиве решений Y соответствует значению времени, возвращаемому в векторе-столбце Т
[T.Y]=solver(@F,tspan,y0, options)	Дает решение, подобное описанному выше, но с параметрами, определяемыми значениями аргумента options, созданного функцией odeset. Обычно используемые параметры включают допустимое значение относительной погрешности RelTol (по умолчанию le-З) и вектор допустимых значений абсолютной погрешности AbsTol (все компоненты по умолчанию равны 1е-6)
[T.Y]=solver(@F,tspan,y0, options,pl,p2...)	Дает решение, подобное описанному выше, передавая дополнительные параметрыpi, р2,... в m-файл F всякий раз, когда он вызывается. Используйте options=[], если никакие параметры не задаются
[T.Y.TE.YE.IE]=solver(@F, tspan,y0,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, вызывая соответствующий решатель из нее

 

В достаточно общем случае вызов солвера для решения задачи Коши про­изводится следующим образом:

[Т, У] = solver(odefun,interval,Y0,options)

где odefun — функция для вычисления вектор-функции правой части сис­темы уравнений, interval — массив из двух чисел, задающий промежуток для решения уравнения, Y0 - заданный вектор начальных значений иско­мой вектор-функции, options - структура для управления параметрами и ходом вычислительного процесса. Солвер возвращает массив т с координа­тами узлов сетки, в которых найдено решение, и матрицу решений Y, каж­дый столбец которой является значением компоненты вектор-функции ре­шения в узлах сетки.

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

Схема решения таких задач в MATLAB состоит из следующих этапов.

1 Приведение дифференциального уравнения к системе дифференциальных уравнений первого порядка (если изначально задана система, то в этом нет необходимости).

2 Написание специальной функции для системы уравнений.

3 Вызов подходящего солвера.

4 Визуализация результата.

Разберем решение дифференциальных уравнений на примере уравнения второго порядка:

y'' + 2y' + 10y = sint.

Соответствующие начальные усло­вия имеют вид:

y(0) = 1, y’(0) = 0.

Теперь исходную задачу надо привести к системе дифференциальных уравнений. Для этого вводят столько вспомо­гательных функции, каков порядок уравнения. В данном случае необходи­мы две вспомогательные функции у1 и у2, определяемые формулами:

y₁ = y, y₂ = y’.

Система уравнений с начальными условиями требуемая для работы такова:

y₁’ = y₂; y₁(0) = 1;

y₂’ = -2y₂ - 10 y₁ + sint; y₂(0) = 0.

 

Второй этап состоит в написании функции для системы дифференциальных уравнений. Функция должна иметь два входных аргумента: переменную t, по которой производится дифференцирование, и вектор, размер которого равен числу неизвестных функций системы. Число и порядок аргументов фиксированы, даже если t явно не входит в систему. Выходным аргументом функции является вектор правой части системы. Текст функции oscil для разбираемого примера приведен в листинге1.

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

Function solvdem

% формирование вектора начальных условий

Y0 = [1; 0];

% вызов солвера от функции oscil, начального и конечного

% времени и вектора начальных условий

[T, Y] = ode113(@oscil, [0 15], Y0);

% вывод графика решения исходного дифференциального уравнения

% (маркеры - точки, линия - сплошная)

Plot(T, Y(:, 1), 'r.-')

% вывод графика производной от решения исходного

% дифференциального уравнения (маркеры - точки, линия - пунктир)

Hold on

Plot(T, Y(:, 2), 'k.:')

% вывод пояcнений на rpaфик

title(' {\ity} \prime\prime+2{\ity} \prime+10{\ity} = sin{\itt}')

xlabel('\itt')

ylabel( '{\ity), (\ity) \prime ')

Legend ('y(t) ', 'v(t)', 4)

Grid on

Hold off

% подфункция вычисления правых частей уравнения

function F = oscil(t, y)

F = [y(2); -2*y(2) - 10*y(1) + sin(t)];

 

В результате применения солвера на экран выводятся графики:

В нашем примере решение задачи Коши для системы дифференциальных уравнений, соответствующей исходной задаче для дифференциального уравнения второго порядка, было получено при помощи солвера ode113, который основан на методе Адамса—Бэшфорта—Милтона. Кроме солвера ode113, MATLAB имеет еще ряд солверов для задачи Коши: ode45, ode23, ode15s, ode23s, ode23t_f ode23tb, интерфейс которых не отличается от ode113.

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

Выбор солвера для решения задачи Коши

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

Очень часто солвер ode45 даст вполне хорошие результаты, им стоит вос­пользоваться в первую очередь. Он основан на формулах Рунге—Кутты четвертого и пятого порядка точности. Солвер ode23 также основан на формулах Рунге—Кутты, но уже более низкого порядка точности. Имеет смысл применять ode23в задачах, содержащих небольшую жесткость, когда требуется получить решение с невысокой степенью точности. Если же тре­буется получить решение нежесткий задачи с высокой точностью то наи­лучший результат даст od113, основанный на методе переменного порядка Адамса-Бэшфорта-Милтона. Солвер ode113 оказывается особенно эф­фективным для нежестких систем дифференциальных уравнении, правые части которых вычисляются по сложным формулам. Все солверы пытаются найти решение с относительной точностью 10^-3 и абсолютной - 10^-6, Хо­рошим тестом качества приближенного решения является увеличение точ­ности вычислений (задание точности вычислении и ряда других параметров описано в следующем разделе).

Если все попытки применения ode45, ode23s, ode113 не приводят к успеху, то возможно, что решаемая система является жесткой. Для решения жестких систем подходит солвер ode15s, основанный на многошаговом методе Гира, который допускает изменение порядка. Еслитребуется решить жесткую за­дачу с невысокой точностью, то хороший результат может дать солвер ode23s, реализующий одношаговый метод Розенброка второго порядка.

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

Управление процессом решения

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

 

 

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

Таблица №3

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 Slope	Вектор начального уклона ур0 ур0 = F(t0,y0)/M(t0, y0)
MaxStep	Максимальный шаг, по умолчанию во всех решателях равен одной десятой интервала tspan
BDF(Backward Differentiation Formulas) [on|{off}]	Указывает, нужно ли использовать формулы обратного дифференцирования (методы Gear) вместо формул численного дифференцирования, используемых в ode15sпо умолчанию
MaxOrder	Максимальный порядок ODE15S [1 | 2 | 3 |4 | 5]

 

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

 

 

Таблица №4

Параметры	ode45	ode23	ode113	ode15s	ode23s	ode23t	ode23tb
RelTol, AbsTol, NormControl	+	+	+	+	+	+	+
OutputFcn, OutputSel, Refine, Stats	+	+	+	+	+	+	+
Events	+	+	+	+	+	+	+
MaxStep, InitialStep	+	+	+	+	+	+	+
Jacobian, JPattern, Vectorized	-	-	-	+	+	+	+
Mass	+	+	+	+	+	+	+
MStateDependence	+	+	+	+	-	+	+
MvPattern	-	-	-	+	-	+	+
MassSingular	-	-	-	+	-	+	-
InitialSlope	-	-	-	+	-	+	-
MaxOrder, BDF	-	-	-	+	-	-	-

 

Способ задания параметров солверов ode45, ode23, odell3, odel5s. ode23s, ode23t и ode23tb аналогичен Tому, который вы применяли при нахождении корней функции или локальных минимумов. Значения параметров записываются в управляющую структуру, которая создается функцией odeset. В общем случае обращение к odeset имеет вид:

Options = odeset (…, ‘вид контроля’, значение, …)

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

значение= odeget (options, вид_контроля)

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

Пример выполнения задания:

Задание: Решите дифференциальное уравнение используя формулы обратного дифференцирования. Решение уравнения представьте в графическом виде. Задайте относительную точность решения и выведите ее значение на экран.

Решение:

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

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

Function DIFURI

options=odeset('BDF','on','RelTol',1e-2);