авторефераты диссертаций БЕСПЛАТНАЯ  БИБЛИОТЕКА

АВТОРЕФЕРАТЫ КАНДИДАТСКИХ, ДОКТОРСКИХ ДИССЕРТАЦИЙ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ

Астрологический Прогноз на год: карьера, финансы, личная жизнь


Математика в высшем образовании

2005 №3

ИННОВАЦИОННЫЕ И ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

И

КОМПЬЮТЕРНЫЕ ПРОДУКТЫ В ПРЕПОДАВАНИИ МАТЕМАТИКИ

УДК 517.91 + 519.67

ОБЫКНОВЕННЫЕ ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ:

MATLAB VERSUS MATHCAD

Ю. Л. Кетков1, А. И. Кузнецов2

1

Нижегородский государственный университет им. Н. И. Лобачевского, Россия, 603950, г. Нижний Новгород, пр. Гагарина, 23, корп. 2;

e-mail: ket@city.ru 2 НИИ Прикладной математики и кибернетики, Россия, 603005, г. Нижний Новгород, ул. Ульянова, 10;

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

Ключевые слова: дифференциальные уравнения, численные методы, жст е кие системы дифференциальных уравнений, MatLab, MathCAD.

1. ВВЕДЕНИЕ Численные методы решения систем линейных и нелинейных обыкновен ных дифференциальных уравнений (ОДУ) один из наиболее важных раз делов “Прикладной математики”. С ними приходится сталкиваться как в учебном процессе при изучении алгоритмов решения таких задач, так и в научно-инженерной деятельности при исследовании сложных динамических систем. Поэтому представляют интерес оценки тех или иных программных и инструментальных средств, доступных на современной вычислительной тех нике. Численные методы не единственный способ исследования динами ческих систем. В инженерной практике широко используются методы ма тематического моделирования, заложенные в 50–60 гг. прошлого столетия и базирующиеся на использовании идей аналоговой вычислительной техники.

Сегодня они нашли свое воплощение в таких инструментальных средствах, как Simulink (составная часть MatLab [1]), GPSS World [2] и др.

Мы сознательно ограничили область исследования только сравнением чи сленных методов, заложенных в соответствующих разделах наиболее мощных инструментов вычислительной математики MatLab и MathCAD. За послед ние годы эти пакеты вс более активно внедряются в вузовскую практику е и, кроме соображений экономического характера, полезно знать, с какими возможностями и ограничениями могут столкнуться преподаватели соответ ствующих дисциплин.

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

y1 = f1 (t, y1, y2,...), y2 = f2 (t, y1, y2,...),....................

yn = fn (t, y1, y2,...). (1) Дифференциальные уравнения высокого порядка, разрешимые относи тельно старшей производной, приводятся к форме (1) путем тривиальных подстановок. Продемонстрируем это на примере линейной модели гармони ческого осциллятора:

2 y + y + y = 0. (2) После подстановок y1 = y и y2 = y уравнение (2) сводится к системе второго порядка:

y1 = y2, y2 = 2 (y2 + y1 ).

Задача Коши. Задача Коши для ОДУ заключается в построении траек тории, удовлетворяющей системе (1) и проходящей через заданную началь ную точку y1 (t0 ), y2 (t0 ),..., yn (t0 ).

Краевая задача. В отличие от задачи Коши, где все n условий заданы в начальной точке, краевая задача должна отыскать такую траекторию, ко торая удовлетворяет n k условиям в начальной точке t0 и k условиям на конце заданного интервала интегрирования.

Неявные ОДУ. В самом общем виде дифференциальное уравнение n-го порядка может быть задано с помощью неявной функции:

F t, y (n), y (n1),..., y, y = 0. (3) Далеко не всегда его удается свести к системе вида (1). И тогда прихо дится изыскивать специальные методы решения уравнений, не разрешенных в явном виде относительно производных. Примером такой задачи является анализ заряда q(t) на обкладках конденсатора в колебательном контуре с индуктивностью, снабженной железным сердечником [3]:

2 + B · (q ) + q 2 = E.

A ln 1 + (q ) (4) Дифференциальные уравнения с запаздывающим аргументом.

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

y (t) = f t, y(t), y(t 1 ), y(t 2 ),..., y(t k ). (5) Математика в высшем образовании 2005 № Одношаговые методы численного решения ОДУ. К этой группе относятся многочисленные схемы типа Рунге – Кутта, базирующиеся на сле дующей схеме решения дифференциального уравнения 1-го порядка [4, 5] и легко обобщаемые на системы вида (1):

y = f (t, y), k1 = h · f (t, y), k2 = h · f (t + 2 h, y + 2,1 · k1 ), k3 = h · f (t + 3 h, y + 3,1 · k1 + 3,2 · k2 ),.........................................

km = h · f (t + m h, y + m,1 · k1 + m,2 · k2 +... + m,m1 · km1 ), y(t + h) = y(t) + 1 · k1 + 2 · k2 +... + m · km. (6) Постоянные коэффициенты i, i,j, i подбираются таким образом, чтобы для произвольной функции f (t, y) формула (6) давала наилучшее приближе ние к отрезку ряда Тейлора:

y(t + h) y(t) + hy (t)/1! + h2y (t)/2! +... + hm y (m) (t)/m!. (7) Количество коэффициентов kj определяет погрешность интегрирования.

Одной из наиболее часто применяемых на практике схем являются формулы Рунге – Кутта 4-го порядка:

k1 = h f (tk, yk ), k2 = h f (tk + 0,5 h, yk + 0,5 k1 ), k3 = h f (tk + 0,5 h, yk + 0,5 k2 ), k4 = h f (tk + h, yk + k3 ), yk+1 = yk + (k1 + 2k2 + 2k3 + k4 )/6.



Довольно часто, с целью контроля за точностью вычислений, формулы Рунге – Кутта применяют одновременно с шагами h и h/2 (при таком подборе шагов удается кое-что сэкономить). И тогда схему интегрирования принято называть формулой интегрирования порядка n/(n + 1).

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

