Solve_ODE (ЭВМ для спецгруппы), страница 2

Преимущество неявной схемы (10) по сравнению с явным методом Эйлера (4) при интегрировании дифференциального уравнения в области убывающего решения (<0) заключается в том, что A‑устойчивость позволяет выбирать шаг, руководствуясь исключительно условием, чтобы локальная ошибка не превышала допустимой погрешности результата. Для явной схемы, где устойчивость имеет место лишь в ограниченной области значений h, может случиться, что шаг будет лимитирован не локальной точностью, а эффектом многократного возрастания малых ошибок, которые сами по себе значительно ниже уровня допустимой погрешности.

2.3. Методы высших порядков

Рассмотренные в предыдущих разделах явная и неявная схемы метода Эйлера соответствуют линейной аппроксимации решения в пределах одного шага. Сохраняя большее число членов в разложении (2), можно построить методы интегрирования, основанные на локальной аппроксимации решений многочленами 2‑й, 3‑й и более высоких степеней. Степень аппроксимирущего многочлена называют порядком метода. Так, метод Эйлера является методом первого порядка.

Чем выше порядок метода, тем меньше ошибка усечения, а следовательно, тем точнее аппроксимируется решение на одном шаге. В общем случае для метода порядка p локальная ошибка соответствует O(h^p+1). Если уровень допустимой погрешности фиксирован, то более высокий порядок позволит интегрировать с бо^льшим шагом, а значит, сократить время решения задачи. Однако, это справедливо лишь при условии устойчивости метода.

Анализ устойчивости схемы численного интегрирования в общем случае представляет собой весьма сложную проблему. Поэтому мы приведем здесь лишь формулировки двух теорем, которые были доказаны Далквистом [1] в 1963 г. [Речь идет только о линейных многошаговых методах]:

1) Ни один явный метод не может быть A‑устойчивым.

2) Неявный метод может обладать A‑устойчивостью только при порядке p2.

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

Следует особо остановиться на формулах численного интегрирования, применяемых в методах высоких порядков. Если исходить непосредственно из разложения в ряд Тейлора (2), то для аппроксимации решения многочленом степени p потребуются производные до порядка p включительно. Однако дифференциальное уравнение вида (1) обеспечивает непосредственное вычисление только первой производной. По этой причине многочлен строят на основании уже вычисленных значений решения и его первой производной в нескольких предыдущих точках. Такие схемы носят название многошаговых или многоточечных. В качестве примера рассмотрим один из наиболее точных методов прогноза и коррекции 4‑го порядка — метод Милна. Для прогноза используется явная 4‑точечная формула

.

Затем проводится коррекция решения однократным применением неявной формулы

.

Локальная ошибка усечения после коррекции оценивается величиной , где  – некоторое значение независимой переменной x в промежутке между x_n3 и x_n+1.

При использовании многоточечной схемы возникает проблема запуска вычислений. Например, чтобы получить решение в очередной точке с помощью метода Милна, необходимо иметь значения x и y в 4‑х предыдущих точках. Постановка задачи Коши дает единственную начальную точку. Поэтому до применения метода Милна необходимо построить решение в нескольких первых точках с помощью какой-либо одноточечной схемы (например, методом Эйлера). Многоточечные схемы типа метода Милна относят к несамостартующим, тогда как метод Эйлера является самостартующим. Существуют также самостартующие методы высоких порядков (например, метод Рунге-Кутта 4‑го порядка).

2.4. Решение систем ОДУ и проблема жесткости

Перейдем теперь к рассмотрению систем дифференциальных уравнений. Пусть, например, дана система двух уравнений

с начальными условиями

Общий подход к построению схемы численного интегрирования системы ОДУ ничем не отличается от случая одного уравнения. Точно так же, как раньше, используя заданные начальные условия для точки x₀, вычисляют компоненты решения y₁ и y₂ в точке x₁=x₀+h, пользуясь, например, той же формулой Эйлера (4) для каждого из уравнений в отдельности. Однако при анализе ошибок и устойчивости решения появляются некоторые особенности.

Нужно учитывать, что функции y₁(x) и y₂(x) могут вести себя по‑разному, так что локальная и глобальная ошибки у каждой компоненты решения будут своими, причем локальная ошибка зависит только от свойств данной компоненты (определяется видом функции F_i), а в глобальную ошибку вносят вклад ошибки других компонент, присутствующих среди аргументов F_i. Поскольку шаг численного интегрирования один и тот же для всех уравнений системы, то его величина, очевидно, должна лимитироваться наименее устойчивой компонентой (или компонентой, имеющей наибольшую локальную ошибку, если метод вполне устойчив).

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

(11)

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

Известно (см., например, [2]), что общее решение системы (11) можно записать в виде

(12)

