Численное решение задачи коши для обыкновенных дифференциальных уравнений

Лекция. Численное решение задачи Коши для обыкновенных дифференциальных уравнений

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

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

При определенных начальных условиях:

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

Для простейшей реакции разложения вещества скорость реакции пропорциональна первой степени концентрации (реакция первого порядка) и потому можно записать

где k называется константой скорости реакции, а — начальная концентрация.

Это уравнение в силу своей простоты легко решается (без привлечения численных методов):

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

Пусть порядок реакции будет равен n. Кроме того, рассмотрим реакцию с выделением тепла (тепловой эффект равен H). Если реактор помещен в термостат с температурой , то количество тепла, передаваемое от него к термостату, можно считать пропорциональным разности температуры T в реакторе и температуры . Если принять, что скорость выделения тепла пропорциональна скорости реакции, то для температуры получим дифференциальное уравнение

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

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

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

Обозначая концентрации веществ A, M, B, C, D и E через С1, С2, C3, C4, C5 и С6, можно записать систему дифференциальных кинетических уравнений для этой схемы реакций в виде:

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

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

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

На следующем рисунке приведены результаты расчет концентраций (их зависимости от времени) для всех веществ.

Обратите внимание на то, что результирующая концентрация продуктов (веществ B, D, E) равна в сумме концентрации исходного вещества A (0.4+0.4+0.2=1), что является дополнительной проверкой выполнения общего баланса вещества.

Как видите, запись системы дифференциальных уравнений в MathCad мало отличается от лычной математической записи на бумаге. Надо только эти уравнения и начальные условия заключить внутрь блока Given …..Odesolve.

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

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

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

Постановка задачи

Требуется численно решить задачу Коши (с начальными условиями) для системы обыкновенных дифференциальных уравнений:

Все последующее изложение будем вести для одного уравнения:

Все сказанное без труда (и практически без каких-либо изменений) переносится на случай нескольких уравнений.

Метод Эйлера

Итак, требуется найти решение задачи(10).

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

Кривая изображена пунктирной линией, поскольку она нам заранее неизвестна.

Основная идея численных методов решения задачи Коши содержится в методе Эйлера.

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

,

что соответствует построению линейной функции, касательной к искомой функции y(x) в точке ().

Поэтому, если взять достаточно близко к и обозначить , то

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

Теперь решение (хотя и приближенно) известно в точке (). Поэтому можно повторить все предыдущие рассуждения и найти решение в следующей точке

И вообще

Можно повторять эти вычисления, получив последовательность значений , которая будет приближенно представлять решение данной задачи Коши.

Это и есть метод Эйлера.

Очевидно, чем меньше h (который называется шагом интегрирования), тем точнее будет решение.

Более того, можно заметить, что в формуле Тейлора ошибка e на каждом шаге имеет порядок . Если на всем отрезке интегрирования длиной L выполнить n шагов, то суммарная ошибка E будет суммой ошибок на всех отдельных шагах, то есть

То есть метод Эйлера имеет точность порядка O(h). Точную оценку получить затруднительно (она зависит от вида функции в правой части дифференциального уравнения), однако, полученная качественная информация оказывается очень полезной.

Уточнения метода Эйлера (методы второго порядка точности)

Чем обусловлена ошибка метода Эйлера?

Если бы производная от искомой функции не изменялась бы на отрезке , то касательная в точке () точно бы совпадала с интегральной кривой, и метод Эйлера дал бы точное решение. Значит, ошибка связана с тем, что указанная производная на отрезке не постоянна. Поэтому естественно попытаться как-то учесть это изменение. Самым простым способом учесть такое изменение является попытка вычислить некоторое среднее значение этой производной на отрезке .

Это можно сделать несколькими способами.

а) вычислим среднее арифметическое значение производной в точках () и ():

Эти рассуждения проиллюстрированы на рисунке 5.

Рисунок 5

б) Можно также в качестве среднего значения производной на отрезке вычислить ее значение в промежуточной (например, средней) точке отрезка, то есть

Эти рассуждения проиллюстрированы на рисунке 6.

Рисунок 6.

Можно показать, что оба последних метода имеют точность порядка , то есть если уменьшить величину шага h в 10 раз, то погрешность решения уменьшится в 100 раз.

Метод Рунге-Кутта четвертого порядка

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