yk+1 = yk + h(0 f k + 1 fk1 + 2 fk2 +...). (8) На начальной стадии (этап разгона) значения правых частей в точках t0, t0 + h, t0 + 2h,... вычисляются по одной из одношаговых схем типа Рун ге – Кутта, зато скорость вычислений на последующих шагах интегрирования увеличивается примерно в 4 раза. Восьми- и девятиточечные схемы Адамса хорошо зарекомендовали себя при расчете траекторий космических аппара тов.

Ю. Л. Кетков, А. И. Кузнецов 2. СИСТЕМЫ ОДУ С ГЛАДКИМ РЕШЕНИЕМ В этом разделе будут рассматриваться средства решения задачи Коши для систем ОДУ вида (1), описывающих “хорошие” динамические системы.

Под “хорошими” имеются в виду системы с достаточно гладкими решения ми амплитуды колебаний неизвестных функций yi (t) и их основные гармо ники не очень сильно отличаются друг от друга. Примером системы такого типа является хорошо изученная линейная модель гармонического осцилля тора, описываемого уравнением (2).

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

u:=name fun(y0, t0, t1, M, D).

Здесь:

name fun имя функции интегрирования (см. табл. 1);

y0 вектор-столбец начальных значений в момент t0;

t1 конец интервала интегрирования (допускается t1 t0);

M количество точек внутри интервала интегрирования, в которых нуж но получить решение;

D векторная функция двух аргументов (независимой переменной t и вектора y), описывающая правые части системы ОДУ.

Функция name fun возвращает матрицу с результатами интегрирования:

е первый столбец u 0 содержит массив значений независимой переменной е t, остальные столбцы u 1, u 2,... содержат компоненты искомой вектор функции y.

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

Рис.1. График зависимости u1 от u На рис. 1 представлен график построенной функции y0 (t). Обратите внимание на использование символа присвоения (:=) в каждом из опера торов приведенного примера (он включается в текст сочетанием клавиш Shift + :).

Математика в высшем образовании 2005 № Функции для решения ОДУ в таком формате получили название функций командной строки. Они были сохранены в MathCAD’е для преемственности со старыми версиями и для обращения из программ, где нельзя использовать многострочную форму записи ОДУ. Еще одним специальным применением функций командной строки является задача о получении значения неизвест ной функции только в конце интервала интегрирования.

Таблица № Имя функции Пояснение 1 Rkxed Метод Рунге – Кутта с фиксированным шагом 2 Rkadapt Метод Рунге – Кутта с переменным шагом 3 Bulstoer Метод Булирша – Штера Более современный вариант не требует приведения ОДУ к каноническо му виду и построен по типу Given (Дано) – Odesolve (Реши систему ОДУ).

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

u := Odesolve ([vy, ] t, t1 [, step]).

Здесь:

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

t независимая переменная;

t1 конец интервала интегрирования (в этом случае началом интервала по умолчанию считается t0=0). Вместо этого может быть задан двухэлемент ный вектор-столбец начала t0 и конца t1 интервала.

step необязательный параметр, задающий количество равномерно от стоящих узлов на интервале интегрирования (по умолчанию step = 100), в которых строится решение.

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

В наборе этой программы участвуют два специальных символа: знак “жирного” равенства (набирается сочетанием клавиш Ctrl + = ) и знак производной (набирается сочетанием клавиш Ctrl + F7 ). Воз можен набор исходного дифференциального уравнения и в другой нотации производных, принятой в математике:

Ю. Л. Кетков, А. И. Кузнецов При наборе уравнений в этом формате используются либо соответствую щие кнопки в окне Calculus, либо клавишные комбинации Shift + / для обозначения первой производной и Shift + Ctrl + / для старших производных.

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

Для оценки точности результата интегрирования можно воспользовать ся тем фактом, что при заданных начальных условиях уравнение (2) имеет аналитическое решение z(t):

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

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

Таблица step rkxed rkadapt 10 2e-1 7e- 100 9e-5 9e- 1000 2e-7 2e- 10000 1e-13 5e- 100000 7e-13 5e- 1000000 8e-12 1e- Из этой таблицы следует, что дробление интервала интегрирования на 100–1000 отрезков обеспечивает точность, вполне приемлемую для решения большинства практических задач.

Версия MatLab-7 предлагает для решения “хороших” систем типа (1) то же три процедуры (см. табл. 3), выбор которых и подбор соответствующих управляющих параметров полностью возлагаются на пользователя.

Математика в высшем образовании 2005 № Таблица № Имя функции Пояснение 1 ode23 Метод Рунге – Кутта 2/3 порядка в модификации Богацки – Шампина 2 ode45 Методы Рунге – Кутта 4/5 порядка в модификации Дорманда – Принца 3 ode113 Метод Адамса в модификации Башфорта – Моултона Обращение к любой из функций odexxx может быть задано в одном из трех форматов:

[t, y] = odexxx(fun, [t0 t1], y0);

[t, y] = odexxx(fun, [t0 t1], y0, options);

[t, y] = odexxx(fun, [t0 t1], y0, options, p1, p2,...).

Здесь fun указатель на функцию вычисления правых частей, к которой про цедуры интегрирования обращаются либо с аргументами (t,y), либо с аргу ментами (t,y,p1,p2,... );





[t0,t1] интервал интегрирования (допускается t1t0);

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

options список управляющих параметров и их значений.

В простейшем случае для исследования поведения линейного гармони ческого осциллятора на интервале времени от 0 до 10 программа на языке MatLab может состоять из двух пользовательских функций:

Головная программа (pendulum.m):

global p1 p omega=0.5;

beta=0.2;

p1=-1/(omega^2);

p2=beta;

[t,y]=ode45(@dif1,[0 10],[0 1]);

plot(t,y(:,1),t,y(:,2)) grid on Функция вычисления правых частей (dif1.m):

function dy=dif1(t,y) global p1 p dy(1,1)=y(2);

dy(2,1)=p1*(p2*y(2)+y(1));

