Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3791 / 1086 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Специальности: Программист, Математик

Лекция 10: Численные методы решения жестких систем обыкновенных дифференциальных уравнений

9.3. Решение линейных ЖС ОДУ и вычисление матричной экспоненты

Рассмотрим один из наиболее простых численных методов решения жестких линейных систем ОДУ [9.10], основанный на представлении решения в явном виде

u_{n + 1} = e^{B\tau }u_{n}

для линейных систем жестких ОДУ вида

$ \dot {u} = Bu, {u}(0) = u_0 . $

Его численная реализация связана с вычислением матричной экспоненты. О свойствах матричной экспоненты (и в целом функций от матриц) можно прочитать в [9.11]. Использование для этого разложения в ряд Тейлора

$  e^{A{\tau}} = {E} + \sum\limits_{k = 1}^{N}{\frac{\tau^k}{k!}{B}^{k}}  $

представляется непригодным, так как простые оценки показывают, что по вычислительным затратам этот алгоритм сопоставим с методом Эйлера с шагом \tau || B ||  < 1:

$  \frac{{{\tau}^{k}}}{{k!}}\|{B}^{k}\|   \approx  \left({\frac{e\tau\|{B}\|}{k}}\right)^{k},  \mbox{ так как }  k!  \sim  \left({\frac{k}{e}}\right)^{k} . $

Следовательно, для того, чтобы k - й член разложения ряда Тейлора был, по крайней мере, порядка O(1), необходимо выполнение условия

$  \frac{{e\tau\|{B}\|}}{k} \sim  O(1), \mbox{ или } \tau\|{B}\|  \sim  
O(1), $

что соответствует условию численного интегрирования с шагом \tau || B || < 1\}. Количество членов ряда при этом N  \sim  \tau\|{B}\| \gg 1, что неприемлемо для решения.

Представим экспоненту в следующем эквивалентном виде:

$  e^{{B}{\tau}} = (e^{{\tau}A/2^p})^{2^p} \approx {\left[{E + \frac{{\tau}}{2^p}B + \ldots + \frac{1}{k!}{\left({\frac{{\tau}}{2^p}}\right)}^k B^k}\right]}^{2^p}.$

При этом значение параметра p выбирают таким, чтобы

$  \frac{{\tau\|{B}\|}}{{2^{p}}} \ll 1 $

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

$  \frac{{\tau\|{B}\|}}{{k \cdot 2^{p}}} \sim  O(1)  \mbox{ и }  N \sim  2^{- p} \tau\|{B}\|, $

что вполне приемлемо при соответствующем выборе параметра p.

В этом алгоритме сначала вычисляют матрицу

$  D = E + \sum\limits_{k = 1}^N{\left({\frac{\tau}{2^p}}\right)^k B^k}, 
  $

затем D^{2^p} путем последовательных перемножений. При таком способе вычисления матрицы также имеется опасность, связанная с тем, что при некоторых p, \lambda _{0}, \Lambda _{0} слагаемые, соответствующие мягкой части спектра, окажутся много меньше единицы (и слагаемых, соответствующих жесткой части). В этом случае предпочтение отдается представлению

$  e^{A{\tau}} = \left[{\exp  \left({- \frac{{\tau}}{2^p}{B}}\right)^{- 1}}\right]^{2^{p}}, $

а последовательность вычислений имеет вид

$  D = E + \frac{{\tau}}{2^p}B; C = D^{- 1};  C^{\prime} = {(C^2 )}^p .  $

При этом на мягкой части спектра, поскольку {\tau}|\lambda_j | \ll 1, имеем

$  \left({1 - \frac{{\tau}}{{2^{p}}} |\lambda_j|}\right)^{- 2^{p}}
 \approx  e^{\lambda {\tau}} ;  $

на жесткой части, поскольку теперь p небольшие и {\tau}|{\Lambda_j}|/2^p \gg 1, получим оценку

$  \left|{\left({\frac{1}{1 - {\tau}|\Lambda_j | \cdot 2^p}}\right)^{2^p}}\right| \ll 1,   $

что уже приемлемо для вычислений.