Линейная краевая задача

Формулировка

Пусть имеем следующую ОДУ 1-го порядка:

\[ \frac{d \mathbf{y}}{d x}=\hat{\mathbf{A}}(x) \cdot \mathbf{y}+\mathbf{f}(x) \]

где \(x \in[0, X], \; \mathbf{y}, \mathbf{f}-n\)-мерные векторы, \(\; \hat{\mathbf{A}}(x) \in \mathbb{R}^{n \times n}\)

Для существования единственного решения необходимо также задать \(n\) граничных условий:

\[ \hat{\mathbf{L}} \cdot \mathbf{y}(0)+\hat{\mathbf{R}} \cdot \mathbf{y}(X)=\mathbf{g}, \quad \hat{\mathbf{L}}, \hat{\mathbf{R}} \in \mathbb{R}^{n \times n}, \quad \mathbf{g} \in \mathbb{R}^n \]

По сути, здесь просто записано \(n\) уравнений на граничные точки в хитром виде. Как правило, условия пишутся гораздо проще.

Пример. Сведём следующее ДУ к виду выше

\[\begin{split} \begin{aligned} & y^{\prime \prime}+q(x) y^{\prime}-p(x) y=f(x), \quad 0<x \leq 1 \\ & y(0)=Y_0, \quad y(1)=Y_1 \end{aligned} \end{split}\]

Для того чтобы свести дифференциальное уравнение к матричному виду, мы можем использовать следующее представление:

\[\begin{split} \mathbf{y}^{\prime}=\left(\begin{array}{c} y^{\prime} \\ y^{\prime \prime} \end{array}\right), \quad \mathbf{y}=\left(\begin{array}{c} y \\ y^{\prime} \end{array}\right), \quad \hat{\mathbf{A}}(x)=\left(\begin{array}{cc} 0 & 1 \\ -p(x) & q(x) \end{array}\right), \quad \mathbf{f}(x)=\left(\begin{array}{c} 0 \\ f(x) \end{array}\right) \end{split}\]

Тогда исходное уравнение можно записать в виде

\[\begin{split} \mathbf{y}^{\prime}=\hat{\mathbf{A}}(x) \mathbf{y}+\mathbf{f}(x), \quad \mathbf{y}(0)=\left(\begin{array}{c} Y_0 \\ ? \end{array}\right), \quad \mathbf{y}(1)=\left(\begin{array}{c} Y_1 \\ ? \end{array}\right) \end{split}\]

Видим, что у нас есть два краевых условия. Т.е. задача корректно поставлена.

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

Обусловленность краевых условий

C точки зрения вычислений также важно, чтобы матрица была хорошо обусловленной. Рассмотрим случай, когда матрица \(\hat{\mathbf{A}}(x)\) имеет \(n\) действительных собственных значений, из них:

1. \(n_1\) очень больших положительных: \(\max \lambda_i \cdot L \gg 1\);

2. \(n_2\) очень больших отрицательных: \(\min \lambda_i \cdot L \ll -1\);

3. Все остальные умеренной величины \(\left|\lambda_i\right| \cdot L \simeq 1\).

Это означает, что точное решение представляется в виде \(n_1\) сильно затухающих экспонент и \(n_2\) сильно возрастающих экспонент.

Для корректности разрешения системы уравнений краевых условий необходимо, чтобы не менее \(n_1\) условий были поставлены слева и не менее \(n_2\) условий справа. Это обусловлено тем, что наши методы, как правило, неустойчивы для быстровозрастающих экспонент, но если их начальные условия заданы «справа», то идя влево они как бы убывают.

Для слабо возрастающих или слабо убывающих компонент краевые условия могут быть поставлены на любом из концов отрезка.

Методы решения линейной краевой задачи

Фундаментальная система решений 🥇

Т.к. задача линейная

\[ \frac{d \mathbf{y}}{d x}=\hat{\mathbf{A}}(x) \cdot \mathbf{y}+\mathbf{f}(x) \]

\[ \hat{\mathbf{L}} \cdot \mathbf{y}(0)+\hat{\mathbf{R}} \cdot \mathbf{y}(X)=\mathbf{g} \]

То общее решение является суммой частного неоднородного и общего однородного (ФСР):

\[ \mathbf{y}(x)=\overline{\mathbf{y}}(x)+\sum_{k=1}^n \alpha_k \mathbf{y}^k \]

Здесь \(\mathbf{y}^k(x)\) есть полная фундаментальная система решений однородной задачи с заданными нами начальными данными, например:

\[ \frac{d \mathbf{y}^k}{d x}=\mathbf{A}(x) \cdot \mathbf{y}^k, \quad k=1,2, \ldots, n \]

