Введение в вычислительную математику
[+]
Московский физико-технический институт
Опубликован: 25.10.2007 | Доступ: свободный | Студентов: 3976 / 1246 | Оценка: 4.50 / 4.33 | Длительность: 24:00:00
ISBN: 978-5-9556-0065-9
Темы: Алгоритмы и дискретные структуры, Математика
Специальности: Программист, Математик
Теги: beta, алгебра, алгоритмы, арифметика, вычисления, дифференциальные уравнения, курсы, линейная алгебра, матрица Якоби, метод трапеций, метод Эйлера, норма вектора, поиск, программирование, разделенные разности, разложение в ряд, системы нелинейных уравнений, спектральные задачи, таблицы бутчера, численное решение уравнений, элементы

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

A
|
версия для печати
< Лекция 9 || Лекция 10: 12345678910 || Лекция 11 >

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, $

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

Дальше >>
< Лекция 9 || Лекция 10: 12345678910 || Лекция 11 >

Вопросы и ответы

вопросов: 1
Sanjar Ibrohimov
Sanjar Ibrohimov
ответить