Рассмотрим линейную систему ОДУ первого порядка
$$$ \frac{du}{dt} = Au + f, u \in R^n, t \in [0, L] $$$с краевыми условиями
Ru(0) + Su(L) = q,
где u, f, q — n - мерные векторы, A(t), R(t), S(t) — матрицы размера n x n.
Для приближенного решения задачи введем расчетную сетку $$\{t_n\}_{n = 0}^{N}$$ и за приближенное решение примем сеточную функцию $$\{u_n\}_{n = 0}^{N}.$$ Рассмотрим методы построения приближенного решения.
Метод построения фундаментальных решений аналогичен известному по курсу дифференциальных уравнений способу построения общего решения системы линейных уравнений первого порядка. Решение представляется в виде
$$$ u(t) = \bar{u}(t) + \sum\limits_{k = 1}^n {\alpha_k u^k}. $$$Здесь uk(t) есть полная фундаментальная система решений
однородной задачи
с начальными данными, например,
uk(0) ={0, ..., 0, 1, 0, ..., 0}T,
где единица стоит на k месте, т.е. в качестве начальных данных
используются векторы uk(0) = Ek. Важно, чтобы решения однородной задачи составляли систему линейно независимых функций. Каждая такая функция ищется численно как решение соответствующей задачи Коши, используя методы, описанные в лекциях 8 и 9.
Пусть $$$ \bar{u}(t) $$$ — частное решение неоднородной системы
$$$ \frac{d\bar{u}}{dt} = A\bar{u} + f $$$с нулевыми начальными условиями $$$ \bar{u}(0) = 0 $.$$ Тогда неопределенные коэффициенты $$\alpha_k$$ находятся из краевых условий
$$\begin{gather*} Ru(0) + Su(L) = q, \mbox{ или }\\ R[\sum\limits_{k = 1}^n{\alpha_k{u}^k (0)}] + S\bar {u}(L) + S[\sum\limits_{k = 1}^n{\alpha_k{u}^k(L)}] = q. \end{gather*}$$Последнее соотношение представляет собой СЛАУ относительно коэффициентов $$\alpha_{k}$$ размерности n.
Полную фундаментальную систему решений однородной задачи можно получить,
используя, например, схему второго порядка точности (
Частное решение получаем аналогично:
$$\begin{gather*} \frac{\bar{u}_{n + 1} - \bar{u}_n}{\tau } = A_{n + 1/2} \frac{\bar{u}_{n + 1} + \bar{u}^{\prime\prime}_n}{2} + f_{n + 1/2}, \\ f_{n + 1/2} = F\left({t_n + \frac{\tau }{2}}\right), \\ \bar{u}_0 (0) = 0, \mbox{или}\\ \bar{u}_{n + 1} = \left({E - \frac{\tau}{2}A}\right)^{- 1}\left[{\left({E + \frac{\tau}{2}A}\right)u_n + \tau f_{n + 1/2}}\right], \bar{u}_0 = 0. \end{gather*}$$Приведем пример, когда решение задачи методом построения фундаментальных решений не проходит. Рассмотрим систему уравнений
$$\begin{gather*} \frac{du}{dt} = {av} + f, \frac{dv}{dt} = bu + g, \\ u(0) = U_0, v(1) = 0 , \end{gather*}$$Подобные системы ОДУ моделируют, например, процессы прохождения излучения или потоков элементарных частиц (гамма - излучение, потоки нейтронов) через разные среды (атмосфера, защита ядерных реакторов) в приближении оптически толстого слоя. Коэффициенты $$a, b \sim 50$$ характерны для защиты реакторов. Найдем общее решение этой задачи с помощью метода фундаментальных систем в виде линейной комбинации двух решений однородных ОДУ:
$$\left( \begin{array}{l} u \\ v \\ \end{array} \right) = \left( \begin{array}{l} {\bar {u}} \\ {\bar {v}} \\ \end{array} \right) + \alpha_1 \left( \begin{array}{l} u^1 \\ v^1 \\ \end{array} \right) + \alpha_2 \left( \begin{array}{l} u^2 \\ v^2 \\ \end{array} \right),$$где (u1, v1) и (u2, v2) — решения двух однородных систем (в дальнейшем для простоты будем полагать, что $$$ \bar {u} = \bar {v} = 0 $.$$) Тогда
а коэффициенты $$\alpha_1$$ и $$\alpha_2$$ находятся из краевых условий.
Общее решение такой системы, как известно, представляет собой сумму двух экспонент
$$\left( \begin{array}{l} u^1 \\ v^1 \\ \end{array} \right) = C_1 \left( \begin{array}{l} a_1 \\ b_1 \\ \end{array} \right)e^{\lambda_1 t} + C_2 \left( \begin{array}{l} a_2 \\ b_2 \\ \end{array} \right)e^{\lambda_2 t}$$(аналогичным образом решение представляется и для u2, v2 ), где ai, bi, i = 1, 2 находятся из решения задачи, Ci; i = 1, 2 — произвольные постоянные; $$\lambda _{1}$$ и $$\lambda _{2}$$ — корни характеристического уравнения$$\left| \begin{array}{cc}
{- \lambda} {a} \\
{b} {- \lambda} \\
\end{array} \right| = 0.$$
Решая U0e - 50t (в случае задачи с защитой ректора U0 — падающий поток нейтронов; задача защиты — значительное его ослабление, приблизительно в e50 раз).
Слагаемые с быстрорастущими экспонентами должны взаимно уничтожиться.
Получение же численного решения - весьма трудная задача, поскольку численное решение имеет большую и быстро возрастающую погрешность. Пусть $$u^{M} = u(1 + \varepsilon ) \sim e^{50}t(1 + 10^{ - 12}),$$ т.е. начальная погрешность имеет порядок 10 - 12. Она возрастает при вычислениях примерно в e50t раз — даже при умеренных t это очень большое число.
Рассмотрим метод
где $$\alpha$$ и $$\beta$$ — пока неизвестные функции (прогоночные коэффициенты), для которых необходимо получить дифференциальные уравнения.
Продифференцируем это соотношение
$$$ \dot {u} = \dot {\alpha} v + \alpha \dot {v} + \dot {\beta} $$$и подставим в него уравнения системы $$$ \dot {u} = {av} + f, \dot {v} = {bu} + g $.$$ В результате получим, что $$$ {av} + f = \dot {\alpha} v + \alpha ({bu} + g) + \dot {\beta} $.$$ Подставим в полученное соотношение уравнение $$u = \alpha v + \beta .$$ Тогда $${av} + f = \dot {\alpha} v + b\alpha ^2 v + b\beta \alpha + \alpha g + \dot {\beta} v + \beta .$$
После приведения подобных членов имеем равенство
$$$ v(\dot {\alpha} + b\alpha ^2 - a) + (\dot {\beta} + b\beta \alpha + \alpha g - f) = 0. $$$Приравнивая к нулю коэффициенты при v и единице, получим два
дифференциальных уравнения для прогоночных коэффициентов:
Дополним их начальными условиями. Левое краевое условие вида u (0) =
U0 запишем в виде прогоночного соотношения $$u(0) = \alpha (0) v(0) + \beta(0),$$ полагая $$\alpha (0) = 0, \beta(0) = U_0.$$ Таким образом, получаем начальные данные для двух задач Коши для $$\alpha (t)$$ и $$\beta (t),$$ которые могут быть решены численно.
Теперь разрешим правое краевое условие. На правой границе отрезка интегрирования имеем условие v(1) = 0 и прогоночное соотношение при t = 1: $$u (1) = \alpha (1)v(1) + \beta (1),$$ откуда получаем $$u(1) =
\beta (1).$$
Далее воспользуемся уравнением $$$ \dot {v} = {bu} + g $,$$ подставив в него прогоночное соотношение $$u = \alpha v + \beta,$$ получим дифференциальное уравнение для v:
Интегрируя эту задачу справа налево, попутно определяем u(t):
Метод
Рассмотрим теперь общую постановку
где $$u, q_i , D_i, f \in R^n, \left\|{D_i}\right\| \simeq O(1), i = 1, \ldots,.$$
Здесь $$\mathbf{A}$$ — постоянная матрица размером N
x N (дальнейшие рассуждения будут справедливы и для систем уравнений с переменными коэффициентами). Определители систем линейных алгебраических уравнений, которыми являются краевые условия на обоих концах интервала интегрирования, полагаются отличными от нуля.
Определение. Рассматриваемая
Отношение $${{\Lambda_0 }/{\lambda_0 }} \gg 1$$ является параметром, характеризующим жесткость системы. В дальнейшем будем полагать $$\Lambda _0 (b - a) \gg 1,$$ $$\lambda_0 (b - a) \simeq O(1).$$
Общее решение такой системы имеет вид
$${u}(t) = \sum\limits_{i = 1}^{N_1 }{\gamma_i^1 e^{\Lambda_i^1t}{{\Omega }}_i^1 + \sum\limits_{i = N_1 + 1}^{N_2 }{\gamma_i^2 e^{\Lambda_i^2 t}{{\Omega }}_i^2 + \sum\limits_{i = N_2 + 1}^{N}{\gamma_i^3 {e}^{\lambda_it}} }}{{ \omega }}_i,$$где $$\Omega_i^1, \Omega_i^2, \omega_i$$ есть собственные векторы матрицы A, соответствующие трем частям спектра. Понятна качественная структура этого решения, содержащего как левый, так и правый пограничные слои. Будем полагать, что количество собственных значений в каждой из трех частей спектра не изменяется. Особенность || f || и || q || - нормы правых частей в системе ОДУ и краевых условиях, соответственно. Для численного решения задачи
можно использовать такую же схему второго порядка, как и ранее. Рассматриваем класс вычислительно корректных задач, для которых $$C = O(1) \ll \exp(\Lambda_0(b - a)).$$
В дальнейшем будем полагать величину $$||A|| (b - a) \approx \Lambda _{0}(b - a)$$ большой, а величину $$exp(||A|| (b - a)) \approx exp(\Lambda _{0}(b - a))$$ — очень большой. В приложениях такие задачи встречаются наиболее часто. Важно отметить и то, что не все возможные постановки задач для жесткой системы приводят к вычислительно корректным алгоритмам. Показывается, что необходимыми (и почти достаточными) условиями корректности являются следующие неравенства: $$k \ge N_1, (N - k) \ge N_2 - N_1,$$ т.е. число краевых условий на левом конце отрезка интегрирования не должно быть меньше быстро убывающих вправо решений, на правом конце — не меньше быстро убывающих влево решений. В противном случае
Проблемы, которые возникают при численном решении eLt, как известно, приводит к потере точности и накоплению вычислительных ошибок. Вторая проблема состоит в следующем. Для вычисления коэффициентов $$\alpha_i,$$ входящих в общее решение неоднородной системы (а оно состоит из суммы частного решения неоднородной системы и общего решения однородной $$u (t) = {\bar{u}}_0 + \sum\limits_{i = 1}^{N}{\alpha_i}u_i$$ ), приходится решать плохо обусловленную СЛАУ $$(D_i, \bar{u}(b) + \sum\limits_i {\alpha_i}u_i(b)) = q_i, i = k + 1, \ldots,.$$
где коэффициенты g, h, s, вообще
говоря, являются функциями независимого переменного t.
Для разностной аппроксимации рассматриваемой краевой задачи введем равномерную разностную сетку $$\{t_n\}_0^{N}, t_n = n\tau , \tau = (b - a)/N$$ и определим на этой сетке сеточную функцию $$\{u_n\}_0^{N}.$$
Коэффициент g, вообще говоря, может не иметь первой производной. Такая задача может возникнуть, например, в случае расчета установившегося распределения температуры в задаче стационарной теплопроводности с контактным разрывом. Представим разностную задачу в виде
Контактный разрыв при этом помещается в узел с номером n.
В этом случае фактически аппроксимируется тепловой поток через границы ячейки разностной сетки, для уравнения теплопроводности получается консервативная разностная схема — подробнее смотри в лекциях, посвященных
Выписанные выше соотношения определяют простейшую разностную схему. Под разностной схемой здесь и ниже понимается совокупность разностных уравнений для определения значений сеточной функции внутри расчетной области, дополненная соответствующими начальными и граничными условиями для этой сеточной функции.
Для определения значений сеточной функции получается СЛАУ с трехдиагональной матрицей
- b0u0 + c0u1 = d0, anun - 1 - bnun + cnun - 1 = dn , n = 1, ..., N - 1, anuN - 1 - bnun = dn,
где $$\begin{gather*} a_n = \frac{{g_{n - 1/2}}}{{\tau ^2 }} - \frac{{h_n}}{{2\tau }}, c_n = \frac{{g_{n + 1/2}}}{{\tau ^2 }} + \frac{{h_n}}{{2\tau }}, \\ b_n = a_n + c_n - h_n , d_n = f_n , b_0 = \frac{A}{\tau } - \frac{B}{2}, c_0 = \frac{A}{\tau } + \frac{B}{2}, d_0 = D \end{gather*}$$
Эта СЛАУ представима в каноническом виде A'u = D,
где A' — матрица
u, d есть векторы - столбцы
D = (d0, d1, ..., dn)T.
Трехдиагональные матрицы часто возникают при численном решении краевых задач как для обыкновенных дифференциальных уравнений, так и для уравнений в частных производных. Ранее матрица подобной структуры встречалась при построении сплайна Шонберга (лекция 6). Характерная особенность таких матриц заключается в том, что при большой размерности матрица имеет ленточную структуру — все элементы вне ленты (n - го порядка требует O(N3) арифметических действий и O(N2) ячеек памяти. В численных методах большую роль играют экономичные алгоритмы, в которых количество арифметических операций пропорционально количеству неизвестных — O(N).
К таким алгоритмам относится метод трехточечной разностной
Решение ищется в виде прогоночного соотношения
un - 1 = pnun + qn, n = 1, ..., N,
где pn и qn — прогоночные коэффициенты, подлежащие определению.
Левое краевое условие также записывается в виде прогоночного соотношения
$$$ u_0 = \frac{c_0 }{d_0 }u_1 - \frac{d_0 }{b_0}, $$$где p1 = c0/ b0, q1 = - d0/ b0 (заметим, что если A > 0 и B < 0, то b0 > c0, 0 < p < 1 ).
Получим рекуррентные формулы, позволяющие последовательно вычислить p2, q2, p3, q3 и т.д. вплоть до pn, qn.
Подставив равенство u n - 1 = pnun + qn в уравнение anun - 1 - bnun + cnu n + 1 = dn, получим an(pnun + qn) - bnun + cnun + 1 = dn, или
Сравнивая эту запись со стандартным видом прогоночного соотношения
un = pn + 1un + 1 + qn + 1,
видим, что для прогоночных коэффициентов должны выполняться равенства
$$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n}, \quad q_{n + 1} = \frac{a_n q_n - d_n}{b_n - a_n p_n}. $$$Эти формулы определяют прямой ход прогонки.
Из краевого условия на правом конце отрезка интегрирования anuN - 1 -
bNuN = dN и прогоночного соотношения uN - 1 = pNuN + qN находим величину uN.
Далее последовательно вычисляются остальные неизвестные un n = N - 1, ..., 1 un - 1 = pnun + qn. Это — обратный ход алгоритма прогонки.
Теперь, после того как описан прогоночный алгоритм, исследуем его на устойчивость.
Для этого рассмотрим вычисление прогоночного коэффициента pn (т.е. этап прямого хода прогонки). В идеальной арифметике этот коэффициент равен $$$ p_{n + 1} = \frac{c_n}{b_n - a_np_n} $,$$ в конечноразрядной арифметике — $$p_{n + 1}^{M} = p_{n + 1} + \Delta_{n + 1} ,$$ где $$\Delta _{n + 1}$$ — погрешность, связанная с округлениями на всех предшествующих этапах вычислений. Полагая $$\Delta _{n + 1}$$ малой, исследуем изменение этой погрешности с ростом n. Для этого запишем соотношение между $$\Delta _{n}$$ и $$\Delta _{n + 1}$$ в виде
где $$\varepsilon _{n}$$ — погрешность вычислений правой части и машинного представления коэффициентов an, bn, cn. Следовательно, $$\Delta _{n + 1}$$ складывается из двух составляющих — локальной погрешности $$\varepsilon _{n}$$ и наследственной $$\Delta _{n + 1}.$$
Полагая $$\Delta_n \ll p_n$$ при $$\Delta _{n} > 0, p_{n} > 0$$ и опуская члены порядка $$O(\Delta ^{2}),$$ получим оценку
$$$ p_{n + 1} + \Delta_{n + 1} \approx \frac{c_n}{b_n - a_n p_n} + \frac{c_n a_n}{{(b_n - a_n p_n)}^2 }\Delta_n + \varepsilon_n . $$$Отсюда, учитывая, что $$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n} $,$$ получим требуемую оценку для эволюции погрешности
$$$ {\Delta_{n + 1} = \frac{a_n}{c_n}p_{n + 1}^2 \Delta_n + \varepsilon_n, n = 0, \ldots , N - 1.} $$$Докажем следующую теорему.
Теорема. Если выполнены условия диагонального преобладания $$\left|{b_n}\right| \ge \left|{a_n}\right| + \left|{c_n}\right|$$ и хотя бы для одной строки матрицы системы имеет место строгое диагональное преобладание $$\left({\left|{b_n}\right| > \left|{a_n}\right| + \left|{c_n}\right|}\right).$$ Пусть, кроме того, 0 < p1 < 1. Тогда алгоритм прогонки устойчив.
Доказательство.
Докажем утверждение теоремы для случая, когда во всех строках матрицы выполнено условие строгого диагонального преобладания. На случай нестрогих неравенств обобщение доказательства очевидно.
0 < p1 < 1. Для определенности положим an > 0, bn > 0, cn > 0.Тогда $$$ p_{n + 1} = \frac{c_n}{b_n - a_np_n} > \frac{c_n}{b_n - a_n} > 0 $.$$
Кроме того, $$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n} < \frac{c_n}{a_n + c_n - a_n p_n} = \frac{c_n}{c_n + a_n (1 - p_n)} \le 1 $,$$ откуда следует 0 < pn + 1 < 1.
Для этого вспомним выражения для коэффициентов линейной системы, полученной при разностной аппроксимации исходного уравнения второго порядка.
$$\begin{gather*} a_n = \frac{g_{n - 1/2}}{\tau ^2 } - \frac{h_n}{2\tau } \approx \frac{g_{n - 1/2}}{\tau ^2 }, g_{n - 1/2} = g(t_n - \tau /2), \\ c_n = \frac{g_{n + 1/2}}{\tau^2} + \frac{h_n}{2\tau} \approx \frac{g_{n + 1/2}}{\tau^2}, \\ \frac{a_n}{c_n} \approx \frac{g(t - \tau/2)}{s(t + \tau/2)} = 1 + c\tau. \end{gather*}$$Вернемся к выражению для эволюции погрешности (10.2). С учетом полученных оценок имеем $$$ \Delta_{n + 1} \le \frac{a_n}{c_n}\Delta_n + \varepsilon_n $.$$ Так как $$$ \frac{a_n}{c_n} \approx 1 + c\tau $,$$ то $$\Delta_{n + 1} \le (1 + c\tau )\Delta_n + \varepsilon_n,$$ где c — константа Липшица для функции g(t). Считаем, что на каждом шаге ошибка округления не превосходит предельного значения, т.е. $$0 \le \varepsilon_n \le \varepsilon.$$
Из цепочки неравенств
$$\begin{gather*} \Delta_2 \le (1 + c\tau )^2 \Delta_0 + \varepsilon \left[{1 + (1 + c\tau )}\right], \\ \Delta_3 \le (1 + c\tau )^3 \Delta_0 + \varepsilon \left[{1 + (1 + c\tau ) + (1 + c\tau )^2 }\right], \end{gather*}$$будет следовать оценка $$$ \Delta_n \le (1 + c\tau )^{n} \Delta_0 + \frac{{(1 + c\tau )^{n}}}{{c\tau }}\varepsilon $.$$ Используя известные из математического анализа неравенства, последнюю оценку можно записать в виде
$$$ \Delta_n \le e^{c\tau n} \Delta_0 + e^{c\tau n} \frac{\varepsilon }{{c\tau }}, c\tau n = ct . $$$При $$ct \sim O(1)$$ погрешности не накапливаются (для краевых задач это реальное условие), поскольку в расчетах $$\tau \gg \varepsilon,$$ где $$\varepsilon$$ — машинный эпсилон (подробнее в лекции 1).
Аналогичное утверждение доказывается и для второго прогоночного коэффициента $$$ q_{n + 1} = \frac{a_n q_n - d_n}{b_n - a_n p_n} $;$$ алгоритм устойчив при выполнении тех же условий.
Покажем устойчивость обратного хода прогонки.
При обратном ходе вычисления проводятся по формулам un = pn + 1un + 1 + qn + 1, откуда, учитывая, что $$u_n^{M} = u_n + \Delta_n$$ и $$u_{n + 1}^{M} = u_{n + 1} + \Delta_{n + 1},$$ получим $$\Delta _{n} = p_{n + 1}\Delta _{n + 1} + \varepsilon _{n},$$ где $$\Delta _{n}$$ — наследственная погрешность, $$\varepsilon _{n}$$ — погрешность округления на n шаге. Очевидно, что обратный ход прогонки устойчив при выполнении условия 0 < pn < 1, или 0 < p1 < 1.
При численном решении краевых задач для обыкновенных дифференциальных уравнений 4 - го порядка возникают СЛАУ с пятидиагональной матрицей вида:
c0u0 - d0u1 + e0u2 = f0, - b1 u0 + c1 u1 - d1 u2 + e1 u3 = f1, anun - 2 - bnun - 1 + cnun - dnun + 1 + enun + 2 = fn, n = 2, ... , N - 2, aN - 1 uN - 3 - bN - 1 uN - 2 + cN - 1 uN - 1 - dN - 1 un = fN - 1, an uN - 2 - bn uN - 1 + cn un = fn .
Алгоритм решения таких систем — пятиточечная прогонка, формулы которой
выводятся аналогично формулам для трехточечной прогонки. Приведем их окончательный вид. В прогоночном соотношении появятся три коэффициента (p, q, r):
un = pn + 1 un + 1 - qn + 1un + 2 + rn + 1, uN - 1 = pn un + rn, un = rN + 1.
Прогоночные коэффициенты находятся по формулам
$$p_{n + 1} = [d_{n} + q_{n}(a_{n} p_{n - 1} - b_{n})]/D_{n} , n = 2, 3, \dots , N - 1, \setminus \setminus \\ p_{1} = d_{0} /c_{0}, p_{2} = (d_{1} - q_{1} b_{1} )/D_{1}, \\ r_{n + 1} = [f_{n} - a_{n} r_{n - 1} - r_{n} (a_{n} p_{n - 1} - b_{n})]/D_{n}, n = 2, 3, \dots , N, \\ r_{1} = f_{0} /c_{0}, r_{2} = (f_{1} + b_{1} r_{1} )/D_{1}, \\ q_{n + 1} = e_{n} /D_{n} , n = 1, 2, 3, \dots , N - 2, q_{1} = e_{0} /c_{0}, \\ D_{n} = c_{n} - a_{n} q_{n - 1} + p_{n} (a_{n} p_{n - 1} - b_{n}), n = 2, 3, \dots , N, \\ \Delta _{1} = c_{1} - b_{1} p_{1} .$$Достаточными условиями устойчивости пятиточечной прогонки являются диагональное преобразование и неравенства
$$\left|{p_n}\right| + \left|{q_n}\right| \le 1, n = 1, \ldots , N - 1, \left|{p_n}\right| \le 1.$$Представим прогоночные соотношения для системы уравнений второго порядка
- B0u0 - C0u1 = D0,
Anun - 1 - Bnun + Cnun + 1 = Dn, n = 1, 2, 3, ..., N - 1,
An{u} - Bnun = Dn,
где un и Dn — векторы, Cn, An, Bn — матрицы.
Обратный ход прогонки выполняется в соответствии с формулами
un = Pn + 1un + 1 + qn + 1, un = qN+ 1.
Здесь Pn + 1 — матрица, qn + 1 — вектор.
Окончательный вид формул для прямого хода
Этот алгоритм называется методом || Pn || < 1 для $$1 \le n \le N$$ ; матрицы B_0 и (Bn - AnPn) невырождены.
Рассмотрим систему нелинейных ОДУ:
$$\left\{\begin{array}{l} \frac{d{u}}{dt} = f(u);\quad t \in [0, L], \\ F\left[{u(0), {u}(L)}\right] = 0, \quad {u}, f, F \in R^n. \\ \end{array}\right. $$$Введем пока неизвестный вектор $$\alpha$$ размерности n такой, что
и решим соответствующую задачу Коши, например, методом Рунге - Кутты; получим решение $$u(t, \alpha ).$$ Используя краевые условия, получаем нелинейную систему алгебраических уравнений
$$F(\alpha , u (L, \alpha)) = 0, \quad\mbox{или}\quad F(\alpha ) = 0,$$которая решается численно методом Ньютона или простых итераций. Отметим, что левая часть данной системы при использовании
Рассмотрим простейшую разностную аппроксимацию краевой задачи для ОДУ второго порядка
$$$ \frac{d^2u}{dt^2} = f(u), u(0) = U_1, u(L) = U_2 $$$следующего вида:
$$$ \frac{u_{n - 1} - 2u_n + u_{n + 1}}{\tau^2} = f(u_n), n = 1, \ldots , N - 1, u_0 = U_1, u_n = U_2, \tau = L/N . $$$Зададим некоторые начальные приближения $$u_n^0 = \varphi (t_n)$$ к искомой функции и построим итерационный процесс, воспользовавшись линейным приближением функции правой части
$$\begin{gather*} f(u_n^{i + 1}) \approx f(u_n^{i}) + f^{\prime}_u (u_n^{i})(u_n^{i + 1} - u_n^{i})), \\ \frac{{u_{n - 1}^{i + 1} - 2u_n^{i + 1} + u_{n + 1}^{i + 1}}}{{\tau ^2 }} = f(u_n^{i}) + f^{\prime} (u_n^{i})(u_n^{i + 1} - u_n^{i}), u_n^0 = \varphi (t_n), i = 0, 1, \ldots \end{gather*}$$где i является итерационным индексом, а начальное приближение $$\varphi (t_n)$$ удовлетворяет граничным условиям $$\varphi(0) = U_1, \varphi(L) = U_2.$$
На первой итерации, т.е. при i = 0, получаем первое
приближение к искомому решению, решая методом трехточечной
Полученное первое приближение вновь подставляем в линеаризованную
разностную задачу и
Для простоты выше рассматривались краевые задачи, для которых задавались граничные условия первого рода, т.е. на границе рассматриваемой области задавалось значение функции. Тогда трудностей с аппроксимацией граничных условий не возникало. Несколько сложнее обстоит дело при задании граничных условий второго и третьего рода (условий на производные и смешанные условия). Ограничимся описанием способов аппроксимации лишь для граничных условий второго рода.
Рассмотрим вначале линейную задачу с постоянными коэффициентами
$$\begin{gather*} \frac{d^2 u}{dt^2} + h\frac{du}{dt} + su = f, \\ u^{\prime}(0) = a, u^{\prime}(1) = b, \end{gather*}$$и соответствующую ей разностную задачу
$$$ \frac{u_{n + 1} - 2u_n + u_{n - 1}}{\tau ^2 } + h\frac{u_{n + 1} - u_{n - 1}}{2\tau } + su}_n = f_n , $$$заменив при этом производные в граничных условиях по формулам односторонего дифференцирования первого порядка
$$u_{1} - u_{0} = \tau a, u_{n} - u_{ N - 1} = \tau b.$$В результате получили разностную схему для аппроксимации дифференциальной задачи. Несмотря на то, что во внутренних узлах разностные формулы приближают дифференциальные уравнения со вторым порядком аппроксимации, решение по этой схеме получается только с первым порядком из - за понижения порядка аппроксимации в граничных точках. Естественно было бы повысить порядок схемы до второго во всех точках, включая граничные.
Наиболее распространены три способа.
0 и N + 1, находящиеся формально за пределами рассматриваемой области. Пусть в них тоже определено значение сеточной функции. Тогда узел с индексом 0 выступает как внутренний узел, и в нем может быть записано соотношение$$$ \frac{{u_1 - 2u_0 + u_{- 1}}}{{\tau ^2 }} + h\frac{{u_1 - u_{- 1}}}{{2\tau }} + su}_0 = f_0 $,$$
при этом можно для граничного условия использовать формулу для вычисления производной со вторым порядком аппроксимации по формуле центральных разностей $$$ \frac{{u_1 - u_{- 1}}}{\tau } = a $.$$ Из этих двух уравнений осталось только исключить значение в фиктивном узле (фиктивной ячейке).
Такой метод аппроксимации граничных условий очень хорошо зарекомендовал себя при решении нелинейных задач для уравнений в частных производных. В таком случае значения сеточной функции в фиктивных точках не исключаются из системы, а тоже находятся численно.
На правом конце отрезка формулы аналогичны. Не составляет труда обобщить этот подход и для граничных условий третьего рода.
$$$ \frac{{u_1 - u_0 }}{\tau } = u^{\prime} (0) + \frac{\tau }{2}u^{\prime\prime}(0) + o(\tau) $.$$
Используя само дифференциальное уравнение, можно выразить значение второй производной функции в окрестности границы:
$$$ \frac{{d^2 u(0)}}{{dt^2 }} = f(0) - h\frac{{du(0)}}{dt} - su}(0), $$$откуда следует
$$$ \frac{{u_1 - u_0 }}{\tau } = a + \frac{\tau }{2}(f(0) - h\frac{{du(0)}}{dt} - su}(0)) + o(\tau ). $$$Заменив в последнем соотношении производную на конечную разность, получим
$$$ \frac{{u_1 - u_0 }}{\tau } = a + \frac{\tau }{2}(f_0 - h\frac{{u_1 - u_0 }}{\tau } - su}_0 ) . $$$Осталось только привести это соотношение к виду, удобному для использования
в
В случае линейных уравнений с постоянными коэффициентами два последних способа приводят к одному и тому же выражению для значения сеточной функции в граничном узле. Для нелинейных уравнений эти способы будут различаться, но каждый из них приводит к повышению порядка аппроксимации граничных условий и, следовательно, разностной схемы.
Краевые задачи на собственные значения достаточно часто встречаются в физических приложениях. Например, это задача определения собственных колебаний струны, сводящаяся к ОДУ вида
$$$ \frac{d}{dt} [k(t)\frac{du}{dt}] + \lambda r(t)u = 0. $$$В приведенном уравнении краевые условия зависят от способа закрепления струны.
Это и задачи собственных колебаний упругого стержня (ОДУ четвертого порядка), нахождения энергетических уровней атома водорода, вычисление критических нагрузок в теории стержней и оболочек и др.
В задачах на собственные значения добавляется еще один пристрелочный параметр — $$\lambda,$$ поэтому эти задачи часто решаются
Если не учитывать правое краевое условие, то получим задачу Коши. Ее численное интегрирование приводит к некому значению на правом конце, зависящему от $$\lambda$$ и, вообще говоря, не равному u2. Варируя параметр $$\lambda,$$ можно добиться выполнения правого краевого условия с некоторой заданной точностью. При этом, разумеется, используются методы численного нахождения корней алгебраического уравнения, обычно
Второй пример — краевая задача на собственные значения для ОДУ второго порядка с нулевыми краевыми условиями:
$$\begin{gather*} \frac{d^2 u}{dt^2 } + {b^{\prime}}(t)\frac{du}{dt} + [{c^{\prime}}(t) + \lambda ]u = 0, \\ u(0) = u(L) = 0. \end{gather*}$$Поскольку это уравнение второго порядка с неизвестным параметром $$\lambda$$ (собственное значение дифференциального оператора), то для его решения требуется третье условие. Однако в силу линейности и однородности задачи решение определяется с точностью до произвольного постоянного множителя, что и является неявным заданием третьего условия. Его можно задать, например, следующим образом:
$$$ \frac{du(a)}{dt} = 1. $$$Трудности при использовании
В качестве тестового примера для сравнения с результатами численного расчета удобно использовать модельную краевую задачу на собственные значения:
$$$ \frac{d^2 u}{dt^2 } + \lambda u = 0, u(0) = u(L) = 0, $$$имеющую точное решение
$$$ \lambda_k = \left({\frac{{\pi k}}{L}}\right)^2, u_k (t) = \sin \left({\frac{{\pi k t}}{L}}\right), k = 1, 2, \ldots $$$Непосредственной подстановкой показывается, что решениями соответствующей разностной задачи на собственные значения
$$$ \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }} + \lambda u_n = 0, n = 1, 2, \ldots , N - 1, \tau N = L, u_0 = u_n = 0, $$$являются собственные значения и собственные функции
$$$ \lambda_k = \frac{4}{\tau^2} \sin^2 \frac {\pi k\tau }{2N}, \quad u^{\tau}_j = \sin \frac{\pi {kt}_n}{L}, \quad k, n = 1, 2, \ldots , N - 1, $$$откуда видно, что
$$$ \lim \limits_{\tau \to 0} \frac{4}{\tau^2} \sin \frac{\pi k\tau}{2L} = (\frac{\pi k}{L})^2 = \lambda_k, $$$т.е. имеет место сходимость решения разностной задачи к решению дифференциальной задачи $$\mathop {\lambda_k^\tau \to \lambda_k}\limits_{\tau \to 0} .$$
Рассмотрим разностную задачу
$$$ \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }} = - f_n , n = 1, \ldots , N - 1, u_0 = u_n = 0, \tau N = L . $$$Построим ее решение в виде разложения по базису из собственных функций разностного оператора:
$$$ {P}^\tau (u) = \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }}. $$$Этот оператор имеет полную ортонормированную систему собственных функций $$$ \omega_k (t_n) = \sqrt {\frac{2}{L}} \sin \frac{{\pi k(\tau n)}}{L}, k = 1, \ldots , N - 1, \tau n = t_n $,$$ которые соответствуют собственным значениям оператора $$P^{\tau }(u)$$: $$$ \lambda_k = \frac{4}{\tau^2 }\sin ^2 \frac{\pi k\tau }{2L} $.$$ Будем искать решение в виде
$$u_n = u(t_n) = \sum\limits_{k = 1}^{N - 1}c_k \omega_k(t_n), \quad n = 1, \ldots , N - 1,$$где ck — пока неизвестные
Подставим выражение для un, fn в виде сумм Фурье в исходное разностное уравнение
или
$$\sum\limits_{k = 1}^{N - 1}{c_k}{P}^\tau \left[{\omega_k(t_n)}\right] = - \sum\limits_{k = 1}^{N - 1}{\hat {f_k}}\omega_k (t_n),$$откуда, с учетом соотношения
$$P^{\tau }(\omega _{k}) = - \lambda \omega _{k},$$получим
$$$ - \sum\limits_{k = 1}^{N - 1}{c_k} (\lambda_k \omega_k) = - \sum\limits_{k = 1}^{N - 1}{\hat f_k}\omega_k , c_k \lambda_k = \hat f_k , c_k = \frac{{\hat f_k}}{{\lambda_k}}. $$$Таким образом, получено решение разностного уравнения в виде суммы Фурье.
Несложные арифметические подсчеты показывают, что O(N) арифметических действий для численного решения этой же разностной задачи (из них 2(N - 1) умножений и N - 1 деление) оказывается более экономичным, чем O(N2) действий ( 2(N - 1)2 умножений и N - 1 деление). Преимущества
$${y^{\prime\prime} = f(x, y), y(0) = a, y(1) = b, }$$
где нелинейная функция f не зависит явно от первой производной y'x. В 1924 году Б.Нумеров предложил следующий метод аппроксимации задачи (10.3):
где введено обозначение fk = f(xk, yk).
В чем заключается отличие метода Нумерова от аппроксимации вида (10.5):
$${y_{n + 1} - 2y_n + y_{n - 1} = h^2 f_n?}$$Описать алгоритмы численного решения нелинейных алгебраических систем (10.4) и
(10.5). В случае (10.4) и f = f(x) это — алгоритм прогонки.
Решение. Выпишем главный член погрешности аппроксимации разностного уравнения (10.4). Для этого подставим в разностное уравнение проекцию на сетку точного решения задачи (10.3). Следует отметить, что конкретный вид решения не важен, достаточно только, чтобы решение существовало. Предположим также, что оно четырежды непрерывно дифференцируемо.
Раскладывая проекции точного решения в правой части (10.4) в ряд Тейлора до четвертого порядка включительно, убедимся, что все нечетные производные взаимно уничтожатся, а четные дадут следующее выражение для главного члена погрешности аппроксимации:
$$$ y_{n + 1} - 2y_n + y_{n - 1} = y^{\prime\prime}(x_n) + \frac{1}{12}h^2 \frac{{d^4 u}}{{dx^4 }}(x_n) + O(h^4 ). $$$Из уравнения (10.3) следует, что во всех внутренних точках области выполняется равенство
$$$ y^{(4)} = \frac{{d^2 }}{{dx^2 }}f(x, y), $$$переходя в последнем равенстве к разностной аппроксимации правой части, можно учесть явно главный член погрешности аппроксимации в (10.5). Приводя подобные слагаемые в (10.4), получаем аппроксимацию Нумерова
$$$ y_{n + 1} - 2y_n + y_{n - 1} = \frac{h^2 }{12}(f_{n + 1} + 10f_n + f_{n - 1}), $$$которая приближает исходную задачу во внутренних точках сеточной области с четвертым порядком.
Такая идея разумного распоряжения правой частью для неоднородных и нелинейных задач приводит к компактным (возможно, название не слишком удачное — термины "компакт" и "компактный" уже давно заняты в
математике под совсем другое!) разностным схемам — схемам повышенного
порядка точности на нерасширенном шаблоне. Действительно, и в элементарном
В случае, когда правая часть явно зависит от первой производной, либо не
получается компактной схемы, либо схема перестает быть экономичной. Действительно, чтобы не ухудшить порядок аппроксимации, необходимо вычислять значение первой производной в соответствующих узлах со вторым порядком. Если во всех точках использовать формулу с центральной разностью, то расширится шаблон схемы — в каждом шаге элементарных вычислений должно теперь участвовать пять точек. Если для точки с индексом n использовать формулу с центральной разностью, а для точек с индексами n + 1 и n - 1 — соответствующие формулы для односторонней производной, то для каждой точки шаблона придется вычислять заново значения функции f. При использовании классического вида аппроксимации Нумерова при каждом элементарном вычислении производится лишь однократное обращение к функции вычисления правой части — лишь для fn + 1, значения fn и fn - 1 уже получены при вычислениях для точки с индексом n - 1. Так как время вычислений, как правило, определяется в таких задачах именно количеством обращений к правой части, то вычисления замедляются в три раза.
y'' + p(x)y' + q(x)y = f(x)
где f, q, p — заданные функции,
используется разностная схема.
Предложить модификацию
Решение. Рассмотрим следующую краевую
y'' + p(x)y' + q(x)y = f(x)
с условиями периодичности
y(p)(0) = y(p)(1),
где p принимает значения 0 и 1. Отметим, что пока период решения считается известным. Отрезок [0, 1] возможно рассматривать без ограничения общности, так как любой отрезок можно перевести в единичный неособым линейным преобразованием, при этом вид уравнения существенно не изменится (функции p и q умножатся на постоянный множитель).
Краевая задача с условиями периодичности решается методом циклической
Очевидно, что оно аппроксимирует дифференциальную задачу со вторым порядком на равномерной сетке. На случай неравномерной сетки рассматриваемый метод легко обобщается.
В силу периодичности дискретное уравнение должно выполняться во всех точках
сетки, включая граничные. Кроме того, в силу граничного условия при p = 0, y0 = yN + 1, система сеточных соотношений примет вид:
где $$a_k = 1 - 0, 5p_k h, b_k = 2 - q_k h^2, c_k = 1 + 0, 5p_k h, \varphi_k = f_k h^2 .$$ Матрица системы линейных уравнений получается "почти трехдиагональной" — от трехдиагональной ее отличают всего два элемента в углах матрицы.
Обобщение стандартных прогоночных соотношений (для трехточечной прогонки) на периодический случай будет иметь вид
$$y_{n - 1} = \alpha_n y_n + \beta_n + \gamma y_n.$$Из приведенного выше соотношения для y0 сразу получаем, что $$\alpha_1 = {{c_0 }/ {b_0, }}$$ $$\beta_1 = - {{f_0 }/{b_0 }},$$ $$\gamma _{1} = a_{0}/ b_{0}.$$ Теперь несложно получить рекуррентную зависимость для прогоночных коэффициентов:
По приведенным выше формулам получаются значения коэффициентов для всех
уравнений с номерами меньше, чем n. Подставим теперь прогоночные
соотношения в последнее уравнение линейной системы. В итоге с учетом введенных выше обозначений получаем
а это соотношение, сгруппировав члены, можно переписать как:
$$y_{n} = \mu _{n}y_{0} + \eta _{n},$$где введены обозначения
$$$ \mu_n = \frac{- c_n}{a_n(\alpha_n + \gamma_n) - b_n}, \eta_n = \frac{\varphi_n - a_n \beta_n}{a_n (\alpha_n + \gamma_n) - b_n}. $$$Теперь выражение для значения сеточной функции yn - 1 подставляем в прогоночные соотношения. Получается выражение, связывающее yn - 1 с y0:
Отсюда получаем следующие рекуррентные соотношения:
$$\mu_{n - 1} = \alpha_n \mu_n + \gamma_n \mu_n , \eta_{n - 1} = \beta_n + \alpha_n \eta_n + \gamma_n \eta_n.$$Отметим, что эти коэффициенты вычисляются в обратном порядке — аналог обратного хода прогонки. Последнее соотношение приводит к явному выражению для y0. В результате получаем
Теперь информации для определения значений искомой функции во всех точках сетки (еще один ход прогонки) достаточно.
Алгоритм периодической прогонки был предложен А.А.Абрамовым.
$${y^{\prime\prime} = e^{y}, \quad y (0) = a, \quad y (1) = b, }$$
$${y^{\prime\prime} = - e^{y}, \quad y (0) = a , \quad y (1) = b.}$$a = 1 и сделав замену y'x = p(y).a = 1 и различных b. Что происходит при 0 < b < 1, 499719998? При b > 1,499719998? [10.6, c. 110].$$\varepsilon y^{\prime\prime} = (y^{\prime})^2, \quad y(0) = 1, y(0) = 0, 0 < \varepsilon \ll 1.$$
y'x = p(y).где $$q \in N$$ — натуральное число, $$0 < \varepsilon \ll 1.$$
Получить численное решение задачи в случаях:
u(x) = x2, A > 1, B> 1,u(x) = x2, A = 1, B > 1 или A > 1, B = 1,u(x) = | x |, A > 1, B > 1,u(x) = | x |, A = 1, B > 1.Что происходит с решением при увеличении q? (В численных
расчетах задать $$\varepsilon = 10^{ - 2}, 10^{ - 3}, 10^{ - 4}).$$
Теоретически задача исследована в [10.8, c. 170 - 171]. В случае u(x) = | x | появляется внутренний пограничный слой — узкая область в окрестности x = 0, где y отличны от | x |.
$$\begin{gather*} \varepsilon y^{\prime\prime} = y - y^3, \\ y(0) = A, \quad y(1) = B, \left|A\right| < \sqrt 2, \quad \left|B\right| < \sqrt 2, \quad 0 < \varepsilon \ll 1. \end{gather*}$$
Решением этой задачи являются так называемые пиковые, или пичковые,
структуры. В [10.8, с. 171 - 174] исследованы свойства решений и приведены графики восьми линейно независимых решений. Там же показано, что при фиксированных A и B существует по четыре линейно независимых решения, таких, что
при всех x, за исключением точек $$x_i = {i/n}, i = \overline{1, n -
1}, n \in N (n \ge 2),$$ где $$\lim\limits_{\varepsilon \to 0}y(x, \varepsilon ) = \sqrt 2 .$$
Найти численно такие структуры для n = 2, n = 3, выбирая
соответствующее начальное приближение при линеаризации по Ньютону. (Положить $$\varepsilon = 10^{ - 2}, 10^{ - 3}, 10^{ - 4}$$ ).
имеет решение с внутренним пограничным слоем в точке x = 1/2
[10.8, c. 175]. Исследовать, как его толщина зависит от параметра $$\varepsilon.$$
Какое начальное приближение надо использовать при решении задачи методом линеаризации?
Известно, что при A = 0, В = 0 данная краевая задача имеет еще два решения, кроме тривиального ( $$y \equiv 0$$ ). Найти их численно.
Всего у этой задачи счетное множество решений.
в зависимости от вида функции a(x). Рассмотреть поведение решения
при $$\varepsilon \to 0.$$ Удалось ли получить пограничный слой типа
всплеска?
y'' = - k2y, y (0) = y (1) = 0.
n — положительное целое число).k. С какими трудностями пришлось столкнуться? Как можно улучшить используемый алгоритм?y'(0) = y (1) = 0.$$\begin{gather*} \frac{d^2 y}{{dr}^2 } + \frac{1}{r}\frac{dy}{dr} = - k^2 y, \\ y(0) = 1, y(1) = 0 \end{gather*}$$
Собственными функциями этой задачи являются цилиндрические функции Бесселя нулевого порядка, а собственными значениями задачи будут нули этих функций:
k1 = 2, 404826 , k2 = 5, 520078 , k3 = 8, 653728 , k4 = 11, 791538.
Показать, что подстановка $$$ \tilde {y} = \frac{y}{\sqrt{r}} $$$ приводит уравнение к виду, для которого можно применять метод Нумерова. Решить численно
Найти численно все уровни энергии частицы в потенциальной яме с потенциалом $$U(x) = - 2{\mathop{\mathrm{sech}}\nolimits}^2 x$$ и соответствующие им функции распределения.
Указание: уровни энергии есть собственные значения $$\lambda _{k}$$ уравнения Шредингера $$y'' + (\lambda - U(x))y = 0$$ с условиями $$y(+ \infty ) = y(- \infty ) = 0,$$ а соответствующие им собственные функции и есть функции распределения.
Данный результат важен для решения уравнения Кортвега - Де Фриза, поэтому обсуждается в [10.9, c. 11 - 17].
x = 0. При этом собственные функции системы будут четными или нечетными относительно x = 0, причем четность или нечетность чередуется с ростом квантового числа (энергии). Проверить этот эффект численно. Каким способом в этом случае можно в два раза сократить объем вычислений при расчете собственных значений для потенциалов?Рассмотреть движение частицы в поле с потенциалом Тоды ([10.10, c. 87 - 89]):
$$$ \ddot {x} = 1 - e^{x} $$$Это уравнение можно трактовать как движение частицы в поле с потенциалом U(x) = ex - x.
Найти численно все периодические решения, удовлетворяющие следующим граничным условиям:
$$\begin{gather*} x(0) = x(120) = 0, \\ \dot {x}(0) = \dot {x}(120) = A \end{gather*}$$и дополнительному условию $$A \ge 10.$$
Как период колебаний зависит от A? Сколько решений получается?
Рассмотрим линейную систему ОДУ первого порядка
$$$ \frac{du}{dt} = Au + f, u \in R^n, t \in [0, L] $$$с краевыми условиями
Ru(0) + Su(L) = q,
где u, f, q — n - мерные векторы, A(t), R(t), S(t) — матрицы размера n x n.
Для приближенного решения задачи введем расчетную сетку $$\{t_n\}_{n = 0}^{N}$$ и за приближенное решение примем сеточную функцию $$\{u_n\}_{n = 0}^{N}.$$ Рассмотрим методы построения приближенного решения.
Метод построения фундаментальных решений аналогичен известному по курсу дифференциальных уравнений способу построения общего решения системы линейных уравнений первого порядка. Решение представляется в виде
$$$ u(t) = \bar{u}(t) + \sum\limits_{k = 1}^n {\alpha_k u^k}. $$$Здесь uk(t) есть полная фундаментальная система решений
однородной задачи
с начальными данными, например,
uk(0) ={0, ..., 0, 1, 0, ..., 0}T,
где единица стоит на k месте, т.е. в качестве начальных данных
используются векторы uk(0) = Ek. Важно, чтобы решения однородной задачи составляли систему линейно независимых функций. Каждая такая функция ищется численно как решение соответствующей задачи Коши, используя методы, описанные в лекциях 8 и 9.
Пусть $$$ \bar{u}(t) $$$ — частное решение неоднородной системы
$$$ \frac{d\bar{u}}{dt} = A\bar{u} + f $$$с нулевыми начальными условиями $$$ \bar{u}(0) = 0 $.$$ Тогда неопределенные коэффициенты $$\alpha_k$$ находятся из краевых условий
$$\begin{gather*} Ru(0) + Su(L) = q, \mbox{ или }\\ R[\sum\limits_{k = 1}^n{\alpha_k{u}^k (0)}] + S\bar {u}(L) + S[\sum\limits_{k = 1}^n{\alpha_k{u}^k(L)}] = q. \end{gather*}$$Последнее соотношение представляет собой СЛАУ относительно коэффициентов $$\alpha_{k}$$ размерности n.
Полную фундаментальную систему решений однородной задачи можно получить,
используя, например, схему второго порядка точности (
Частное решение получаем аналогично:
$$\begin{gather*} \frac{\bar{u}_{n + 1} - \bar{u}_n}{\tau } = A_{n + 1/2} \frac{\bar{u}_{n + 1} + \bar{u}^{\prime\prime}_n}{2} + f_{n + 1/2}, \\ f_{n + 1/2} = F\left({t_n + \frac{\tau }{2}}\right), \\ \bar{u}_0 (0) = 0, \mbox{или}\\ \bar{u}_{n + 1} = \left({E - \frac{\tau}{2}A}\right)^{- 1}\left[{\left({E + \frac{\tau}{2}A}\right)u_n + \tau f_{n + 1/2}}\right], \bar{u}_0 = 0. \end{gather*}$$Приведем пример, когда решение задачи методом построения фундаментальных решений не проходит. Рассмотрим систему уравнений
$$\begin{gather*} \frac{du}{dt} = {av} + f, \frac{dv}{dt} = bu + g, \\ u(0) = U_0, v(1) = 0 , \end{gather*}$$Подобные системы ОДУ моделируют, например, процессы прохождения излучения или потоков элементарных частиц (гамма - излучение, потоки нейтронов) через разные среды (атмосфера, защита ядерных реакторов) в приближении оптически толстого слоя. Коэффициенты $$a, b \sim 50$$ характерны для защиты реакторов. Найдем общее решение этой задачи с помощью метода фундаментальных систем в виде линейной комбинации двух решений однородных ОДУ:
$$\left( \begin{array}{l} u \\ v \\ \end{array} \right) = \left( \begin{array}{l} {\bar {u}} \\ {\bar {v}} \\ \end{array} \right) + \alpha_1 \left( \begin{array}{l} u^1 \\ v^1 \\ \end{array} \right) + \alpha_2 \left( \begin{array}{l} u^2 \\ v^2 \\ \end{array} \right),$$где (u1, v1) и (u2, v2) — решения двух однородных систем (в дальнейшем для простоты будем полагать, что $$$ \bar {u} = \bar {v} = 0 $.$$) Тогда
а коэффициенты $$\alpha_1$$ и $$\alpha_2$$ находятся из краевых условий.
Общее решение такой системы, как известно, представляет собой сумму двух экспонент
$$\left( \begin{array}{l} u^1 \\ v^1 \\ \end{array} \right) = C_1 \left( \begin{array}{l} a_1 \\ b_1 \\ \end{array} \right)e^{\lambda_1 t} + C_2 \left( \begin{array}{l} a_2 \\ b_2 \\ \end{array} \right)e^{\lambda_2 t}$$(аналогичным образом решение представляется и для u2, v2 ), где ai, bi, i = 1, 2 находятся из решения задачи, Ci; i = 1, 2 — произвольные постоянные; $$\lambda _{1}$$ и $$\lambda _{2}$$ — корни характеристического уравнения$$\left| \begin{array}{cc}
{- \lambda} {a} \\
{b} {- \lambda} \\
\end{array} \right| = 0.$$
Решая U0e - 50t (в случае задачи с защитой ректора U0 — падающий поток нейтронов; задача защиты — значительное его ослабление, приблизительно в e50 раз).
Слагаемые с быстрорастущими экспонентами должны взаимно уничтожиться.
Получение же численного решения - весьма трудная задача, поскольку численное решение имеет большую и быстро возрастающую погрешность. Пусть $$u^{M} = u(1 + \varepsilon ) \sim e^{50}t(1 + 10^{ - 12}),$$ т.е. начальная погрешность имеет порядок 10 - 12. Она возрастает при вычислениях примерно в e50t раз — даже при умеренных t это очень большое число.
Рассмотрим метод
где $$\alpha$$ и $$\beta$$ — пока неизвестные функции (прогоночные коэффициенты), для которых необходимо получить дифференциальные уравнения.
Продифференцируем это соотношение
$$$ \dot {u} = \dot {\alpha} v + \alpha \dot {v} + \dot {\beta} $$$и подставим в него уравнения системы $$$ \dot {u} = {av} + f, \dot {v} = {bu} + g $.$$ В результате получим, что $$$ {av} + f = \dot {\alpha} v + \alpha ({bu} + g) + \dot {\beta} $.$$ Подставим в полученное соотношение уравнение $$u = \alpha v + \beta .$$ Тогда $${av} + f = \dot {\alpha} v + b\alpha ^2 v + b\beta \alpha + \alpha g + \dot {\beta} v + \beta .$$
После приведения подобных членов имеем равенство
$$$ v(\dot {\alpha} + b\alpha ^2 - a) + (\dot {\beta} + b\beta \alpha + \alpha g - f) = 0. $$$Приравнивая к нулю коэффициенты при v и единице, получим два
дифференциальных уравнения для прогоночных коэффициентов:
Дополним их начальными условиями. Левое краевое условие вида u (0) =
U0 запишем в виде прогоночного соотношения $$u(0) = \alpha (0) v(0) + \beta(0),$$ полагая $$\alpha (0) = 0, \beta(0) = U_0.$$ Таким образом, получаем начальные данные для двух задач Коши для $$\alpha (t)$$ и $$\beta (t),$$ которые могут быть решены численно.
Теперь разрешим правое краевое условие. На правой границе отрезка интегрирования имеем условие v(1) = 0 и прогоночное соотношение при t = 1: $$u (1) = \alpha (1)v(1) + \beta (1),$$ откуда получаем $$u(1) =
\beta (1).$$
Далее воспользуемся уравнением $$$ \dot {v} = {bu} + g $,$$ подставив в него прогоночное соотношение $$u = \alpha v + \beta,$$ получим дифференциальное уравнение для v:
Интегрируя эту задачу справа налево, попутно определяем u(t):
Метод
Рассмотрим теперь общую постановку
где $$u, q_i , D_i, f \in R^n, \left\|{D_i}\right\| \simeq O(1), i = 1, \ldots,.$$
Здесь $$\mathbf{A}$$ — постоянная матрица размером N
x N (дальнейшие рассуждения будут справедливы и для систем уравнений с переменными коэффициентами). Определители систем линейных алгебраических уравнений, которыми являются краевые условия на обоих концах интервала интегрирования, полагаются отличными от нуля.
Определение. Рассматриваемая
Отношение $${{\Lambda_0 }/{\lambda_0 }} \gg 1$$ является параметром, характеризующим жесткость системы. В дальнейшем будем полагать $$\Lambda _0 (b - a) \gg 1,$$ $$\lambda_0 (b - a) \simeq O(1).$$
Общее решение такой системы имеет вид
$${u}(t) = \sum\limits_{i = 1}^{N_1 }{\gamma_i^1 e^{\Lambda_i^1t}{{\Omega }}_i^1 + \sum\limits_{i = N_1 + 1}^{N_2 }{\gamma_i^2 e^{\Lambda_i^2 t}{{\Omega }}_i^2 + \sum\limits_{i = N_2 + 1}^{N}{\gamma_i^3 {e}^{\lambda_it}} }}{{ \omega }}_i,$$где $$\Omega_i^1, \Omega_i^2, \omega_i$$ есть собственные векторы матрицы A, соответствующие трем частям спектра. Понятна качественная структура этого решения, содержащего как левый, так и правый пограничные слои. Будем полагать, что количество собственных значений в каждой из трех частей спектра не изменяется. Особенность || f || и || q || - нормы правых частей в системе ОДУ и краевых условиях, соответственно. Для численного решения задачи
можно использовать такую же схему второго порядка, как и ранее. Рассматриваем класс вычислительно корректных задач, для которых $$C = O(1) \ll \exp(\Lambda_0(b - a)).$$
В дальнейшем будем полагать величину $$||A|| (b - a) \approx \Lambda _{0}(b - a)$$ большой, а величину $$exp(||A|| (b - a)) \approx exp(\Lambda _{0}(b - a))$$ — очень большой. В приложениях такие задачи встречаются наиболее часто. Важно отметить и то, что не все возможные постановки задач для жесткой системы приводят к вычислительно корректным алгоритмам. Показывается, что необходимыми (и почти достаточными) условиями корректности являются следующие неравенства: $$k \ge N_1, (N - k) \ge N_2 - N_1,$$ т.е. число краевых условий на левом конце отрезка интегрирования не должно быть меньше быстро убывающих вправо решений, на правом конце — не меньше быстро убывающих влево решений. В противном случае
Проблемы, которые возникают при численном решении eLt, как известно, приводит к потере точности и накоплению вычислительных ошибок. Вторая проблема состоит в следующем. Для вычисления коэффициентов $$\alpha_i,$$ входящих в общее решение неоднородной системы (а оно состоит из суммы частного решения неоднородной системы и общего решения однородной $$u (t) = {\bar{u}}_0 + \sum\limits_{i = 1}^{N}{\alpha_i}u_i$$ ), приходится решать плохо обусловленную СЛАУ $$(D_i, \bar{u}(b) + \sum\limits_i {\alpha_i}u_i(b)) = q_i, i = k + 1, \ldots,.$$
где коэффициенты g, h, s, вообще
говоря, являются функциями независимого переменного t.
Для разностной аппроксимации рассматриваемой краевой задачи введем равномерную разностную сетку $$\{t_n\}_0^{N}, t_n = n\tau , \tau = (b - a)/N$$ и определим на этой сетке сеточную функцию $$\{u_n\}_0^{N}.$$
Коэффициент g, вообще говоря, может не иметь первой производной. Такая задача может возникнуть, например, в случае расчета установившегося распределения температуры в задаче стационарной теплопроводности с контактным разрывом. Представим разностную задачу в виде
Контактный разрыв при этом помещается в узел с номером n.
В этом случае фактически аппроксимируется тепловой поток через границы ячейки разностной сетки, для уравнения теплопроводности получается консервативная разностная схема — подробнее смотри в лекциях, посвященных
Выписанные выше соотношения определяют простейшую разностную схему. Под разностной схемой здесь и ниже понимается совокупность разностных уравнений для определения значений сеточной функции внутри расчетной области, дополненная соответствующими начальными и граничными условиями для этой сеточной функции.
Для определения значений сеточной функции получается СЛАУ с трехдиагональной матрицей
- b0u0 + c0u1 = d0, anun - 1 - bnun + cnun - 1 = dn , n = 1, ..., N - 1, anuN - 1 - bnun = dn,
где $$\begin{gather*} a_n = \frac{{g_{n - 1/2}}}{{\tau ^2 }} - \frac{{h_n}}{{2\tau }}, c_n = \frac{{g_{n + 1/2}}}{{\tau ^2 }} + \frac{{h_n}}{{2\tau }}, \\ b_n = a_n + c_n - h_n , d_n = f_n , b_0 = \frac{A}{\tau } - \frac{B}{2}, c_0 = \frac{A}{\tau } + \frac{B}{2}, d_0 = D \end{gather*}$$
Эта СЛАУ представима в каноническом виде A'u = D,
где A' — матрица
u, d есть векторы - столбцы
D = (d0, d1, ..., dn)T.
Трехдиагональные матрицы часто возникают при численном решении краевых задач как для обыкновенных дифференциальных уравнений, так и для уравнений в частных производных. Ранее матрица подобной структуры встречалась при построении сплайна Шонберга (лекция 6). Характерная особенность таких матриц заключается в том, что при большой размерности матрица имеет ленточную структуру — все элементы вне ленты (n - го порядка требует O(N3) арифметических действий и O(N2) ячеек памяти. В численных методах большую роль играют экономичные алгоритмы, в которых количество арифметических операций пропорционально количеству неизвестных — O(N).
К таким алгоритмам относится метод трехточечной разностной
Решение ищется в виде прогоночного соотношения
un - 1 = pnun + qn, n = 1, ..., N,
где pn и qn — прогоночные коэффициенты, подлежащие определению.
Левое краевое условие также записывается в виде прогоночного соотношения
$$$ u_0 = \frac{c_0 }{d_0 }u_1 - \frac{d_0 }{b_0}, $$$где p1 = c0/ b0, q1 = - d0/ b0 (заметим, что если A > 0 и B < 0, то b0 > c0, 0 < p < 1 ).
Получим рекуррентные формулы, позволяющие последовательно вычислить p2, q2, p3, q3 и т.д. вплоть до pn, qn.
Подставив равенство u n - 1 = pnun + qn в уравнение anun - 1 - bnun + cnu n + 1 = dn, получим an(pnun + qn) - bnun + cnun + 1 = dn, или
Сравнивая эту запись со стандартным видом прогоночного соотношения
un = pn + 1un + 1 + qn + 1,
видим, что для прогоночных коэффициентов должны выполняться равенства
$$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n}, \quad q_{n + 1} = \frac{a_n q_n - d_n}{b_n - a_n p_n}. $$$Эти формулы определяют прямой ход прогонки.
Из краевого условия на правом конце отрезка интегрирования anuN - 1 -
bNuN = dN и прогоночного соотношения uN - 1 = pNuN + qN находим величину uN.
Далее последовательно вычисляются остальные неизвестные un n = N - 1, ..., 1 un - 1 = pnun + qn. Это — обратный ход алгоритма прогонки.
Теперь, после того как описан прогоночный алгоритм, исследуем его на устойчивость.
Для этого рассмотрим вычисление прогоночного коэффициента pn (т.е. этап прямого хода прогонки). В идеальной арифметике этот коэффициент равен $$$ p_{n + 1} = \frac{c_n}{b_n - a_np_n} $,$$ в конечноразрядной арифметике — $$p_{n + 1}^{M} = p_{n + 1} + \Delta_{n + 1} ,$$ где $$\Delta _{n + 1}$$ — погрешность, связанная с округлениями на всех предшествующих этапах вычислений. Полагая $$\Delta _{n + 1}$$ малой, исследуем изменение этой погрешности с ростом n. Для этого запишем соотношение между $$\Delta _{n}$$ и $$\Delta _{n + 1}$$ в виде
где $$\varepsilon _{n}$$ — погрешность вычислений правой части и машинного представления коэффициентов an, bn, cn. Следовательно, $$\Delta _{n + 1}$$ складывается из двух составляющих — локальной погрешности $$\varepsilon _{n}$$ и наследственной $$\Delta _{n + 1}.$$
Полагая $$\Delta_n \ll p_n$$ при $$\Delta _{n} > 0, p_{n} > 0$$ и опуская члены порядка $$O(\Delta ^{2}),$$ получим оценку
$$$ p_{n + 1} + \Delta_{n + 1} \approx \frac{c_n}{b_n - a_n p_n} + \frac{c_n a_n}{{(b_n - a_n p_n)}^2 }\Delta_n + \varepsilon_n . $$$Отсюда, учитывая, что $$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n} $,$$ получим требуемую оценку для эволюции погрешности
$$$ {\Delta_{n + 1} = \frac{a_n}{c_n}p_{n + 1}^2 \Delta_n + \varepsilon_n, n = 0, \ldots , N - 1.} $$$Докажем следующую теорему.
Теорема. Если выполнены условия диагонального преобладания $$\left|{b_n}\right| \ge \left|{a_n}\right| + \left|{c_n}\right|$$ и хотя бы для одной строки матрицы системы имеет место строгое диагональное преобладание $$\left({\left|{b_n}\right| > \left|{a_n}\right| + \left|{c_n}\right|}\right).$$ Пусть, кроме того, 0 < p1 < 1. Тогда алгоритм прогонки устойчив.
Доказательство.
Докажем утверждение теоремы для случая, когда во всех строках матрицы выполнено условие строгого диагонального преобладания. На случай нестрогих неравенств обобщение доказательства очевидно.
0 < p1 < 1. Для определенности положим an > 0, bn > 0, cn > 0.Тогда $$$ p_{n + 1} = \frac{c_n}{b_n - a_np_n} > \frac{c_n}{b_n - a_n} > 0 $.$$
Кроме того, $$$ p_{n + 1} = \frac{c_n}{b_n - a_n p_n} < \frac{c_n}{a_n + c_n - a_n p_n} = \frac{c_n}{c_n + a_n (1 - p_n)} \le 1 $,$$ откуда следует 0 < pn + 1 < 1.
Для этого вспомним выражения для коэффициентов линейной системы, полученной при разностной аппроксимации исходного уравнения второго порядка.
$$\begin{gather*} a_n = \frac{g_{n - 1/2}}{\tau ^2 } - \frac{h_n}{2\tau } \approx \frac{g_{n - 1/2}}{\tau ^2 }, g_{n - 1/2} = g(t_n - \tau /2), \\ c_n = \frac{g_{n + 1/2}}{\tau^2} + \frac{h_n}{2\tau} \approx \frac{g_{n + 1/2}}{\tau^2}, \\ \frac{a_n}{c_n} \approx \frac{g(t - \tau/2)}{s(t + \tau/2)} = 1 + c\tau. \end{gather*}$$Вернемся к выражению для эволюции погрешности (10.2). С учетом полученных оценок имеем $$$ \Delta_{n + 1} \le \frac{a_n}{c_n}\Delta_n + \varepsilon_n $.$$ Так как $$$ \frac{a_n}{c_n} \approx 1 + c\tau $,$$ то $$\Delta_{n + 1} \le (1 + c\tau )\Delta_n + \varepsilon_n,$$ где c — константа Липшица для функции g(t). Считаем, что на каждом шаге ошибка округления не превосходит предельного значения, т.е. $$0 \le \varepsilon_n \le \varepsilon.$$
Из цепочки неравенств
$$\begin{gather*} \Delta_2 \le (1 + c\tau )^2 \Delta_0 + \varepsilon \left[{1 + (1 + c\tau )}\right], \\ \Delta_3 \le (1 + c\tau )^3 \Delta_0 + \varepsilon \left[{1 + (1 + c\tau ) + (1 + c\tau )^2 }\right], \end{gather*}$$будет следовать оценка $$$ \Delta_n \le (1 + c\tau )^{n} \Delta_0 + \frac{{(1 + c\tau )^{n}}}{{c\tau }}\varepsilon $.$$ Используя известные из математического анализа неравенства, последнюю оценку можно записать в виде
$$$ \Delta_n \le e^{c\tau n} \Delta_0 + e^{c\tau n} \frac{\varepsilon }{{c\tau }}, c\tau n = ct . $$$При $$ct \sim O(1)$$ погрешности не накапливаются (для краевых задач это реальное условие), поскольку в расчетах $$\tau \gg \varepsilon,$$ где $$\varepsilon$$ — машинный эпсилон (подробнее в лекции 1).
Аналогичное утверждение доказывается и для второго прогоночного коэффициента $$$ q_{n + 1} = \frac{a_n q_n - d_n}{b_n - a_n p_n} $;$$ алгоритм устойчив при выполнении тех же условий.
Покажем устойчивость обратного хода прогонки.
При обратном ходе вычисления проводятся по формулам un = pn + 1un + 1 + qn + 1, откуда, учитывая, что $$u_n^{M} = u_n + \Delta_n$$ и $$u_{n + 1}^{M} = u_{n + 1} + \Delta_{n + 1},$$ получим $$\Delta _{n} = p_{n + 1}\Delta _{n + 1} + \varepsilon _{n},$$ где $$\Delta _{n}$$ — наследственная погрешность, $$\varepsilon _{n}$$ — погрешность округления на n шаге. Очевидно, что обратный ход прогонки устойчив при выполнении условия 0 < pn < 1, или 0 < p1 < 1.
При численном решении краевых задач для обыкновенных дифференциальных уравнений 4 - го порядка возникают СЛАУ с пятидиагональной матрицей вида:
c0u0 - d0u1 + e0u2 = f0, - b1 u0 + c1 u1 - d1 u2 + e1 u3 = f1, anun - 2 - bnun - 1 + cnun - dnun + 1 + enun + 2 = fn, n = 2, ... , N - 2, aN - 1 uN - 3 - bN - 1 uN - 2 + cN - 1 uN - 1 - dN - 1 un = fN - 1, an uN - 2 - bn uN - 1 + cn un = fn .
Алгоритм решения таких систем — пятиточечная прогонка, формулы которой
выводятся аналогично формулам для трехточечной прогонки. Приведем их окончательный вид. В прогоночном соотношении появятся три коэффициента (p, q, r):
un = pn + 1 un + 1 - qn + 1un + 2 + rn + 1, uN - 1 = pn un + rn, un = rN + 1.
Прогоночные коэффициенты находятся по формулам
$$p_{n + 1} = [d_{n} + q_{n}(a_{n} p_{n - 1} - b_{n})]/D_{n} , n = 2, 3, \dots , N - 1, \setminus \setminus \\ p_{1} = d_{0} /c_{0}, p_{2} = (d_{1} - q_{1} b_{1} )/D_{1}, \\ r_{n + 1} = [f_{n} - a_{n} r_{n - 1} - r_{n} (a_{n} p_{n - 1} - b_{n})]/D_{n}, n = 2, 3, \dots , N, \\ r_{1} = f_{0} /c_{0}, r_{2} = (f_{1} + b_{1} r_{1} )/D_{1}, \\ q_{n + 1} = e_{n} /D_{n} , n = 1, 2, 3, \dots , N - 2, q_{1} = e_{0} /c_{0}, \\ D_{n} = c_{n} - a_{n} q_{n - 1} + p_{n} (a_{n} p_{n - 1} - b_{n}), n = 2, 3, \dots , N, \\ \Delta _{1} = c_{1} - b_{1} p_{1} .$$Достаточными условиями устойчивости пятиточечной прогонки являются диагональное преобразование и неравенства
$$\left|{p_n}\right| + \left|{q_n}\right| \le 1, n = 1, \ldots , N - 1, \left|{p_n}\right| \le 1.$$Представим прогоночные соотношения для системы уравнений второго порядка
- B0u0 - C0u1 = D0,
Anun - 1 - Bnun + Cnun + 1 = Dn, n = 1, 2, 3, ..., N - 1,
An{u} - Bnun = Dn,
где un и Dn — векторы, Cn, An, Bn — матрицы.
Обратный ход прогонки выполняется в соответствии с формулами
un = Pn + 1un + 1 + qn + 1, un = qN+ 1.
Здесь Pn + 1 — матрица, qn + 1 — вектор.
Окончательный вид формул для прямого хода
Этот алгоритм называется методом || Pn || < 1 для $$1 \le n \le N$$ ; матрицы B_0 и (Bn - AnPn) невырождены.
Рассмотрим систему нелинейных ОДУ:
$$\left\{\begin{array}{l} \frac{d{u}}{dt} = f(u);\quad t \in [0, L], \\ F\left[{u(0), {u}(L)}\right] = 0, \quad {u}, f, F \in R^n. \\ \end{array}\right. $$$Введем пока неизвестный вектор $$\alpha$$ размерности n такой, что
и решим соответствующую задачу Коши, например, методом Рунге - Кутты; получим решение $$u(t, \alpha ).$$ Используя краевые условия, получаем нелинейную систему алгебраических уравнений
$$F(\alpha , u (L, \alpha)) = 0, \quad\mbox{или}\quad F(\alpha ) = 0,$$которая решается численно методом Ньютона или простых итераций. Отметим, что левая часть данной системы при использовании
Рассмотрим простейшую разностную аппроксимацию краевой задачи для ОДУ второго порядка
$$$ \frac{d^2u}{dt^2} = f(u), u(0) = U_1, u(L) = U_2 $$$следующего вида:
$$$ \frac{u_{n - 1} - 2u_n + u_{n + 1}}{\tau^2} = f(u_n), n = 1, \ldots , N - 1, u_0 = U_1, u_n = U_2, \tau = L/N . $$$Зададим некоторые начальные приближения $$u_n^0 = \varphi (t_n)$$ к искомой функции и построим итерационный процесс, воспользовавшись линейным приближением функции правой части
$$\begin{gather*} f(u_n^{i + 1}) \approx f(u_n^{i}) + f^{\prime}_u (u_n^{i})(u_n^{i + 1} - u_n^{i})), \\ \frac{{u_{n - 1}^{i + 1} - 2u_n^{i + 1} + u_{n + 1}^{i + 1}}}{{\tau ^2 }} = f(u_n^{i}) + f^{\prime} (u_n^{i})(u_n^{i + 1} - u_n^{i}), u_n^0 = \varphi (t_n), i = 0, 1, \ldots \end{gather*}$$где i является итерационным индексом, а начальное приближение $$\varphi (t_n)$$ удовлетворяет граничным условиям $$\varphi(0) = U_1, \varphi(L) = U_2.$$
На первой итерации, т.е. при i = 0, получаем первое
приближение к искомому решению, решая методом трехточечной
Полученное первое приближение вновь подставляем в линеаризованную
разностную задачу и
Для простоты выше рассматривались краевые задачи, для которых задавались граничные условия первого рода, т.е. на границе рассматриваемой области задавалось значение функции. Тогда трудностей с аппроксимацией граничных условий не возникало. Несколько сложнее обстоит дело при задании граничных условий второго и третьего рода (условий на производные и смешанные условия). Ограничимся описанием способов аппроксимации лишь для граничных условий второго рода.
Рассмотрим вначале линейную задачу с постоянными коэффициентами
$$\begin{gather*} \frac{d^2 u}{dt^2} + h\frac{du}{dt} + su = f, \\ u^{\prime}(0) = a, u^{\prime}(1) = b, \end{gather*}$$и соответствующую ей разностную задачу
$$$ \frac{u_{n + 1} - 2u_n + u_{n - 1}}{\tau ^2 } + h\frac{u_{n + 1} - u_{n - 1}}{2\tau } + su}_n = f_n , $$$заменив при этом производные в граничных условиях по формулам односторонего дифференцирования первого порядка
$$u_{1} - u_{0} = \tau a, u_{n} - u_{ N - 1} = \tau b.$$В результате получили разностную схему для аппроксимации дифференциальной задачи. Несмотря на то, что во внутренних узлах разностные формулы приближают дифференциальные уравнения со вторым порядком аппроксимации, решение по этой схеме получается только с первым порядком из - за понижения порядка аппроксимации в граничных точках. Естественно было бы повысить порядок схемы до второго во всех точках, включая граничные.
Наиболее распространены три способа.
0 и N + 1, находящиеся формально за пределами рассматриваемой области. Пусть в них тоже определено значение сеточной функции. Тогда узел с индексом 0 выступает как внутренний узел, и в нем может быть записано соотношение$$$ \frac{{u_1 - 2u_0 + u_{- 1}}}{{\tau ^2 }} + h\frac{{u_1 - u_{- 1}}}{{2\tau }} + su}_0 = f_0 $,$$
при этом можно для граничного условия использовать формулу для вычисления производной со вторым порядком аппроксимации по формуле центральных разностей $$$ \frac{{u_1 - u_{- 1}}}{\tau } = a $.$$ Из этих двух уравнений осталось только исключить значение в фиктивном узле (фиктивной ячейке).
Такой метод аппроксимации граничных условий очень хорошо зарекомендовал себя при решении нелинейных задач для уравнений в частных производных. В таком случае значения сеточной функции в фиктивных точках не исключаются из системы, а тоже находятся численно.
На правом конце отрезка формулы аналогичны. Не составляет труда обобщить этот подход и для граничных условий третьего рода.
$$$ \frac{{u_1 - u_0 }}{\tau } = u^{\prime} (0) + \frac{\tau }{2}u^{\prime\prime}(0) + o(\tau) $.$$
Используя само дифференциальное уравнение, можно выразить значение второй производной функции в окрестности границы:
$$$ \frac{{d^2 u(0)}}{{dt^2 }} = f(0) - h\frac{{du(0)}}{dt} - su}(0), $$$откуда следует
$$$ \frac{{u_1 - u_0 }}{\tau } = a + \frac{\tau }{2}(f(0) - h\frac{{du(0)}}{dt} - su}(0)) + o(\tau ). $$$Заменив в последнем соотношении производную на конечную разность, получим
$$$ \frac{{u_1 - u_0 }}{\tau } = a + \frac{\tau }{2}(f_0 - h\frac{{u_1 - u_0 }}{\tau } - su}_0 ) . $$$Осталось только привести это соотношение к виду, удобному для использования
в
В случае линейных уравнений с постоянными коэффициентами два последних способа приводят к одному и тому же выражению для значения сеточной функции в граничном узле. Для нелинейных уравнений эти способы будут различаться, но каждый из них приводит к повышению порядка аппроксимации граничных условий и, следовательно, разностной схемы.
Краевые задачи на собственные значения достаточно часто встречаются в физических приложениях. Например, это задача определения собственных колебаний струны, сводящаяся к ОДУ вида
$$$ \frac{d}{dt} [k(t)\frac{du}{dt}] + \lambda r(t)u = 0. $$$В приведенном уравнении краевые условия зависят от способа закрепления струны.
Это и задачи собственных колебаний упругого стержня (ОДУ четвертого порядка), нахождения энергетических уровней атома водорода, вычисление критических нагрузок в теории стержней и оболочек и др.
В задачах на собственные значения добавляется еще один пристрелочный параметр — $$\lambda,$$ поэтому эти задачи часто решаются
Если не учитывать правое краевое условие, то получим задачу Коши. Ее численное интегрирование приводит к некому значению на правом конце, зависящему от $$\lambda$$ и, вообще говоря, не равному u2. Варируя параметр $$\lambda,$$ можно добиться выполнения правого краевого условия с некоторой заданной точностью. При этом, разумеется, используются методы численного нахождения корней алгебраического уравнения, обычно
Второй пример — краевая задача на собственные значения для ОДУ второго порядка с нулевыми краевыми условиями:
$$\begin{gather*} \frac{d^2 u}{dt^2 } + {b^{\prime}}(t)\frac{du}{dt} + [{c^{\prime}}(t) + \lambda ]u = 0, \\ u(0) = u(L) = 0. \end{gather*}$$Поскольку это уравнение второго порядка с неизвестным параметром $$\lambda$$ (собственное значение дифференциального оператора), то для его решения требуется третье условие. Однако в силу линейности и однородности задачи решение определяется с точностью до произвольного постоянного множителя, что и является неявным заданием третьего условия. Его можно задать, например, следующим образом:
$$$ \frac{du(a)}{dt} = 1. $$$Трудности при использовании
В качестве тестового примера для сравнения с результатами численного расчета удобно использовать модельную краевую задачу на собственные значения:
$$$ \frac{d^2 u}{dt^2 } + \lambda u = 0, u(0) = u(L) = 0, $$$имеющую точное решение
$$$ \lambda_k = \left({\frac{{\pi k}}{L}}\right)^2, u_k (t) = \sin \left({\frac{{\pi k t}}{L}}\right), k = 1, 2, \ldots $$$Непосредственной подстановкой показывается, что решениями соответствующей разностной задачи на собственные значения
$$$ \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }} + \lambda u_n = 0, n = 1, 2, \ldots , N - 1, \tau N = L, u_0 = u_n = 0, $$$являются собственные значения и собственные функции
$$$ \lambda_k = \frac{4}{\tau^2} \sin^2 \frac {\pi k\tau }{2N}, \quad u^{\tau}_j = \sin \frac{\pi {kt}_n}{L}, \quad k, n = 1, 2, \ldots , N - 1, $$$откуда видно, что
$$$ \lim \limits_{\tau \to 0} \frac{4}{\tau^2} \sin \frac{\pi k\tau}{2L} = (\frac{\pi k}{L})^2 = \lambda_k, $$$т.е. имеет место сходимость решения разностной задачи к решению дифференциальной задачи $$\mathop {\lambda_k^\tau \to \lambda_k}\limits_{\tau \to 0} .$$
Рассмотрим разностную задачу
$$$ \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }} = - f_n , n = 1, \ldots , N - 1, u_0 = u_n = 0, \tau N = L . $$$Построим ее решение в виде разложения по базису из собственных функций разностного оператора:
$$$ {P}^\tau (u) = \frac{{u_{n - 1} - 2u_n + u_{n + 1}}}{{\tau ^2 }}. $$$Этот оператор имеет полную ортонормированную систему собственных функций $$$ \omega_k (t_n) = \sqrt {\frac{2}{L}} \sin \frac{{\pi k(\tau n)}}{L}, k = 1, \ldots , N - 1, \tau n = t_n $,$$ которые соответствуют собственным значениям оператора $$P^{\tau }(u)$$: $$$ \lambda_k = \frac{4}{\tau^2 }\sin ^2 \frac{\pi k\tau }{2L} $.$$ Будем искать решение в виде
$$u_n = u(t_n) = \sum\limits_{k = 1}^{N - 1}c_k \omega_k(t_n), \quad n = 1, \ldots , N - 1,$$где ck — пока неизвестные
Подставим выражение для un, fn в виде сумм Фурье в исходное разностное уравнение
или
$$\sum\limits_{k = 1}^{N - 1}{c_k}{P}^\tau \left[{\omega_k(t_n)}\right] = - \sum\limits_{k = 1}^{N - 1}{\hat {f_k}}\omega_k (t_n),$$откуда, с учетом соотношения
$$P^{\tau }(\omega _{k}) = - \lambda \omega _{k},$$получим
$$$ - \sum\limits_{k = 1}^{N - 1}{c_k} (\lambda_k \omega_k) = - \sum\limits_{k = 1}^{N - 1}{\hat f_k}\omega_k , c_k \lambda_k = \hat f_k , c_k = \frac{{\hat f_k}}{{\lambda_k}}. $$$Таким образом, получено решение разностного уравнения в виде суммы Фурье.
Несложные арифметические подсчеты показывают, что O(N) арифметических действий для численного решения этой же разностной задачи (из них 2(N - 1) умножений и N - 1 деление) оказывается более экономичным, чем O(N2) действий ( 2(N - 1)2 умножений и N - 1 деление). Преимущества
$${y^{\prime\prime} = f(x, y), y(0) = a, y(1) = b, }$$
где нелинейная функция f не зависит явно от первой производной y'x. В 1924 году Б.Нумеров предложил следующий метод аппроксимации задачи (10.3):
где введено обозначение fk = f(xk, yk).
В чем заключается отличие метода Нумерова от аппроксимации вида (10.5):
$${y_{n + 1} - 2y_n + y_{n - 1} = h^2 f_n?}$$Описать алгоритмы численного решения нелинейных алгебраических систем (10.4) и
(10.5). В случае (10.4) и f = f(x) это — алгоритм прогонки.
Решение. Выпишем главный член погрешности аппроксимации разностного уравнения (10.4). Для этого подставим в разностное уравнение проекцию на сетку точного решения задачи (10.3). Следует отметить, что конкретный вид решения не важен, достаточно только, чтобы решение существовало. Предположим также, что оно четырежды непрерывно дифференцируемо.
Раскладывая проекции точного решения в правой части (10.4) в ряд Тейлора до четвертого порядка включительно, убедимся, что все нечетные производные взаимно уничтожатся, а четные дадут следующее выражение для главного члена погрешности аппроксимации:
$$$ y_{n + 1} - 2y_n + y_{n - 1} = y^{\prime\prime}(x_n) + \frac{1}{12}h^2 \frac{{d^4 u}}{{dx^4 }}(x_n) + O(h^4 ). $$$Из уравнения (10.3) следует, что во всех внутренних точках области выполняется равенство
$$$ y^{(4)} = \frac{{d^2 }}{{dx^2 }}f(x, y), $$$переходя в последнем равенстве к разностной аппроксимации правой части, можно учесть явно главный член погрешности аппроксимации в (10.5). Приводя подобные слагаемые в (10.4), получаем аппроксимацию Нумерова
$$$ y_{n + 1} - 2y_n + y_{n - 1} = \frac{h^2 }{12}(f_{n + 1} + 10f_n + f_{n - 1}), $$$которая приближает исходную задачу во внутренних точках сеточной области с четвертым порядком.
Такая идея разумного распоряжения правой частью для неоднородных и нелинейных задач приводит к компактным (возможно, название не слишком удачное — термины "компакт" и "компактный" уже давно заняты в
математике под совсем другое!) разностным схемам — схемам повышенного
порядка точности на нерасширенном шаблоне. Действительно, и в элементарном
В случае, когда правая часть явно зависит от первой производной, либо не
получается компактной схемы, либо схема перестает быть экономичной. Действительно, чтобы не ухудшить порядок аппроксимации, необходимо вычислять значение первой производной в соответствующих узлах со вторым порядком. Если во всех точках использовать формулу с центральной разностью, то расширится шаблон схемы — в каждом шаге элементарных вычислений должно теперь участвовать пять точек. Если для точки с индексом n использовать формулу с центральной разностью, а для точек с индексами n + 1 и n - 1 — соответствующие формулы для односторонней производной, то для каждой точки шаблона придется вычислять заново значения функции f. При использовании классического вида аппроксимации Нумерова при каждом элементарном вычислении производится лишь однократное обращение к функции вычисления правой части — лишь для fn + 1, значения fn и fn - 1 уже получены при вычислениях для точки с индексом n - 1. Так как время вычислений, как правило, определяется в таких задачах именно количеством обращений к правой части, то вычисления замедляются в три раза.
y'' + p(x)y' + q(x)y = f(x)
где f, q, p — заданные функции,
используется разностная схема.
Предложить модификацию
Решение. Рассмотрим следующую краевую
y'' + p(x)y' + q(x)y = f(x)
с условиями периодичности
y(p)(0) = y(p)(1),
где p принимает значения 0 и 1. Отметим, что пока период решения считается известным. Отрезок [0, 1] возможно рассматривать без ограничения общности, так как любой отрезок можно перевести в единичный неособым линейным преобразованием, при этом вид уравнения существенно не изменится (функции p и q умножатся на постоянный множитель).
Краевая задача с условиями периодичности решается методом циклической
Очевидно, что оно аппроксимирует дифференциальную задачу со вторым порядком на равномерной сетке. На случай неравномерной сетки рассматриваемый метод легко обобщается.
В силу периодичности дискретное уравнение должно выполняться во всех точках
сетки, включая граничные. Кроме того, в силу граничного условия при p = 0, y0 = yN + 1, система сеточных соотношений примет вид:
где $$a_k = 1 - 0, 5p_k h, b_k = 2 - q_k h^2, c_k = 1 + 0, 5p_k h, \varphi_k = f_k h^2 .$$ Матрица системы линейных уравнений получается "почти трехдиагональной" — от трехдиагональной ее отличают всего два элемента в углах матрицы.
Обобщение стандартных прогоночных соотношений (для трехточечной прогонки) на периодический случай будет иметь вид
$$y_{n - 1} = \alpha_n y_n + \beta_n + \gamma y_n.$$Из приведенного выше соотношения для y0 сразу получаем, что $$\alpha_1 = {{c_0 }/ {b_0, }}$$ $$\beta_1 = - {{f_0 }/{b_0 }},$$ $$\gamma _{1} = a_{0}/ b_{0}.$$ Теперь несложно получить рекуррентную зависимость для прогоночных коэффициентов:
По приведенным выше формулам получаются значения коэффициентов для всех
уравнений с номерами меньше, чем n. Подставим теперь прогоночные
соотношения в последнее уравнение линейной системы. В итоге с учетом введенных выше обозначений получаем
а это соотношение, сгруппировав члены, можно переписать как:
$$y_{n} = \mu _{n}y_{0} + \eta _{n},$$где введены обозначения
$$$ \mu_n = \frac{- c_n}{a_n(\alpha_n + \gamma_n) - b_n}, \eta_n = \frac{\varphi_n - a_n \beta_n}{a_n (\alpha_n + \gamma_n) - b_n}. $$$Теперь выражение для значения сеточной функции yn - 1 подставляем в прогоночные соотношения. Получается выражение, связывающее yn - 1 с y0:
Отсюда получаем следующие рекуррентные соотношения:
$$\mu_{n - 1} = \alpha_n \mu_n + \gamma_n \mu_n , \eta_{n - 1} = \beta_n + \alpha_n \eta_n + \gamma_n \eta_n.$$Отметим, что эти коэффициенты вычисляются в обратном порядке — аналог обратного хода прогонки. Последнее соотношение приводит к явному выражению для y0. В результате получаем
Теперь информации для определения значений искомой функции во всех точках сетки (еще один ход прогонки) достаточно.
Алгоритм периодической прогонки был предложен А.А.Абрамовым.
$${y^{\prime\prime} = e^{y}, \quad y (0) = a, \quad y (1) = b, }$$
$${y^{\prime\prime} = - e^{y}, \quad y (0) = a , \quad y (1) = b.}$$a = 1 и сделав замену y'x = p(y).a = 1 и различных b. Что происходит при 0 < b < 1, 499719998? При b > 1,499719998? [10.6, c. 110].$$\varepsilon y^{\prime\prime} = (y^{\prime})^2, \quad y(0) = 1, y(0) = 0, 0 < \varepsilon \ll 1.$$
y'x = p(y).где $$q \in N$$ — натуральное число, $$0 < \varepsilon \ll 1.$$
Получить численное решение задачи в случаях:
u(x) = x2, A > 1, B> 1,u(x) = x2, A = 1, B > 1 или A > 1, B = 1,u(x) = | x |, A > 1, B > 1,u(x) = | x |, A = 1, B > 1.Что происходит с решением при увеличении q? (В численных
расчетах задать $$\varepsilon = 10^{ - 2}, 10^{ - 3}, 10^{ - 4}).$$
Теоретически задача исследована в [10.8, c. 170 - 171]. В случае u(x) = | x | появляется внутренний пограничный слой — узкая область в окрестности x = 0, где y отличны от | x |.
$$\begin{gather*} \varepsilon y^{\prime\prime} = y - y^3, \\ y(0) = A, \quad y(1) = B, \left|A\right| < \sqrt 2, \quad \left|B\right| < \sqrt 2, \quad 0 < \varepsilon \ll 1. \end{gather*}$$
Решением этой задачи являются так называемые пиковые, или пичковые,
структуры. В [10.8, с. 171 - 174] исследованы свойства решений и приведены графики восьми линейно независимых решений. Там же показано, что при фиксированных A и B существует по четыре линейно независимых решения, таких, что
при всех x, за исключением точек $$x_i = {i/n}, i = \overline{1, n -
1}, n \in N (n \ge 2),$$ где $$\lim\limits_{\varepsilon \to 0}y(x, \varepsilon ) = \sqrt 2 .$$
Найти численно такие структуры для n = 2, n = 3, выбирая
соответствующее начальное приближение при линеаризации по Ньютону. (Положить $$\varepsilon = 10^{ - 2}, 10^{ - 3}, 10^{ - 4}$$ ).
имеет решение с внутренним пограничным слоем в точке x = 1/2
[10.8, c. 175]. Исследовать, как его толщина зависит от параметра $$\varepsilon.$$
Какое начальное приближение надо использовать при решении задачи методом линеаризации?
Известно, что при A = 0, В = 0 данная краевая задача имеет еще два решения, кроме тривиального ( $$y \equiv 0$$ ). Найти их численно.
Всего у этой задачи счетное множество решений.
в зависимости от вида функции a(x). Рассмотреть поведение решения
при $$\varepsilon \to 0.$$ Удалось ли получить пограничный слой типа
всплеска?
y'' = - k2y, y (0) = y (1) = 0.
n — положительное целое число).k. С какими трудностями пришлось столкнуться? Как можно улучшить используемый алгоритм?y'(0) = y (1) = 0.$$\begin{gather*} \frac{d^2 y}{{dr}^2 } + \frac{1}{r}\frac{dy}{dr} = - k^2 y, \\ y(0) = 1, y(1) = 0 \end{gather*}$$
Собственными функциями этой задачи являются цилиндрические функции Бесселя нулевого порядка, а собственными значениями задачи будут нули этих функций:
k1 = 2, 404826 , k2 = 5, 520078 , k3 = 8, 653728 , k4 = 11, 791538.
Показать, что подстановка $$$ \tilde {y} = \frac{y}{\sqrt{r}} $$$ приводит уравнение к виду, для которого можно применять метод Нумерова. Решить численно
Найти численно все уровни энергии частицы в потенциальной яме с потенциалом $$U(x) = - 2{\mathop{\mathrm{sech}}\nolimits}^2 x$$ и соответствующие им функции распределения.
Указание: уровни энергии есть собственные значения $$\lambda _{k}$$ уравнения Шредингера $$y'' + (\lambda - U(x))y = 0$$ с условиями $$y(+ \infty ) = y(- \infty ) = 0,$$ а соответствующие им собственные функции и есть функции распределения.
Данный результат важен для решения уравнения Кортвега - Де Фриза, поэтому обсуждается в [10.9, c. 11 - 17].
x = 0. При этом собственные функции системы будут четными или нечетными относительно x = 0, причем четность или нечетность чередуется с ростом квантового числа (энергии). Проверить этот эффект численно. Каким способом в этом случае можно в два раза сократить объем вычислений при расчете собственных значений для потенциалов?Рассмотреть движение частицы в поле с потенциалом Тоды ([10.10, c. 87 - 89]):
$$$ \ddot {x} = 1 - e^{x} $$$Это уравнение можно трактовать как движение частицы в поле с потенциалом U(x) = ex - x.
Найти численно все периодические решения, удовлетворяющие следующим граничным условиям:
$$\begin{gather*} x(0) = x(120) = 0, \\ \dot {x}(0) = \dot {x}(120) = A \end{gather*}$$и дополнительному условию $$A \ge 10.$$
Как период колебаний зависит от A? Сколько решений получается?
Для получения официальных документов о завершении программы дополнительного профессионального образования (удостоверения о повышении квалификации, дипломов о профессиональной переподготовке и MBA) необходимо предоставить:
Внимание! Вы можете не заказывать доставку бумажной версии официального документы, а скачать его в электронном виде и распечатать самостоятельно. Информация о выданном документе в течение 1 месяца загружается в Федеральную информационную систему «Федеральный реестр сведений о документах об образовании и (или) о квалификации, документах об обучении» - ФИС ФРДО.
Доступ на новый сайт осуществляется с использованием адреса электронной почты, который был указан вами при регистрации на "старом". Мы постарались перенести все ваши данные с прежнего ресурса, однако не исключена вероятность потери части информации.
При возникновении проблемы со входом, воспользуйтесь функцией сброса пароля
Если вы обнаружите несоответствия, пожалуйста, сообщите нам.