\[ \mathbf{y}^k(0)=(0, \ldots, 0,1,0, \ldots, 0)^{\mathrm{T}} \]

Т.е. для их нахождения нужно решить \(n\) задач Коши, которые мы уже умеем решать. Например, Рунге-Куттой.

Выберем также такое частное решение, что оно удовлетворяет задаче Коши с нулевыми начальными условиями:

\[ \frac{d \overline{\mathbf{y}}}{d x}=\mathbf{A} \overline{\mathbf{y}}+\mathbf{f}(x) \]

\[ \overline{\mathbf{y}}(0) = \vec{0} \]

Осталось лишь подобрать коэффициенты \(\alpha_k\), чтобы полученное нами общее решение удовлетворяло краевым условиям. Просто подставим наше общее решение в уравнение на краевые условия:

\[ \mathbf{R}\left[\sum_{k=1}^n \alpha_k \mathbf{y}^k(0)\right]+\mathbf{L} \overline{\mathbf{y}}(X)+\mathbf{L}\left[\sum_{k=1}^n \alpha_k \mathbf{y}^k(X)\right]=\mathbf{g} \]

Для однозначной разрешимости этого уравнения при произвольных краевых условиях, достаточно линейной независимости векторов при варьируемых \(\alpha_k\):

\[\begin{split} \operatorname{det}\left(\begin{array}{c} \mathbf{R y}^1(0)+\mathbf{S y ^ { 1 } ( L )} \\ \cdots \\ \mathbf{R y}^n(0)+\mathbf{S y ^ { n } ( L )} \end{array}\right) \neq 0 \end{split}\]

Здесь каждый вектор в строке необходимо дополнительно транспонировать, чтобы получилась матрица.

Это условие не обязано всегда выполняться, то есть даже в линейном случае могут встретиться такие постановки задачи, что решение либо не существует, либо не единственное, если матрица СЛАУ для коэффициентов вырожденная.

Более того, если в спектре матрицы \(\hat{\mathbf{A}}(x)\) есть большие положительные собственные числа, мы не сможем решить задачу Коши из-за неустойчивости. Поэтому вполне себе можно находить общие решения однородного решая краевую задачу для них каким-нибудь другим способом (см. ниже).

Таким образом, метод ФСР прост и хорош, но есть нюансы.

Метод дифференциальной прогонки 🥇

Проиллюстрируем этот метод на конкретном примере. Все выкладки легко перенести на произвольный линейный случай, но почему-то в литературе никто этого не сделал ;\

Имеем уравнение из примера:

\[\begin{split} \begin{aligned} & y^{\prime \prime}+q(x) y^{\prime}-p(x) y=f(x), \quad 0<x \leq 1 \\ & y(0)=Y_0, \quad y(1)=Y_1 \end{aligned} \end{split}\]

Перепишем это уравнение в виде следующей системы (естественно, это можно сделать не единственным способом)