где ₁ и ₂ – собственные значения матрицы A, составленной из коэффициентов a_ij системы (11), а (c₁₁, c₂₁) и (c₁₂, c₂₂) – собственные векторы, соответствующие ₁ и ₂. Поскольку собственные векторы матрицы определяются с точностью до постоянного множителя, то для однозначного выбора решения необходимо учесть начальные условия. Если система (11) интегрируется методом Эйлера, то поведение глобальных ошибок будет зависеть от величин ₁, ₂ и шага h, причем характер этой зависимости качественно соответствует тому, что наблюдается в случае одного уравнения (выражения (7) и (10)).

Что произойдет, если величины ₁ и ₂ будут различаться на несколько порядков? Пусть, например, ₁=1, ₂=1000. Решение (12) в этом случае содержит две составляющие: медленно ( ) и быстро меняющуюся ( ). Для явной схемы Эйлера условие устойчивости решения приводит к тому, что шаг интегрирования лимитируется «быстрой» составляющей (h<0.002). В то же время уже при смещении от точки x₀ на расстояние x=0.01 «быстрая» составляющая уменьшится по сравнению с начальным уровнем более чем в 2·10⁴ раз, т. е. практически исчезнет, тогда как «медленная» составляющая к этому моменту успеет измениться всего на 1%. Таким образом, за исключением очень узкого начального участка вид решения целиком определяется «медленной» составляющей, т. е. именно эта составляющая характеризует масштаб по оси x для процесса, описываемого системой уравнений (11). Например, интервал, на котором функции y₁(x) и y₂(x) изменятся в несколько раз, соответствует x>1. Чтобы получить решение на таком интервале, потребуется сделать не менее 1000 шагов. Если мы попытаемся увеличить размер шага, то присутствие в уравнениях членов, отвечающих ₂, немедленно приведет к неустойчивости. Хотя значение почти везде неотличимо от нуля, самые незначительные ошибки в этой составляющей, появившиеся, например, в результате округлений, лавинообразно усиливаются и через несколько шагов делают дальнейшее решение бессмысленным. Внешне это выглядит как колебания вычисляемых значений y_i со все возрастающей амплитудой, пока не произойдет переполнение (превышение предельной величины чисел, представимых на ЭВМ). В то же самое время локальная ошибка усечения для метода Эйлера определяется величиной второй производной (главный член отбрасываемой части разложения решения в ряд Тейлора), которая на всем интервале интегрирования за исключением короткого начального участка соответствует «медленной» составляющей. Иными словами, если бы не было неустойчивости, то шаг интегрирования ограничивался бы только допустимой локальной погрешностью решения, т. е. почти везде определялся бы не «быстрой», а «медленной» составляющей.

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

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

, (13)

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

Именно жесткость создает основные трудности при численном решении дифференциальных уравнений химической кинетики [3]. Форма кинетических уравнений такова, что матричные элементы якобиана (13) содержат либо константы скоростей элементарных стадий, либо произведения константы скорости на концентрации одного или двух реагентов (в случае би- и тримолекулярных стадий, соответственно), причем диагональные элементы отрицательны либо равны нулю. В первом приближении можно считать, что диапазон собственных значений матрицы Якоби J примерно соответствует диапазону значений констант скоростей [приведенных к первому порядку]. Во многих реальных химических процессах с участием активных промежуточных частиц (атомов, радикалов) различие констант на 10-15 порядков является вполне обычным, а это означает, что при интегрировании до времен, характерных для медленных реакций (определяющих продолжительность процесса в целом), потребуется 10¹⁰-10¹⁵ шагов, если применять стандартные численные методы. Такой объем вычислений недоступен даже для самых быстродействующих современных ЭВМ, не говоря уже о накоплении ошибок округления из‑за огромного количества арифметических операций.

В качестве примера можно привести результаты, полученные в работе [4]. Исследовалась модель химической стадии радиолиза водного раствора метанола, которая включала 37 элементарных стадий с константами скорости в диапазоне от 10³ до 10¹¹ М^1с^1. При попытке интегрирования соответствующей системы ОДУ методом Рунге-Кутта оказалось, что максимальная величина шага, при которой еще сохраняется устойчивость решения, не превышает 10^9 с. Для получения решения в интервале 0t4·10^5 с время счета на ЭВМ с быстродействием порядка 1 Мфлоп/с (БЭСМ‑6) составило около 5 минут. Однако практический интерес представляло моделирование хода процесса в течение минут и даже часов, что оказалось просто неосуществимым, так как, например, получение решения в интервале до 100 секунд потребовало бы времени счета порядка 10⁷ минут или 20 лет!

2.5. Методы интегрирования жестких уравнений

Наиболее подходящим для интегрирования жестких систем мог бы быть метод, обладающий A‑устойчивостью. Однако, согласно теореме Далквиста [1], порядок такого метода не может быть выше 2, а в этом случае шаг все равно приходится делать весьма малым из‑за больших ошибок усечения. Хотя до 1970 года устойчивые неявные схемы 1‑го порядка и применялись на практике, их эффективность была недостаточной.