Учебные материалы по математике | Численное решение задачи коши для обыкновенных дифференциальных уравнений | Matematiku5
Вузы по математике Готовые работы по математике Как писать работы по математике Примеры решения задач по математике Решить задачу по математике online

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


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

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

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

(1)

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

(2)

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

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

(3)

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

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

(4)

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

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

(5)

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

(6)

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

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

(7)

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

(8)

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

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

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

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

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

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

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

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

 

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

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

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

(9)

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

(10)

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

Метод Эйлера

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

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

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

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

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

,

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

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

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

(11)

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

И вообще

(12)

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

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

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

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

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

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

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

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

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

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

(13)

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

Рисунок 5

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

(14)

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

Рисунок 6.

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

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

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

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

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

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

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

(15)

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

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

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

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

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

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

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

(16)

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

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

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

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

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

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

(17)

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

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

В краевых задачах (в отличие от задачи Коши, в которой все дополнительные условия задаются в начальной точке) дополнительные условия задаются на концевых точках отрезка интегрирования (краях отрезка — отсюда и название «краевые задачи»).

Иногда такие задачи называют задачами с граничными условиями или просто «граничными задачами».

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

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

Дифференциальное движение материальной точки описывается уравнениями Ньютона

Здесь V{u, v) –вектор скорости, — вектор сил.

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

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

2)  в начальной точке известна только величина скорости (не известен угол a), но задана точка B, в которую должно попасть материальное тело.

В первом случае, зная все условия в начальной точке, можно решать задачу «пошагово» (например, методом Эйлера), получив после выполнения всех шагов искомую траекторию

Во втором случае такая «пошаговая» процедура невозможна. Чтобы попасть в заданную точку B, надо определить неизвестный угол a, но для этого надо решить задачу для всей траектории. В реальной жизни эта задача решается метод «пристрелки», когда многократно подбирается этот угол (корректируя «недолет–перелет»).

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

Пример краевой задачи

Рассмотрим прямоточный химический реактор (в виде трубы постоянного сечения).

На вход в реактор со скоростью v подается реагент с начальной концентрацией .

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

,

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

Приняв для диффузионного потока закон Фика (диффузионный поток пропорционален градиенту концентрации,

,

(коэффициент пропорциональности D называется коэффициентом диффузии), и считая скорость v постоянной, получим

.

В качестве граничных условий примем на левом конце (x=0):

,

А при будем считать, что реакция полностью закончилась, то есть

.

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

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

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

Положив

, ,

И подставляя это в дифференциальное уравнение, получим

,

Или

(18)

.

Здесь введены обозначения

и .

Обратите внимание, что величины P и R являются безразмерными (как и все дифференциальное уравнение).

Граничные условия в безразмерном виде запишутся как

(19)

и

(20)

.

Такая безразмерная форма не только более удобна для численного решения (по крайней мере, число безразмерных параметров — их в задаче всего два — меньше, чем число исходных физических параметров, x, D, v, C0, D, k), но и более физична.

Общие соображения о численном решении краевых задач

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

,

Здесь i — номер узла разностной схемы.

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

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

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

или или

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

В таких случаях мы советуем использовать пакеты типа COMSOL Multiphysics, в которых решение задачи выглядит более естественно, требуя от пользователя только понимания постановки задачи.

Ниже показано, как решить задачу (18-20) в COMSOL Multiphysics.

1. Сначала естественно выбирается физическое приложение, к которому относится наша задача. Это стационарная (steady-state) одномерная (1D) задача из области Конвекции и диффузии (Convection and Diffusion).

2. Затем для удобства можно задать величину (опция меню Options->Constants) параметров P и R уравнения

3.Затем в меню Draw создается область решения. В данном случае это всего лишь отрезок прямой линии [0;1].

4. В окне установки физических параметров уравнения (Subdomain Settings) требуется ввести значения параметров нашей конкретной задачи

Обратите внимание на подсказку в верхней части окна, где в общем виде записано уравнение вашей физической задачи. Сопоставьте его со своим уравнением (18) и введите коэффициенты, как показано на рисунке.

4. Теперь надо задать граничные условия. На левом конце зададим начальную концентрацию (в безразмерном виде начальная концентрация равна 1), а на левом конце зададим условие свободно выхода продукта (convection flux)

4. После этого достаточно нажать кнопку (Solve) и вы увидите решение своей задачи

Замечание: В качестве семестровой работы можно выбрать следующее:

а) найти задачу химической кинетики для реальной системы реакций и решить ее в MathCad.

б) Выбрать (по согласованию с преподавателем) функции MathCad решения дифференциальных уравнений, не рассмотренных в этой лекции (можно обратиться к электронным книгам, поставляемым с пакетом, например к книге ODE Solve Blocks)

б) выбрать любой из примеров модуля Chemical Engineering пакета COMSOL Multiphysics и внести в него изменения по согласованию с преподавателем.

Перечень контрольных вопросов

1.  Различие между краевой задачей и задачей Коши для обыкновенных дифференциальных уравнений.

2.  Умение составлять систему кинетических уравнений для произвольной системы реакций.

3.  Основная идея численного решения задачи Коши (метод Эйлера).

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

5.  Понятие многошаговых методов.

6.  Понятие неявных методов.

7.  Общие принципы выбора метода для различных задач. Выбор метода решения в MathCad.

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

9.  Общее понятие о разностных методах краевой задачи.

10.  (Не обязательно) Решение одномерной краевой задачи массопереноса с диффузий в пакете COMSOL.

Наташа

Автор

Наташа — контент-маркетолог и блогер, но все это не мешает ей оставаться адекватным человеком. Верит во все цвета радуги и не верит в теорию всемирного заговора. Увлекается «нефрохиромантией» и тайно мечтает воссоздать дома Александрийскую библиотеку.

Распродажа дипломных

 Скидка 30% по промокоду Diplom2020

А ты боишься COVID-19?

 Пройди опрос и получи промокод