Дополнительные параметры p1 и p2 объявлены глобальными только для того, чтобы сократить затраты на их передачу в функцию вычисления пра вых частей (собственно, роль параметра p2 с таким же успехом могла играть Ю. Л. Кетков, А. И. Кузнецов переменная beta). Две последние строки в головной программе предназначе ны для визуализации графиков функций y1(t) и y2(t). Указатель на функцию вычисления правых частей (@dif1) может быть задан и в виде строки с име нем функции (‘dif1’).

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

set1 = odeset(’Stats’,’on’,’RelTol’,1e-4);

В приведенном выше примере параметру с именем Stats, управляющему выдачей итоговой информации о процессе интегрирования, присвоено значе ние on (включено) и параметру RelTol, управляющему относительной точно стью интегрирования, присвоено значение 1e-4. По умолчанию итоговая вы дача отключена (Stats = o), а относительная точность RelTol = 1e-3. Сфор мированное таким образом значение set1 должно быть указано в качестве четвертого параметра в обращении к функции odexxx:

[t,y]=ode45(@dif1,[0 10],[1 0],set1);

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

Таблица Параметр Значение Пояснение AbsTol Вектор абсолютной i-я компонента задает абсолютную точности {1e-6} точность по координате yi (t) RelTol Значение Относительная точность в первом относительной приближении определяет количество верных точности 1e-3 знаков в каждой координате yi (t) InitialStep Начальный шаг Определяется автоматически, интегрирования и его менять не рекомендуется MaxStep Максимальный шаг интегрирования, {0.1*|t1-t0|} Stats Режим итоговой Итоговая выдача включает количество выдачи {off} удачных попыток по изменению шага интегрирования (successful steps), число неудачных попыток (failed attempts), количество обращений к функции вычисления правых частей (function evaluations) и некоторые сведения, характерные для других процедур интегрирования Events Указание на Указание на m-файл, содержащий функцию, обработчик события отслеживающую поведение yi (t). Таким образом, в частности, можно прервать процесс интегрирования или переключиться на обработку другой системы уравнений Математика в высшем образовании 2005 № Самыми естественными вопросами, возникающими при исследовании конкретной динамической системы, являются следующие:

• каким методом интегрирования воспользоваться;

• какое время будет затрачено на процесс интегрирования;

• какая точность будет достигнута;

• как воспользоваться результатами интегрирования или пронаблюдать за ними?

На первый вопрос большинство авторов книг по MatLab’y рекомендуют начинать с процедуры ode45. Во-первых, среди методов Рунге – Кутта этот метод обладает повышенной точностью. Во-вторых, он действительно непло хо справляется с “хорошими” и “условно хорошими” системами.

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

Сложнее получить ответ на вопрос о точности. Надеяться на значения управляющих переменных AbsTol и RelTol особенно не стоит. Они только в первом приближении отвечают за ту или иную точность. Но когда количество шагов интегрирования становится достаточно большим, то, естественно, на бегают ошибки округления. Хорошо, если при каких-то значениях исходных данных вы знаете аналитическое решение. Тогда можно попытаться распро странить знания о полученной точности и на другие варианты системы. Од нако на практике известен довольно простой прием надо проинтегрировать систему “вперед” от начального момента времени t0 до конечного значения t1, а затем, внеся минимальные изменения в программу, произвести обратное интегрирование от t1 к t0 (естественно, не забыв обновить начальные усло вия). Расхождение между конечным результатом интегрирования “назад” и исходной начальной точкой дает вполне приемлемые сведения о точности интегрирования.

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

Для построения графика зависимости той или иной переменной yi (t) можно воспользоваться всеми массивами, полученными в результате интегрирова ния. По оси x можно отложить время t, а по оси y массив той или иной компоненты y(:,i). Т. к. в этих массивах содержится одинаковое число компо нент (пусть нам пока неизвестное), то график будет построен без каких-либо осложнений. Столь же просто вывести график на фазовой плоскости до статочно по оси x, например, отложить y1, т. е. массив y(:,1), а по оси y y2, т. е. массив y(:,2). А если нам по условиям задачи понадобится использо вать в дальнейших вычислениях значение некоторой функции в момент tk (например, yi (tk))? Тогда следует прибегнуть к одной из функций интерпо ляции, которых в MatLab’e предостаточно. Например, воспользоваться самой простой линейной интерполяцией:

yi = interp1(t, y(:,i),tk) Правда, если значение tk не принадлежит интервалу интегрирования [t0, t1], то результатом “интерполяции” будет значение NaN (Not a Number не число).

Ю. Л. Кетков, А. И. Кузнецов Если в левой части строки с обращением к одной из функций odexxx не за дать выходные параметры, то функция интегрирования сама строит график по всем вычисленным точкам траектории.

Теперь о некоторых сравнительных характеристиках процедур интегри рования, приведенных в табл. 3. Мы проделали все измерения на примере гармонического осциллятора, который в случае отсутствия затухания ( = 0) имеет простое аналитическое решение:

y1 = sin(t/), y2 = cos(t/).

В табл. 5 приведены результаты интегрирования затухающего осциллято ра ( = 0.5, = 0.2). С относительной точностью порядка 1e-3 они дают сходные результаты, но количество вычислений правой части и количество результатов в построенных массивах отличается довольно сильно. Графики функций, построенные по результатам интегрирования, на внешний взгляд, не различаются. Из табл. 5 видно, что для данной задачи наиболее эффек тивной оказалась процедура ode113.

Таблица ode23 ode45 ode y1(10) 0.006226 0.006330 0. y2(10) 0.010901 0.010880 0. Длина массива t 93 125 Вычисление правых частей 301 187 Второй эксперимент был связан с попыткой изменения абсолютной и от носительной точностей для наиболее популярной процедуры ode45. В связи с тем, что амплитуды функций y1, y2 отличаются незначительно, индивидуаль ные относительные точности задавались равными. Результаты эксперимента приведены в табл. 6 (менялась только абсолютная точность) и табл. 7 (меня лась только относительная точность).