\[\begin{split} \left\{\begin{aligned} & \frac{d y}{d x}=a(x) v+b(x) \\ & \frac{d v}{d x}=c(x) y+g(x) \end{aligned} \right . \end{split}\]

И будем искать решение в виде т.н. прогоночного соотношения:

\[ y(x)=\alpha(x) v(x)+\beta(x) \]

где \(\alpha(x)\) и \(\beta(x)\) - пока неизвестные функции (прогоночные коэффициенты), для которых необходимо получить дифференциальные уравнения.

Тогда, решение можно найти, решив следующие 3 задачи Коши в соответствующем порядке:

\[\begin{split} \begin{aligned} & \alpha^{\prime}(x)+c \cdot \alpha^2(x)-a=0, \quad \alpha(0)=0 \\ \\ & \beta^{\prime}(x)+\alpha \cdot \beta (x) \cdot c+\alpha \cdot g-b=0, \quad \beta(0)=Y_0\\ \\ & v^{\prime}(x)=\alpha \cdot c \cdot v(x)+c \cdot \beta+g, \quad v(1)=\cfrac{Y_1 - \beta(1)}{\alpha(1)} \end{aligned} \end{split}\]

Метод очень хорош для нежестких задач.

Доказательство

Продифференцируем это соотношение

\[ y^{\prime}=\alpha^{\prime} v+\alpha v^{\prime}+\beta^{\prime} \]

и подставим в него уравнения системы

\[\begin{split} \left\{\begin{aligned} & y^{\prime}=a v+b \\ & v^{\prime}=c y+g \end{aligned} \right . \end{split}\]

В результате получим, что \(a v+b=\alpha^{\prime} v+\alpha(c y+g)+\beta^{\prime}\). Подставим в полученное равенство прогоночное соотношение \(y=\alpha v+\beta\).

Тогда

\[ a v+b=\alpha^{\prime} v+\alpha^2 c v+c \beta+\alpha g+\beta^{\prime} v+\beta . \]

После приведения подобных членов имеем равенство

\[ v\left(\alpha^{\prime}+c \alpha^2-a\right)+\left(\beta^{\prime}+c \beta \alpha+\alpha g-b\right)=0 . \]

Приравнивая к нулю коэффициенты при \(v\) и единице, получим два дифференциальных уравнения для прогоночных коэффициентов:

\[\begin{split} \left \{\begin{aligned} & \alpha^{\prime}+c \alpha^2-a=0 \\ & \beta^{\prime}+\alpha \beta c+\alpha g-b=0 \end{aligned} \right . \end{split}\]

Дополним их начальными условиями.

• Условия слева

Левое краевое условие вида \(y(0)=Y_0\) запишем в виде прогоночного соотношения \(y(0)=\alpha(0) v(0)+\beta(0)\), полагая \(\alpha(0)=0, \quad \beta(0)=Y_0\).

Таким образом, получаем начальные данные для двух задач Коши для \(\alpha(x)\) и \(\beta(x)\), которые могут быть решены численно.

• Условия справа

Их просто подгоним из равенства граничных условий справа.

Далее воспользуемся уравнением \(v^{\prime}=c y+g\), подставим в него прогоночное соотношение \(y=\alpha v+\beta\), получим дифференциальное уравнение для \(v:\)

\[ v^{\prime}(x)=\alpha \cdot c \cdot v(x)+c \cdot \beta+g, \quad v(1)=\cfrac{Y_1 - \beta(1)}{\alpha(1)} \]

Метод линейной интерполяции*

Рассмотрим смешанную краевую задачу для уравнения

\[ y^{\prime \prime}+p(x) y^{\prime}+q(x) y+f(x)=0, \quad 0<x<l \]

Граничные условия возьмем в смешанном виде

\[\begin{split} \begin{gathered} y(0)=a, \\ y^{\prime}(l)+\alpha y(l)=b . \end{gathered} \end{split}\]

Выберем любые два значения \(\lambda=\lambda_{1}, \;\lambda=\lambda_{2}\) и решим две задачи Коши с начальными условиями:

\[\begin{split} \begin{aligned} &y_{1}(0)=a, \quad y_{1}^{\prime}(0)=\lambda_{1} \\ &y_{2}(0)=a, \quad y_{2}^{\prime}(0)=\lambda_{2} \end{aligned} \end{split}\]

Полученные решения обозначим как \(y=y_{1}(x)\) и \(y=y_{2}(x)\). Найдем соответствующие значения левых частей в граничном условии.

\[\begin{split} \begin{aligned} &y_{1}^{\prime}(l)+\alpha y_{1}(l)=b_{1} \\ &y_{2}^{\prime}(l)+\alpha y_{2}(l)=b_{2} \end{aligned} \end{split}\]

Теперь значение недостающего начального условия \(y^{\prime}(0)=\lambda\) можно найти с помощью линейной интерполяции:

\[ \frac{\lambda-\lambda_{1}}{\lambda_{2}-\lambda_{1}}=\frac{b-b_{1}}{b_{2}-b_{1}} \quad \text {, T.e. } \lambda=\lambda_{1}+\left(\lambda_{2}-\lambda_{1}\right) \frac{b-b_{1}}{b_{2}-b_{1}} . \]

Полученное значение \(\lambda\) и будет являться недостающим начальным условием. Объясняется это линейностью задачи. Как известно, дифференциальное уравнение имеет общее решение

\[ y(x)=c_{1} u_{1}(x)+c_{2} u_{2}(x)+y_{n}(x), \]

где \(u_{1}(x)\) и \(u_{2}(x)\) - линейно-независимые решения однородного уравнения (при \(f(x) \equiv 0)\), а \(y_{n}(x)\) - какое-либо решение неоднородного уравнения (частное решение неоднородного уравнения). Удовлетворяя левому граничному условию, в общем решении останется одна неизвестная постоянная, которая входит в выражение для \(y(x)\) линейным образом. Проведя в плоскости \((\lambda, b)\) прямую, проходящую через две точки \(\left(\lambda_{1}, b_{1}\right)\) и \(\left(\lambda_{2}, b_{2}\right)\) при заданном значении \(b\) мы однозначно найдем точное значение \(\lambda\).

Теперь таблицу значений функции \(y(x)\) (и ее производной) можно найти интерполяцией \(y(x)=y_{1}(x)+\left[y_{1}(x)-y_{1}(x)\right] \cdot \frac{\lambda-\lambda_{1}}{\lambda_{2}-\lambda_{1}}\).

Однако на практике, жертвуя машинным временем, обычно проводят третий расчет задачи Коши с условиями

\[ y(0)=a, \quad y^{\prime}(0)=\lambda . \]

Замечание.

На практике обычно выбирают простейшие значения \(\lambda\), например \(\lambda_{1}=1\) и \(\lambda_{2}=0\). При этом, если само уравнение является однородным, т.е. \(f(x) \equiv 0\), и граничное условие так же однородное, \(y(0)=0\), то решение имеет вид \(y_{2}(x) \equiv 0, b_{2}=0\). Тогда второй расчет (при \(\lambda_{2}=0\) ) нет необходимости производить и формула даст ответ в виде

\[ \lambda=\lambda_{1} \frac{b}{b_{1}} . \]

Поэтому, если предложено решить однородное уравнение, то следует посмотреть, есть ли однородное граничное условие и если оно есть, то начинать расчет следует от этой границы! Эта рекомендация остается в силе и для других методов решения краевых задач.

Жесткость в краевых задачах

Рассмотрим краевую задачу:

\[ \frac{d \mathbf{y}}{d x}=\mathbf{A}(x) \cdot \mathbf{y}+\mathbf{f}, \quad x \in[a, b] \]

с граничными условиями в следующем виде:

\[ \hat{\mathbf{L}} \cdot \mathbf{y}(0)+\hat{\mathbf{R}} \cdot \mathbf{y}(X)=\mathbf{g}, \quad \hat{\mathbf{L}}, \hat{\mathbf{R}} \in \mathbb{R}^{n \times n}, \quad \mathbf{g} \in \mathbb{R}^n \]

Определение 1. Краевая задача является жесткой, если спектр собственных значений матрицы \(\hat{\mathbf{A}}(x)\) можно разделить на три части:

1. Левый жесткий спектр, для которого выполняется \(\operatorname{Re} \Lambda_i^1 \leq-\Lambda_0\), \(\left|\operatorname{Im} \Lambda_i^1\right|<\Lambda_0, \Lambda_0(b-a) \gg 1, \Lambda_0 \gg 1, \quad i=1, \ldots, N_1\)

2. Правый жесткий спектр, для которого \(\operatorname{Re} \Lambda_i^2 \geq \Lambda_0, \quad\left|\operatorname{Im} \Lambda_i^2\right|<\Lambda_0\), \(i=N_1+1, \ldots, N_2\)

3. Мягкий спектр \(\left|\lambda_i\right| \leq \lambda_0, \quad \lambda_0(b-a) \simeq O(1), i=N_2+1, \ldots, N\).

Отношение \(\Lambda_0 / \lambda_0 \sim 1\) называется параметром, характеризующим жесткость системы.

Общее решение жесткой линейной системы ОДУ с постоянными коэффициентами имеет вид

\[ \mathbf{y}(x)=\sum_{i=1}^{N_1} \gamma_i^1 e^{\Lambda_i^1 x} \boldsymbol{\Omega}_i^1+\sum_{i=N_1+1}^{N_2} \gamma_i^2 e^{\Lambda_i^2 x} \boldsymbol{\Omega}_i^2+\sum_{i=N_2+1}^N \gamma_i^3 e^{\lambda_i x} \boldsymbol{\omega}_i \]

где \(\boldsymbol{\Omega}_i^1, \boldsymbol{\Omega}_i^2, \boldsymbol{\omega}_i-\) собственные векторы матрицы \(\mathbf{A}\), соответствующие трем частям спектра.

Понятна качественная структура этого решения, содержащего и левый, и правый пограничные слои. Для решений жестких линейных краевых задач должно быть выполнено неравенство \(\|\mathbf{y}\| \leq \mathrm{C}(\|\mathbf{f}\|+\|\mathbf{g}\|)\), \(\|\mathbf{f}\|\) и \(\|\mathbf{g}\|\) - нормы правых частей в системе ОДУ и краевых условиях соответственно.

Правильно работать только в классе вычислительно корректных задач, для которых \(C=O(1)\ll\exp \left(\Lambda_0(b-a)\right)\).

Для вычислительной корректности (почти) достаточно, чтобы число краевых условий на левом конце отрезка интегрирования не было меньше быстро убывающих вправо решений, на правом конце - не меньше числа быстро убывающих решений влево.

В противном случае краевая задача оказывается вычислительно некорректной, так как \(C=O\left(\exp \left(\Lambda_0(b-a)\right)\right)\).