(м-д Эйлера 1-го порядка)

(1-й м-д Эйлера 2-го порядка)

(2-й м-д Эйлера 2-го порядка)

Теперь можно легко поверить (опуская все выкладки), что метод 4-го порядка имеет вид

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

Многоточечные методы (методы Адамса)

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

Этот недостаток исправляется в так называемых многошаговых методах Адамса.

Чтобы получить представление об этих методах, запишем метод Адамса с двумя предыдущими точками (рисунок 7).

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

В итоге получим

Рисунок 7. Простейший многошаговый метод Адамса.

Уточнения одношаговых методов Рунге-Кутта связаны с добавлением внутренних точек на отрезке . Точно также уточнения многошаговых методов связаны с добавлением нескольких предыдущих точек. Из-за громоздкости вычислений здесь эти формулы не приводятся (их можно н6айти, например, в книге Хемминга [].

Неявные методы

Все предыдущие методы (и методы Рунге-Кутта и методы Адамса) явно вычисляют значение в точке по информации в предыдущих точках (одной или нескольких). По этой причине их называют явными методами.

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

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

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

Можно точно также построить (по аналогии с приведенными ранее рассуждениями) неявные методы Рунге-Кутта более высокого порядка точности. Можно также построить неявные многошаговые методы.

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

Метод «прогноза–коррекции»

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

а) прогноз

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

б) коррекция

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

На практике, достаточно выполнять три-четыре итерации, чтобы получить нужную точность.

Понятие о неустойчивости численных методов

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

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

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

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

Для иллюстрации рассмотрим задачу

Эта задача имеет точно решение .

Однако решение этой задачи при произвольных начальных условиях имеет вид (проверьте, решая уравнение разделением переменных):

Поэтому для начального условия (для любого малого ) получим решение, которое стремится к нулю при , а для начального условия получим неограниченно растущее решение (рисунок 8).

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

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

где вектор сил реакции поверхности на шарик.

Рисунок

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

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

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

Для иллюстрации рассмотрим задачу

,

точное решение которой имеет вид . В указанном выше смысле это решение устойчиво относительно начальных условий.

Применим теперь для решения этой задачи численный метод ():

Для данной задачи . Поэтому метод в данном случае дает

Эта формула порождает некоторую последовательность , поведение которой может быть полностью проанализировано []:

.

Видно, что при (при любом h)_первое слагаемое стремится к нулю, а второе, осциллируя, стремится к бесконечности.

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

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

Понятие о «жестких» дифференциальных уравнениях

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

Понятие жесткости можно проиллюстрировать на примере решения уравнения

Точным решением которого является

Применим (для простоты рассуждений) метод Эйлера, получим последовательность

,

где слагаемое аппроксимирует в точном решении.

Величина быстро убывает с ростом x, и решение, начиная с некоторого x, мало отличается от единицы. При этом интуитивно кажется, что для интегрирования можно взять достаточно большой шаг (поскольку решение почти не меняется). Однако из последнего соотношения видно, что при величина будет расти, свидетельствуя о неустойчивости.

Таким образом, хотя слагаемое при больших x практически не вносит никакого вклада в решение, при численном решении его приходится аппроксимировать очень точно и выбирать h очень малым.

Реально описанная ситуация встречается при решении задач химической кинетики, описывающих систему реакций, характерные времена которых сильно отличаются (такое различие может достигать нескольких порядков) [].

В целом «жесткие» системы требуют применения специальных методов (это, как правило, неявные методы).

В Mathcad для решения «жестких» систем предлагаются методы BDF, Radau, Stiffb, Stiffr.

Метод AdamsBDF сам определяет, является ли система жесткой, и в этом случае вызывается метод BDF, если же система не жесткая, то вызывается обычный метод Adams.

Все эти методы являются «неявными» методами, о которых речь шла выше.

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

Совет: Загляните в Справку MathCad (разделу Calculus and Differntial Equations в шпаргалках (QuickShets)). Там вы можете найти достаточно много тем для семестровых работ. Еще лучше, если у вас к MathCad приложены электронные книги (E-books), среди которых имеется книга DIFFERENTIAL EQUATIONS SOLVE BLOCK, в которой вы можете найти массу интересных инженерных приложений (в том числе и пример с «жесткими» дифференциальными уравнениями — проблему Хайриса)