Таблица AbsTol y1(10) y2(10) Вычисление правых частей 1e-1 0.008309 0.001785 1e-2 0.007452 0.005904 1e-3 0.006422 0.009851 1e-4 0.006320 0.010835 1e-5 0.006331 0.010882 1e-6 0.006330 0.010880 1e-7 0.006331 0.010883 1e-8 0.006330 0.010880 Таблица RelTol y1(10) y2(10) Вычисление правых частей 1e-1 0.007564 0.004812 1e-2 0.006296 0.010402 1e-3 0.006330 0.010880 1e-4 0.006344 0.010907 1e-5 0.006346 0.010908 1e-6 0.006346 0.010908 Математика в высшем образовании 2005 № Для анализа точности интегрирования был проделан эксперимент с неза тухающим осциллятором на достаточно большом периоде интегрирования при значениях параметров AbsTol и RelTol по умолчанию. Интервалы инте грирования оканчивались в моменты времени (2k + 1), где функция y должна принимать значение 1, а функция y2 обращаться в 0. Для k = 1 тести рованию подверглись все три функции ode23, ode45, ode113 (см. три первых объединенных строки табл. 8). Для больших k использовалась только проце дура ode45. Конечно, о точности интегрирования надо судить по функции y1, значение которой приблизительно за 1000 шагов интегрирования откло нилось от аналитики на 6%. Примерно такая же погрешность наблюдалась со смещением точки по оси t, в которой функция y2 обращалась в 0.

Таблица k y1 y2 Шагов Вычислений интегрирования правой части 1 0.990278 0.000992 95 (ode23) 0.999472 0.000887 117 (ode45) 0.998527 0.001655 60 (ode113) 2 0.99911 0.001496 181 4 0.998367 0.002710 309 8 0.996876 0.005134 565 15 0.994215 0.009358 1013 3. ЖЕСТКИЕ СИСТЕМЫ ОДУ В отличие от “хороших” систем ОДУ к “жестким” системам тяготеют ди намические процессы, протекание которых характеризуется одновременным присутствием компонент с резко различающимися частотными и амплитуд ными характеристиками. Такими свойствами обладают многие задачи хими ческой кинетики, физические процессы с резко выраженными нелинейностя ми и т. п. Использование стандартных методов интегрирования для решения жестких ОДУ страдает двумя недостатками:

• либо резко увеличивается объем вычислений из-за автоматического уменьшения шага интегрирования для достижения заданной точности (адаптивный метод начинает подстраиваться под наиболее сильно ос циллирующую или резко меняющуюся функцию);

• либо точность вычислений оставляет желать лучшего.

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

