Гурский Д., Турбина Е. - Вычисления в MathCad 12 (1077322), страница 104
Текст из файла (страница 104)
Аx5l. у :=bulstoer(yO,0,10,10~=9.32674444016026,Sys , 10,0.l)У5 1 = 9.99770030614651last(x ( l > )=9lastG ( l ) ) = 9x :=bulstoer(yO, 0,10,10" 3,Sys, 10000, O.l) у := bulstoer(yO,0,10,10"x 5 j = 9.32674444016026last(x ( l > ) = 351616,Sys, 10000, O.l)У5 { = 9.99770030614651'l a S t(y ( l > ) = 50Как видно, точность решения, полученного методом Булирша-Штера, зависит не отколичества разбиений всего интервала поиска, а от заданной величины TOL, котораяопределяет внутреннее количество шагов каждого промежутка [х;х+Н].14.2.5. Жесткие системы ОДУПонятие о жесткости (Stiff) системы дифференциальных уравнений появилось в середине XX века в связи с решением систем, описывающих протекание катализируемыххимических реакций. Оказалось, что такие, казалось бы, эффективные и абсолютнонадежные численные методы, как алгоритм Рунге-Кутта (да и все другие разработанные на том этапе явные алгоритмы) не способны решать системы подобного рода.
Этонеприятное открытие послужило толчком для развития обширнейшего и важнейшегонаправления в математике численных методов, и сейчас алгоритмы, способные решатьжесткие системы, играют на практике куда большую роль, чем классические явныеметоды.Точного определения жесткой системы ОДУ в общем не существует. Обычно жесткойсчитается система, описывающая сильно разномасштабные процессы. Трудности, связанные с решением таких систем явными численными методами, объясняются тем, чтоусловие устойчивости требует брать шаг слишком мелким. А это означает, что для того,чтобы найти корректный результат, алгоритм придется прокрутить миллионы илидаже миллиарды раз, что может быть непосильной задачей даже для современного компьютера.4 7 0 • Глава 14.
Дифференциальные уравненияРассмотрим следующее простейшее линейное ОДУ:-У(О - -a-y(t)dtПри небольшом значении коэффициента а вычислительный блок Given-Odesolve (реализующий метод Рунге-Кутта 4-го порядка с фиксированным шагом) без проблем решает это уравнение (рис. 14.21). Например:а:=3Given-У(О = -a-y(t)dtУ(0) = 1f:=Odesolve(t,3,10)flt) 0.5 -Рис. 14.21. Верное решение ОДУПри увеличении же значения коэффициента а до 30 в решении появляется неустойчивость (рис. 14.22).123tРис. 14.22. Неустойчивость численного решения ОДУПриведенный пример показывает, что одни и те же уравнение или система уравнениймогут быть жесткими или нет в зависимости от значения (точнее, соотношения) входящих в них параметров.
Как правило, чем более различаются коэффициенты при разных членах уравнений системы ОДУ, тем более сложно решить ее с помощью обычных численных методов. Это связано с тем, что чем жестче система или уравнение, тембольше требуется шагов для устойчивой работы численного алгоритма. Так, например,14.2. Численное решение ОДУ в форме задачи Коши * 4 7 1чтобы получить верный результат для рассматриваемого уравнения при а=30, следуетразбить промежуток не на 10, а на 50 отрезков.
Однако в случае многих жестких систем невозможно прийти к верному решению обычным уменьшением длины шага, поскольку такой подход потребует столь огромного количества оборотов цикла, что задача будет не решаема в вычислительном плане.Разрешена же данная проблема была с помощью так называемых неявных абсолютноустойчивых разностных методов.
Особенностью этих методов является то, что ониустойчивы при любой величине шага, поэтому последний может определяться толькоисходя из соображений точности. В Mathcad жесткие системы ОДУ можно решать двумя способами: с помощью вычислительного блока Given-Odesolve, выбрав в контекстном меню функции Odesolve неявный алгоритм Stiff, либо задействовав встроенныефункции, реализующие наиболее популярные неявные методы.• Stiffb(yO,tO,tl,M,F,J). Метод Булирша-Штера для жестких систем ОДУ;• Stiffr(yO,tO,tl,M,F,3). Метод Розенброка для решения жестких систем ОДУ;LJ Radau(yO, tO, tl, M, F).
Метод RADAU5 (в терминологии справочной системы Mathcad) для жестких систем ОДУ.Данные функции требуют задания следующих параметров:• уО — вектор начальных условий. Задается точно так же, как при решении простыхсистем ОДУ;• tO — начальная точка для переменной;• tl — правая граница для интервала изменения переменной;• М — количество шагов, которое должен использовать численный метод;• F — векторная функция, определяющая систему ОДУ. Задается точно по таким жеправилам, как и в случае использования функций решения нежестких систем;• J — матричная функция размерности Nr(N+l), где N — количество уравнений в системе.
В первом столбце J должны находиться функции производных соответствующих строк F по независимой переменной. В остальных N столбцах должен бытьрасположен якобиан (матрица частных производных всех соответствующих уравнениям функций по всем входящим в них переменным) решаемой системы дифференциальных уравнений. Аналогично F, J должна быть задана как функция двухпеременных.Результатом работы перечисленных функций, как и всех остальных встроенных функций, предназначенных для решения систем ОДУ, является матрица, в первом столбцекоторой содержатся значения переменной, а в остальных — соответствующие им величины искомых функций.Обратите внимание, что функция Radau не требует задания якобиана, однако приее использовании необходимо повысить точность вычислений, уменьшив значениевстроенной константы TOL.
Точность расчета с помощью Stiffr и Stiffb также определяется TOL.Рассмотрим решение жесткой системы ОДУ на примере задачи химической кинетики. Пусть исходное вещество S под действием фермента Е превращается в конечныйпродукт Р. Фермент Е, действуя на вещество S, образует с ним промежуточный комплекс ES, после распада которого происходит образование продукта реакции Р и регенерация фермента Е. Комплекс ES не является абсолютно устойчивым и может разрушаться с образованием исходных реагентов.
Каждая из перечисленных стадий реакции4 7 2 • Глава 14. Дифференциальные уравненияпротекает с определенной скоростью, которая характеризуется константой скорости kn.Схематично все процессы можно представить следующим образом:E +S =к2E S - — Р +ЕПод скоростью реакции понимается изменение (уменьшение) концентрации реагирующих веществ в единицу времени. Чем больше kn, тем соответственно больше и скорость реакции.
Кроме того, нетрудно догадаться, что реакция будет протекать тем быстрее, чем больше концентрация исходных компонентов. Учитывая эти факты, запишемуравнения, описывающие изменение концентрации веществ в реакционной системес течением времени.Концентрация исходного реагента S уменьшается но мере образования промежуточного комплекса ES и восполняется за счет его частичного распада:d[S]= -kl-[E]-[S]+k2-[ES]dtКонцентрация фермента Е уменьшается по мере образования промежуточного комплекса ES.
Противоположный процесс наблюдается при регенерации фермента Е в результате образования продукта Р, а также при разрушении комплекса ES:4Ш] =-kl-[E]-[S] +(k2 +k3)-[ES]dtКонцентрация промежуточного комплекса ES увеличивается за счет связывания фермента Е с веществом S. При образовании конечного продукта Р и при разрушении комплекса концентрация ES в реакционной системе падает.- (k2 + k3)-[ES]= kldtКонечный продукт Р накапливается в системе за счет распада комплекса ES:кdtИтак, мы получили систему уравнений, решение которой представлено в примере 14.20.Чтобы получить задачу Коши, в систему необходимо добавить начальные условия —исходные концентрации фермента Е и вещества S.
Разумеется, в момент начала реакции концентрации промежуточного комплекса ES и продукта Р будут равны нулю.Изучая пример, обратите особое внимание на особенности задания параметра J.Пример 14.20. Решение жесткой системы ОДУ с помощью встроенныхфункцийkl := 10к2:=10'кЗ:= 10к2у 0У0:=0.1оF(t,y) :=к1- У 1 у 2 -(к2+кЗ)-у 014.2. Численное решение ОДУ в форме задачи КошиdtОdxLdtd v-Цdxd0 - Fvdx0 ^dx.— F x,dxdxу2dxx,7dx_At' 3 dx7x,7.4'Fy2 — F x,dxVdx•> Avx,3—Fl27XЕ1x,dxL V-x,yx,dxoУ13dxX^LCO100VоЛx7dt1VdxJ(t,y) :=F(t y)x,Vdx«• 4 7 3-1000000-У 2-1000000-yjVA/_0-ioooooo-y2 -юооооо-yj оJ(t,y)ЮОООООУ2юооооо-yj о1000100•-1000000У 2-1000000-yj-Ю00000-У 2 - 1 0 0 0 0 0 0 - ^J(t,y) :=ЮООООО-У211000si :=Stifib(yO,0,1000,5000, F,J)0юооооо-yj о00s2 := StiffifyO, 0,1000,5000, F, J)4 7 4 • Глава 14.
Дифференциальные уравненияПолученное нами решение абсолютно адекватно физически (рис. 14.23): с течениемвремени концентрация исходного вещества S падает (sl<l>) по мере накопления в реакционной системе конечного продукта Р (sl<4>).
На начальном этапе, когда количество Е и S в системе велико, происходит их интенсивное взаимодействие, вследствиечего концентрация комплекса ES резко возрастает. Затем за счет образования продукта Р и частичного разрушения комплекса концентрация ES постепенно снижается(sl<3>). Процесс превращения S в Р протекает с участием фермента Е с последующейего регенерацией, поэтому концентрация Е в ходе реакции практически не меняется100005001000Рис. 14.23. Решение жесткой системы дифференциальных уравненийОбратите внимание, каким интересным и необычным способом был подсчитан намиякобиан.
Мы использовали ту особенность решаемой системы, что в ней нет никакихспециальных функций, и все искомые элементы входят в нее в первой степени. А этоозначает, что при вычислении производной по определенной неизвестной, последняябудет исключена из уравнения. Поэтому, не сделав никакой фактической ошибки, мыпросто заменили соответствующие элементы вектора у на переменную х. Необходимость этой операции связана с тем, что оператор дифференцирования Mathcad не может вычислять производной по матричному элементу.Безусловно, вы можете сказать, что приведенный в рассматриваемом примере способзадания якобиана чрезмерно громоздок и гораздо легче было бы просто подсчитать егоэлементы в голове. Конечно, в случае таких простых для вычисления производныхфункций, какими являются правые части уравнений химической кинетики, следуетпоступить именно так.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