В качестве примера “жесткой” задачи приведем модель распространения пламени, предложенную Шампином и описанную в книге Моулера “Numerical Computing with MATLAB” (http://www.mathworks.com/moler).

Ю. Л. Кетков, А. И. Кузнецов Пусть y(t) радиус шара, охваченного пламенем. Скорость изменения y(t) состоит из двух слагаемых, одно из которых прямо пропорционально поверхности (y 2 ), охваченной пламенем, а второе обратно пропорционально объему шара (y 3 ):

y = y2 y3.

Начальный радиус шара y(0) = очень мал, и процесс исследуется на интервале времени [0, 2/]. У задачи Шампина имеется аналитическое реше ние, с которым удобно сравнивать решения, полученные различными метода ми интегрирования. Для построения графика аналитического решения при = 0.01 мы воспользовались следующей модификацией программы Молера:

y=dsolve(’Dy=y^2-y^3’,’y(0)=1/100’);

y=simplify(y);

pretty(y);

ezplot(y,0,200) grid on Эта программа с помощью расширения Symbolic Math Toolbox находит аналитическое решение y(t)=1/({\sf lambertw}(99 exp(99 - t)) + 1), где lambertw функция Ламберта, и строит график y(t) на интервале [0,200] (рис. 2).

Рис. 2. График изменения радиуса горящего шара при = Существенно более жесткой система становится при = 0.0001 (рис. 3).

Из приведенных графиков следует, что функция y(t) довольно резко воз растает в районе t = 1/ и сразу после этого выходит на стационар y(t) = 1.

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

Рис. 3. График изменения радиуса горящего шара при = Математика в высшем образовании 2005 № Таблица № Имя функции Пояснение 1 ode15s Метод конечных разностей переменного порядка в соче тании со схемой обратного дифференцирования (метод Гира), обеспечивает достаточно высокую точность ин тегрирования 2 ode23s Модифицированный метод Розенброка 2-го порядка, обеспечивает среднюю точность интегрирования со ско ростью выше, чем ode15s 3 ode23t Метод трапеций для интерполяции производных, реко мендуется для умеренно жестких систем, обеспечивает среднюю точность 4 ode23tb Комбинированный метод Рунге – Кутта в сочетании с правилом трапеций и схемы обратного дифференциро вания, обеспечивает среднюю точность Программа интегрирования для всех исследуемых методов имеет одина ковую структуру:

delta=0.0001;

F=inline(’y^2-y^3’,’t’,’y’) opts=odeset(’Stats’,’on’);

tic;

ode45(F,[0,2/delta],delta,opts);

toc grid on Однако результаты интегрирования отличаются довольно резко. Функция ode45, не предназначенная для работы с жесткими системами, спустя 24 сек(!) выдает график, напоминающий аналитическое решение (рис. 4):

Смущают два обстоятельства. Во-первых, затраченное время и количе ство обращений к вычислению правой части 22741. Наконец, уж очень подозрительна толщина графика в районе скачка и после выхода на стацио нар. Последнее несложно уточнить “под лупой”, использовав режим Zoom in (рис. 5).

Оказывается, что вместо абсолютно прямой линии после скачка (там 2 = y 3 = 1 и y = 0) функция ode45 продолжает выдавать какие-то ко y лебания (из-за них линия и кажется толстой).

Увеличенные графики функции y(t), полученные специальными функ циями интегрирования жестких систем, приведены на рис. 6. Их поведение не совсем точно отражает рекомендации табл. 8, почерпнутые из руковод ства по MATLAB. Самый “точный” метод (функция ode15s) не удержался от перерегулирования в окрестности выхода на стационар и показал время пе реключения, опережающее аналитику примерно на 3% (9680 вместо 10000).

Наиболее близким к аналитическому решению оказались результаты инте грирования с помощью функции ode23s (точность выхода на момент пере ключения порядка 0.2%). Конечно, по результатам интегрирования конкрет ной системы нельзя делать какие-либо обобщающие выводы, но и не следует Ю. Л. Кетков, А. И. Кузнецов Рис.4. Результат работы функции ode Рис. 5. Увеличенный фрагмент графика в окрестности t = Рис. 6. Фрагменты решения жесткой системы в окрестности скачка Математика в высшем образовании 2005 № слепо доверять рекомендациям руководств. Числа, расположенные в нижней строке рис. 6, фиксируют количество обращений к правой части (числитель дроби) и количество решений систем линейных алгебраических уравнений (знаменатель дроби). По ним можно судить о быстродействии той или иной функции интегрирования. По сравнению с функцией ode45 специализирован ные методы интегрирования жестких систем работают примерно в 1000 раз быстрее.

В системе MathCAD для решения умеренно жестких и жестких систем ОДУ рекомендуются функции, список которых приведен в табл. 10. Первые две функции не требуют задания якобиана системы, поэтому явное обра щение к ним может быть выполнено с помощью обращения, описанного в разделе 2:

u:=Radau(y0,t0,t1,M,D) u:=Bulstoer(y0,t0,t1,M,D) Для двух других функций обязательным требованием является задание якобиана системы:

u:=Stiff*(y0,t0,t1,M,D,J) Таблица № Имя функции Пояснение 1 Radau Метод RADAU5 для жестких систем 2 Bulstoer Стандартный метод Булирша-Штера для хороших и умеренно жестких систем 3 Stib Метод Булирша-Штера, модифицированный для жест ких систем 4 Stir Метод Розенброка для жестких систем Мы протестировали каждую из функций интегрирования жестких ОДУ на задаче о распространении пламени в умеренно жестком ( = 102 ) и жест ком ( = 104 ) вариантах. Так как исследуемое уравнение первого порядка, казалось бы естественным отступить от векторного задания функции, вычис ляющей правую часть системы. И для методов Radau и Bulstoer это предпо ложение оказалось справедливым:

Однако специализированные функции Stib и Stir потребовали соблюде ния более жестких правил оформления функций D и J:

Ю. Л. Кетков, А. И. Кузнецов Результаты тестирования сведены в табл. 11. Из таблицы следует, что для достижения хорошей точности решения приходится задавать довольно боль шое количество точек на интервале интегрирования. Разочаровали откро венно слабые возможности метода Radau, который в руководстве по системе MathCAD позиционирован как средство исследования жестких систем.

Таблица 2 = 10 = Radau M=10: грубое решение При всех M отказ от работы M=100: грубоватое решение M=1000: грубоватое решение M=10000: приемлемое решение Bulstoer M=10: очень грубое решение M=10: очень грубое решение M=100: приемлемое решение M=100: переполнение M=1000: хорошее решение M=1000: очень грубое решение M=10000: хорошее решение M=10000: хорошее решение Stib M=10: грубоватое решение M=10: грубоватое решение M=100: приемлемое решение M=100: хорошее решение M=1000: хорошее решение M=1000: отличное решение M=10000: отличное решение M=10000: отличное решение Stir M=10: грубое решение M=10: грубоватое решение M=100: грубоватое решение M=100: хорошее решение M=1000: хорошее решение M=1000: отличное решение M=10000: хорошее решение M=10000: отличное решение 4. НЕЯВНЫЕ ОДУ Неявное дифференциальное уравнение, не разрешенное относительно производной, задается уравнением F (t, y, y ) = 0.

Если оно аналитически неразрешимо относительно старшей производной, то возможность его решения предлагает только система MATLAB 7, в кото рой предусмотрена специальная функция ode15i:

[t,y] = ode15i(fun, [t0 t1], y0, yp0);

[t,y] = ode15i(fun, [t0 t1], y0, yp0, options);

[t,y] = ode15i(fun, [t0 t1], y0, yp0, options, p1,p2,...);

Смысл входных и выходных параметров тот же, что и для других функ ций odexxx. В пояснении нуждается лишь дополнительный аргумент yp начальное значение производной. Для неявного дифференциального уравне ния оно должно задаваться вместе с начальным моментом t0 и начальным значением функции y0 так, чтобы удовлетворялось условие F(t0,y0,yp0)= (его левая часть вычисляется функцией fun). Подобрать такие y0 и yp0 помо гает вспомогательная функция decic, обращение к которой имеет вид:

[y0mod,yp0mod] = decic(odefun,t0,y0,fixed_y0,yp0,fixed_yp0);

Математика в высшем образовании 2005 № Здесь:

y0, yp0 задаваемые пользователем гипотетические начальные значения, которые могут и не удовлетворять уравнению F(t0,y0,yp0)=0;

y0mod, yp0mod модифицированные начальные значения, вычисляемые функцией decic, которые достаточно точно удовлетворяют этому уравнению.

Процессом их вычисления можно управлять с помощью параметров xed y0 и xed yp0. Так, если задано xed y0=1, то y0mod=y0. Если же xed y0=0, то y0mod может отличаться от y0. Аналогичную роль для вы числения yp0mod играет параметр xed yp0. Если мы имеем дело с системой уравнений, когда переменные y и yp являются векторами, xed y0 и xed yp также задаются как вектора и соглашение об их значениях 0 и 1 относится по отдельности к каждой их компоненте. Если любой из этих управляющих векторов задать в виде пустой матрицы [ ], это равносильно тому, что все его компоненты нулевые. По сути дела, функция decic находит приближенное решение нелинейного уравнения.

Применим эти соображения к упоминавшемуся ранее колебательному “контуру с железом” [3, 4]. Катушка этого контура снабжена железным сер дечником, магнитный поток в котором является нелинейной функцией (i) от протекающего через катушку тока i. На практике эту функцию достаточно хорошо аппроксимируют выражением:

(i) = A · arctg (w · i/S) + B · w · i/S.

Здесь A, B и S некоторые константы, зависящие от материала сердечника, а w число витков провода в катушке.

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

2 + Bq + q 2 E = 0.

F (t, q, q ) = A · ln 1 + q Здесь q заряд на обкладках конденсатора;

q = dq/dt ток в контуре;

E значение энергии, зависящее от начальных условий.

Вообще говоря, при заданных начальных значениях заряда q(0) и его производной qp(0) константа E достаточно просто вычисляется [4]. Но мы поступим иначе зададим правильное значение q(0) = 0 и значение E, со ответствующее значению qp(0) = 1. Но вместо правильного значения началь ной производной при обращении к вспомогательной функции decic укажем “неправильное” значение qp(0) = 2.

Головная программа ferrum.m имеет вид global A B E A=1.4;

B=25;

E=A*log(2)+25;

t0=0;

q0=0;

qp0=2;

[q0,qp0]=decic(@ABX,t0,q0,1,qp0,0) Ю. Л. Кетков, А. И. Кузнецов [t,y]=ode15i(@ABX,[0 100],q0,qp0);

plot(t,y) grid on Функция вычисления F (t, q, qp) может быть оформлена так:

function res=ABX(t,q,qp) global A B E res=A*log(1+qp^2)+B*qp^2+q^2-E;

Результаты работы служебной программы decic приведены ниже:

q0 = qp0 = Результаты последующего интегрирования с помощью функции ode15i совпадают с данными, приведенными в [4].

Следует отметить, что попытка задать гипотетическое значение qp(0), меньшее 1, приводит к сбою в работе программы decic, которой не удается подобрать нужное решение с требуемой точностью.

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

[t, y] = ode15i(fun, [t0 t1], y0, yp0, set1), где значение set1 формируется с помощью функции odeset, аргументами ко торой являются пары вида (’Имя’, значение). В нашем случае понадобится параметр ’Events’, а его значением будет указатель на функцию, определя ющую три выходных управляющих параметра. Имена выходных аргументов могут быть произвольными, но их последовательность имеет функциональ ную предопределенность:

function [v,i,d]=qq(t,q) v=q-4;

i=1;

d=1;

Контролируемое событие заключается в том, что значение первого выход ного параметра v обращается в нуль при условии возрастания этой величины (признак d=1) и процесс интегрирования прекращается при наступлении ука занного события (признак i=1).

Обращение к ode15i в нашем случае может выглядеть так:

set1 = odeset(’Events’,@qq);

[t,q] = ode15i(@ABX,[t0 t1],q0,qp0,set1);

Результат работы модифицированной таким образом программы ferrum.m приведен на рис. 7.

Справка по функции ode15i содержит дезинформацию по поводу оформ ления функции qq реакции на событие прерывания процесса интегрирова ния. Там указано, что эта функция должна получать три входных параметра (t,q,qp). Если следовать этому указанию, то MATLAB выдает сообщение о том, что аргументов слишком много.

Математика в высшем образовании 2005 № Рис. 7. Интегрирование до выполнения условия q = 5. РЕШЕНИЕ ОДУ С ЗАПАЗДЫВАЮЩИМ АРГУМЕНТОМ Возможность исследования дифференциальных уравнений с запаздыва ющим аргументом предусмотрена только в системе MatLab-7. Для решения систем ОДУ вида (5) с набором положительных фиксированных задержек 1, 2,..., k предназначена функция dde23, обращение к которой в простей шем случае имеет вид:

res = dde23(@f_r, [?1, ?2, :, ?k], @f_h, [t0, t1]);

здесь:

@f_r указатель на функцию вычисления правых частей уравнения (5), при обращении к которой функция dde23 передает три аргумента (теку щее время t, вектор текущего решения y и вектор запаздывающих значений yd = [y(t 1 ), y(t 2 ),..., y(t k )]);

[1, 2,..., k ] вектор постоянных запаздываний;

@f_h указатель на функцию, генерирующую значения запаздывающих аргументов y(t) при t t0. В частном случае этим параметром может быть вектор-столбец, если начальные значения являются константами;

[t0, t1] интервал интегрирования (t0 t1);

res структура, в которой выдается решение. Е поле res.x содержит е последовательность времен ti, в которых найдены приближенные значения вектора y. Найденное решение y(ti ) выдается в поле res.y.

В состав демонстрационных примеров версии MatLab-7 входит файл ddex1.m, который мы используем в качестве иллюстрации работы функции dde23. Система исходных уравнений имеет вид y1’(t) = y1(t-1) y2’(t) = y1(t-1)+ y2(t-0.2) y3’(t) = y2(t) Для вычисления правых частей можно воспользоваться функцией f r.m (имя функции может быть и другим):

function dydt=f_r(t,y,y_d) y1=y_d(:,1);

y2=y_d(:,2);

dydt=[y1(1);

y1(1)+y2(2);

y(2)];

Ю. Л. Кетков, А. И. Кузнецов Предыстория при t=0 задана тремя одинаковыми константами:

y1(t)= y2(t)= y3(t)= Это позволяет либо воспользоваться вектор-столбцом ones(3,1) в обращении к функции dde23, либо написать файл предыстории f h.m:

function s=f_h(t) s=ones(3,1);

Так как вектор запаздываний представлен двумя компонентами [1, 0.2], то для интегрирования нашей системы на интервале [0,1] можно воспользовать ся одним из двух следующих обращений:

res = dde23(@f_r,[1, 0.2],@f_h,[0, 1]);

res = dde23(@f_r,[1, 0.2],ones(3,1),[0, 1]);

График найденного решения со всеми тремя кривыми y1 (t), y2 (t) и y3 (t) мож но построить с помощью оператора plot(res.x,res.y). Результат его работы при веден на рис. 8. Если нас интересует только одна из графических зависимо стей, например, y2 (t), то это можно сделать следующим образом:

plot(res.x,res.y(2,:)) Рис. 8. Система ОДУ с запаздыванием Набор управляющих параметров функции dde23 изменяется с помощью вспомогательной функции ddeset.

6. РЕШЕНИЕ КРАЕВЫХ ЗАДАЧ Напомним вкратце постановку краевых задач и схему их решения на при мере дифференциального уравнения второго порядка [5]:

y = f (x, y, y ). (9) Краевые условия на отрезке интегрирования [a, b] обычно задают в одной из трех форм:

y(a) =, y(b) =, (10) y (a) =, y (b) =, (11) y (a) + k1 · y(a) =, y (b) + k2 · y(b) =. (12) Математика в высшем образовании 2005 № Интервал интегрирования [a, b] разбивают на n равных отрезков, в концах которых значения производных заменяют их конечно-разностными аналога ми. В результате исходное дифференциальное уравнение сводится к n 1 ал гебраическому уравнению относительно значений искомой функции в узлах y1, y2,..., yn1. Присоединяя к ним граничные условия, в которых при необ ходимости значения производных также заменены разностными формулами, мы получаем n + 1 уравнение относительно значений неизвестной функции.

Решение этой системы алгебраических уравнений дает приближенное реше ние краевой задачи.

В отличие от задачи Коши, которая при заданных начальных условиях, как правило, имеет единственное решение, краевая задача может не иметь решения вообще, иметь конечное число решений или иметь бесконечное число решений. Последнее можно наблюдать на довольно простом примере: пусть уравнению (9) удовлетворяет функция y = g(x) и граничные условия заданы в формате (11). Тогда любая функция вида y = g(x) + const тоже является решением краевой задачи при тех же граничных условиях.

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

sol=bvp4c(@f_e, @f_b, [a b], y_init);

здесь:

f_e указатель на функцию вычисления правых частей системы ОДУ, сведенной к виду (1);

f_b указатель на функцию ограничений, приведенную к формату f y(a), y(b) = 0;

[a b] интервал интегрирования;

y_init начальное приближение к искомому решению;

sol структура с полями sol.x и sol.y, в которых получается решение.

При вызове функции вычисления правых частей ей передаются текущие значения независимой переменной x и вектор текущих значений y. Она, в свою очередь, должна сформировать вектор значений dy/dx в соответствии с исходными ОДУ.

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

Как правило, для этой цели используется вспомогательная функция bvpinit с двумя векторными аргументами:

y_init=bvpinit(linspace(a,b,n),y0);

здесь:

linspace функция, разбивающая интервал [a b] на n равных отрезков и генерирующая вектор узловых точек начальной сетки. Значение n выбира ется пользователем, чаще всего задают n = 10;

y0 вектор начальных констант, определяющий стартовые значения ис комых функций на интервале [a b]. Первая константа y0(1) используется в качестве начального приближения к первой функции y1 во всех внутренних узлах сетки, вторая константа y0(2) в качестве начального приближения Ю. Л. Кетков, А. И. Кузнецов ко второй функции y2 и т. д. В случае единственного решения краевой задачи выбор n и вектора y0 особого влияния на ход решения не оказывают. Однако краевая задача может иметь не одно решение, и тогда приходится перебирать несколько вариантов задания начальных значений.

В качестве первого примера краевой задачи поступим следующим об разом. Предположим, что решаемая задача имеет аналитическое решение y = eax cos(x), которое мы будем рассматривать на интервале [0, 3.5]. За таким решением очень удобно наблюдать, да и граничные условия на функ цию легко задаются. Остается только составить подходящую систему диф ференциальных уравнений второго порядка:

y1 = y = eax cos(x), y1 = y2 = aeax cos(x) eax sin(x), y2 = (a2 1)y1 + 2aeax sin(x).

Выберем a = 0.2 и зададим граничные условия y(0) = 1 и y(3.5) = 0.

Головная программа решения нашей краевой задачи может выглядеть так:

global a a=0.2;

y_init=bvpinit(linspace(0,3.5*pi,10),[0 1]);

sol=bvp4c(@f_e1,@f_b1,y_init);

plot(sol.x,sol.y);

grid on В качестве начального приближения во внутренних точках интервала ин тегрирования выбраны константные функции y1(x) 0 и y2(x) 1, что совсем не похоже на граничные условия.

Функция вычисления вектора правых частей системы ОДУ (f_e1.m) име ет вид function dydx=f_e1(x,y) global a dydx=[y(2);

y(1)*(a^2-1)+2*a*exp(-a*x)*sin(x)];

Функция соблюдения граничных условий (f_b1.m) имеет вид function res=f_b1(ya,yb) res=[ya(1)-1;

yb(1)];

Результат работы головной программы приведен на рис. 9. Никакие изме нения в выборе узлов начальной сетки и в задании начальных приближений на окончательный результат никакого влияния не оказывают.

В качестве примера второй краевой задачи, характерной наличием двух решений, мы воспользуемся материалом из файлов справки системы MatLab 7. Там демонстрируется решение уравнения y + |y| = 0 с граничными усло виями y(0) = 0 и y(4) = 2. В качестве начального приближения интервал интегрирования разбивается на 5 равных отрезков, в узлах которого заданы константные функции y1(x) 0 и y2(x) 1. Головная программа решения второй краевой задачи имеет вид:

Математика в высшем образовании 2005 № Рис. 9. Единственное решение первой краевой задачи Рис. 10. Первое решение краевой Рис. 11. Второе решение краевой задачи задачи y_init = bvpinit(linspace(0,4,5),[0 1]);

sol = bvp4c(@twoode,@twobc,y_init);

plot(sol.x,sol.y);

grid on Функции вычисления вектора правых частей (twoode.m) и соблюдения гра ничных условий (twobc.m) таковы:

function dydx = twoode(x,y) dydx = [y(2);

-abs(y(1))];

function res = twobc(ya,yb) res = [ya(1);

yb(1) + 2];

График первого решения с начальными условиями [0 1] приведен на рис. 10.

Однако, стоит заменить начальные константы на [-10 0], как находится дру гое решение (рис. 11).

Для построения более гладкого графика найденных решений можно вос пользоваться сервисной программой deval, предварительно разбив интервал интегрирования на большее количество равных отрезков (по умолчанию 100) x = linspace(0,4);

y = deval(sol,x);

Ю. Л. Кетков, А. И. Кузнецов Средства, предлагаемые системой MathCAD для решения краевых задач, напоминают подходы для решения систем ОДУ. Более современный подход в формате Given – Odesolve кажется еще более простым, чем аналогичная функ ция bvp4c в системе MatLab. Он даже не требует задания начального при ближения к искомому решению:

Однако таким способом находится только первое решение краевой зада чи, приведенное на рис. 10. Для поиска второго решения приходится слегка изменить начальное условие:

Хотя внешне график нового решения почти ничем не отличается от кри вой, приведенной на рис. 11, но вс же этим методом следует пользоваться е очень осторожно, т. к. есть возможность получить константную погрешность (т. е. не уменьшается с повышением количества шагов метода), которую будет практически невозможно учесть при использовании полученного решения.

Второй способ решения краевых задач, широко рекламируемый в ряде книг [6, 7], носит название метода “пристрелки”. Он заключается в ручном или автоматическом подборе недостающих начальных условий таких, при кото рых с определенной точностью будут выполнены граничные условия на кон це интервала интегрирования. Для уравнений второго порядка управление “пристрелкой” выглядит довольно просто. Сводим краевую задачу к задаче Коши, задав каким-то образом недостающее начальное условие. Находим ин тегральную кривую, явно обратившись к той или иной функции интегрирова ния. Оцениваем разницу между получившимся условием на конце интервала и заданным в краевой задаче. Затем корректируем искусственное начальное условие и, многократно решая задачу Коши, пытаемся минимизировать по лучающуюся разницу. Когда исследуемая система ОДУ имеет более высокий порядок, приходится минимизировать отклонение на конце интервала, управ ляя двумя или более параметрами. И в режиме ручного подбора процедура “пристрелки” может затянуться надолго, да и достигаемая точность может оказаться неудовлетворительной.

Для автоматизации процедуры подбора недостающих начальных условий в MathCAD’e можно воспользоваться функцией командной строки sbval:

S:=sbval(v,t0,t1,D,f0,f1) Здесь:

S вектор, в который функция sbval помещает оптимальные недостаю щие начальные условия;

v исходное приближение для недостающих начальных условий, которое улучшается в процессе итеративного поиска;

Математика в высшем образовании 2005 № t0, t1 начало и конец интервала интегрирования;

D(t, y) вектор-функция, описывающая правые части системы ОДУ;

f 0(t0, v) вектор-функция, которая представляет полный набор началь ных условий. Константы в этом векторе соответствуют условиям, заданным в краевой задаче, а элементы вектора v недостающим начальным условиям.

Последовательность этих компонент вектора v определяет порядок выдачи результатов в векторе S;

f 1(t1, y) вектор-функция разности между приближенным решением y и заданными условиями (функция sbval пытается обратить в 0 эти разности).

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

После определения недостающих начальных условий краевая задача све лась к задаче Коши:

Найденное численное решение можно сравнить с аналитическим y = = 4/(x + 2) 1.

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

7. ВЫВОДЫ Мы отдаем себе отчет в том, что конкретные примеры задач, приведен ные в настоящей статье, не могут быть достаточным основанием для выне сения окончательного вердикта в пользу тех или иных инструментальных Ю. Л. Кетков, А. И. Кузнецов средств численного исследования систем ОДУ. Более того, личные привязан ности авторов по поводу систем MatLab и MathCAD не совпадают. Однако мы считаем, что следующие рекомендации носят достаточно объективный характер:

• для исследования “хороших” систем ОДУ пакеты MatLab и MathCAD предлагают примерно эквивалентные средства интегрирования. Некото рое преимущество системы MathCAD заключается в возможности набора ОДУ в формате, принятом в математической литературе, и в отсутствии необходимости разрешения уравнений относительно старшей производной;

• средства обработки жестких систем ОДУ в пакете MatLab несколько пред почтительнее. Они не требуют обязательного задания якобиана системы;

• Если необходимо получить гладкое решение, то MathCAD позволяет его получить в формате given –odesolve. В MatLab’e в этом случае необходимо обратиться к одной из функций сглаживания (например, deval);

• средства исследования неявных ОДУ, не разрешимых относительно стар шей производной, доступны только в пакете MatLab-7;

• только пакет MatLab включает средства для исследования систем ОДУ с запаздыванием;

• возможности обоих пакетов по решению краевых задач достаточно близ ки, но MathCAD позволяет исследовать системы ОДУ с дополнительными условиями во внутренних точках интервала интегрирования.

ЛИТЕРАТУРА 1. Дьяконов В. П. MATLAB 6/6.1/6.5 + Simulink 4/5. Основы применения. Полное руко водство пользователя. М.: Солон, 2002. 768 с.

2. Кудрявцев Е. М. GPSS World. Основы имитационного моделирования различных си стем. М.: ДМК Пресс, 2004. 320 с.

3. Андронов А. А., Витт А. А., Хайкин С. Э. Теория колебаний. М.: Физматгиз, 1959.

4. Кетков Ю. Л., Кетков А. Ю., Шульц М. М. MATLAB 7: программирование, численные методы. СПб.: БХВ-Петербург, 2005. 752 с.

5. Березин И. С., Жидков Н. П. Методы вычислений. Т. II. М.: Гос. изд. физ.-мат. ли тературы, 1959. 620 с.

6. Кирьянов Д. В. Mathcad 12. СПб.: БХВ-Петербург, 2005. 576 с.

7. Гурский Д. А., Турбина Е. С. Вычисления в Mathcad 12. СПб.: Питер, 2006. 544 с ORDINARY DIFFERENTIAL EQUATIONS:

MATLAB VERSUS MATHCAD Yuli L. Ketkov, A. I. Kuznetsov The article describes the features of two automatic mathematical calculations sys tems MathCAD and MatLab in the eld of solving ordinary dierential equations and systems of the equations. A comparison of abilities of implemented methods for solving various classes of dierential equations is given.

Keywords: dierential equations, numerical methods, sti dierential equations, Mat Lab, MathCAD.



 

Похожие работы:





 
2013 www.netess.ru - «Бесплатная библиотека авторефератов кандидатских и докторских диссертаций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.