Для однозначного определения неизвестной функции \( u(x) \) уравнение (1) дополняется двумя граничными условиями на концах отрезка \( [0, l] \). Задаваться может функция \( u(x) \) (граничное условие первого рода), поток \( w(x) = −k(x) \frac{d u}{dx} (x) \) (граничное условие второго рода) или же их линейная комбинация (граничное условие третьего рода): $$ \begin{equation} \tag{2} u(0) = \mu_1, \quad u(l) = \mu_2, \end{equation} $$ $$ \begin{equation} \tag{3} −k(0) \frac{d u}{dx} (0) = \mu_1, \quad k(l) \frac{d u}{dx} (l) = \mu_2 \end{equation} $$ $$ \begin{equation} \tag{4} −k(0) \frac{d u}{dx} (0) + \sigma_1 u(0) = \mu_1, \quad k(l) \frac{d u}{dx} (l) + \sigma_2 u(l) = \mu_2. \end{equation} $$
Эллиптические уравнения второго порядка, прототипом которых является уравнение (1), используются при моделирование многих физико-механических процессов.
Кроме того,могут рассматриваться задачи с несамосопряженным оператором, когда, например, $$ \begin{equation} \tag{5} - \frac{d }{dx} \left( k(x) \frac{d u}{dx} \right) + v(x) \frac{d u}{dx} + q(x) u = f(x), \quad 0 < x < l. \end{equation} $$
Уравнение конвекции-диффузии (5) является модельным при исследование процессов в механике сплошной среды.
Обозначим через \( \bar{\omega}_h \) равномерную сетку с шагом \( h \) на отрезке \( [0, l] \): $$ \bar{\omega}_h = \{ x\ :\ x = x_i = ih, \ i = 0, 1, \ldots, N, \ Nh = l \}, $$ причем \( \omega_h \) — множество внутренних узлов сетки, \( \partial\omega_h \) — множество граничных узлов.
Будем использовать безындексные обозначения, когда \( u = u_i = u(x_i) \). Для левой разностной производной имеем $$ u_{\bar{x}} = \frac{u_{i} - u_{i-1}}{h} = \frac{d u}{dx}(x_i) - \frac{h}{2} \frac{d^2 u}{dx^2} + O(h^2), $$ т.е. левая разностная производная аппроксимирует производную \( \frac{d u}{dx} \) с первым порядком (погрешность аппроксимации \( O(h) \)) в каждом внутреннем узле при \( u(x) \in C^2(0, l) \). Аналогично для правой разностной производной $$ u_{x} = \frac{u_{i+1} - u_{i}}{h} = \frac{d u}{dx}(x_i) + \frac{h}{2} \frac{d^2 u}{dx^2} + O(h^2), $$ Для трехточечного шаблона \( (x_{i-1}, x_{i}, x_{i+1}) \) можно использовать центральную разностную производную: $$ u_{\mathring{x}} = \frac{u_{i+1} - u_{i-1}}{h} = \frac{d u}{dx}(x_i) + \frac{h^2}{3} \frac{d^3 u}{dx^3} + O(h^3), $$ которая аппроксимирует производную \( \frac{du}{dx} \) со вторым порядком при \( u(x) \in C^3(0, l) \).
Для второй производной \( \frac{d^2 u}{dx^2} \) получим $$ u_{x \bar{x}} = \frac{u_x - u_{\bar{x}}}{h} = \frac{u_{i+1} - 2 u_i + u_{i-1}}{h^2} $$
Этот разностный оператор апроксимирует в узле \( x = x_i \) вторую производную со вторым порядком при \( u(x) \in C^4 (0, l) \).
Для внутренних узлов сетки аппроксимируем дифференциальный оператор $$ \begin{equation} \tag{6} \mathcal{L} u = - \frac{d }{dx} \left( k(x) \frac{d u}{dx} \right) + q(x) u, \quad x \in (0, l), \end{equation} $$ с достаточно гладкими коэффициентами и решением разностным оператором $$ \begin{equation} \tag{7} Ly = - (a y_{\bar{x}})_x + cy, \quad x \in \omega_h. \end{equation} $$
Для аппроксимации со вторым порядком необходимо выбрать коэффициенты разностного оператора так, чтобы $$ \begin{align} \tag{8} \frac{a_{i+1} - a_i}{h} &= \frac{d k}{dx}(x_i) + O(h^2),\\ \tag{9} \frac{a_{i+1} + a_i}{2} &= k(x_i) + O(h^2),\\ \tag{10} c_i &= q(x_i) + O(h^2). \end{align} $$
В соответствии с (10) положим, например, \( c_i = q(x_i) \), а условиям (8) и (9) удовлетворяют, в частности, следующие формулы для определения \( a_i \): $$ \begin{align*} a_{i} &= k_{i-1/2} = k(x_i - 0.5h),\\ a_{i} &= \frac{k_{i} + k_{i-1}}{2}, \\ a_{i} &= 2 \left( \frac{1}{k_i} + \frac{1}{k_{i-1}} \right)^{-1}. \end{align*} $$
Метод формальной замены дифференциальных операторов разностными может использоваться и при аппроксимации граничных условий.
Наиболее просто аппроксимируются граничные условия (2): $$ \begin{equation} \tag{11} y_0 = \mu_1, \quad y_N = \mu_2. \end{equation} $$
Для аппроксимации граничных условий второго и третьего рода со вторым порядком в граничных узлах \( x = x_0 = 0 \) и \( x = x_N = l \) привлекается уравнение (1) — аппроксимация на решениях задачи. В случае уравнения (1) краевые условия (4) аппроксимируются разностными соотношениями: $$ \begin{align*} -a_1 y_{\bar{x},1} + \left( \sigma_1 + \frac{h}{2} q_0 \right) u_0 &= \mu_1 + \frac{h}{2} f_0 \\ a_N y_{\bar{x},N} + \left( \sigma_2 + \frac{h}{2} q_N \right) u_N &= \mu_2 + \frac{h}{2} f_N \end{align*} $$
Итак, мы имеем разностную схему, состоящую из разностного уравнения $$ \begin{equation} \tag{12} - (a y_{\bar{x}})_x + cy = \varphi, \quad x \in \omega_h, \end{equation} $$ дополненного аппроксимацией соответствующих граничных условий.
Для нахождения приближенного решения краевой задачи для обыкновенного дифференциального уравнения необходимо решить соответствующую систему линейных алгебраических уравнений. Для нахождения разностного решения используются традиционные прямые методы линейной алгебры. Излагаемый метод прогонки (алгоритм Томаса), как хорошо известно, является классическим методом Гаусса для матриц специальной ленточной структуры.
Система уравнений (12) с граничными условиями представляет собой частный случай систем линейных алгебраических уравнений с трехдиагональной матрицей, которые имеют вид $$ \begin{equation} \tag{13} - \alpha_{i} y_{i-1} + \gamma_{i} y_{i} - \beta_{i} y_{i-1} = \varphi_{i}, \quad i = 1, 2, \ldots, N-1, \end{equation} $$ $$ \begin{equation} \tag{14} y_0 = \kappa_1 y_1 + \mu_1, \quad y_N = \kappa_2 y_{N-1} + \mu_2. \end{equation} $$ Для численного решения таких систем применяется метод прогонки, представляющий собой рекуррентные формулы для вычисления прогоночных коэффициентов (прямая прогонка): $$ \xi_{i+1} = \frac{\beta_i}{\gamma_i - \alpha_i \xi_i},\quad \vartheta_{i+1} = \frac{\varphi_i + \alpha_i\vartheta_{i}}{\gamma_i - \alpha_i \xi_i}, \quad i = 1, 2, \ldots, N-1, $$ при $$ \xi_1 = \kappa_1, \quad \vartheta_1 = \mu_1 $$
Для решения имеем (обратная прогонка) $$ y_i = \xi_{i+1} y_{i+1} + \vartheta_{i+1}, \quad i = N-1, N-2, \ldots, 0, \quad y_N = \frac{\kappa_2 \vartheta_N + \mu_2}{1 - \kappa_2 \xi_N}. $$
Алгоритм прогонки можно применять, если знаменатели в расчетных формулах не обращаются в нуль. Для возможности применения метода прогонки достаточно потребовать, чтобы для системы уравнений (13), (14) удовлетворяли условиям $$ \alpha_i \ne 0, \quad \beta_i \ne 0, \quad |c_i| \geq |\alpha_i| + |\beta_i|, \quad i = 1, 2, \ldots, N-1, $$ $$ |\kappa_1| \leq 1, \quad |\kappa_2| \leq 1. $$
Написать программу, которая реализует этот метод стрельбы при решении краевой задачи для нелинейного уравнения второго порядка $$ \frac{d^2 u}{dx^2} = f\left( x, u, \frac{du}{dx} \right), \quad 0 < x < l, $$ $$ u(0) = \mu_1, \quad u(l) = \mu_2. $$
Продемонстрировать работоспособность этой программы при решении краевой задачи $$ \frac{d^2 u}{dx^2} = 100 u(u-1), \quad 0 < x < 1, $$ $$ u(0) = 0, \quad u(1) = 2, $$ на различных сетках.
Метод стрельбы основан на решении задач Коши для системы уравнений первого порядка $$ \frac{du}{dx} = v, \quad \frac{dv}{dx} = f\left( x, u, \frac{du}{dx} \right), \quad 0 < x < l $$ с начальными условиями $$ u(0) = \mu_1, \quad v(0) = \theta. $$ Величина \( \theta \) подбирается так, чтобы выполнялось граничное условие \( u(l) = \mu_2 \). Для нахождения параметра пристрелки \( \theta \) можно использовать метод Ньютона для решения соответствующего нелинейного уравнения, а для решения задачи Коши для системы обыкновенных дифференциальных уравнений — метод явный метод Рунге-Кутта четвертого порядка точности.
Написать программу для приближенного решения задачи $$ \frac{d^2 u}{dx^2} = f(x, u), \quad 0 < x < l, $$ $$ u(0) = \mu_1, \quad u(l) = \mu_2 $$ разностным методом на равномерной сетке, используя итерационное решение сеточной задачи методом Ньютона. Использовать эту программу для решения задачи с правой частью \( f(t, u) = - e^{u} \), \( l = 1 \) и граничных условиях \( \mu_1 = \mu_2 = 0 \).
Разностная схема $$ - (a y_{\bar{x}})_x = \varphi, \quad x \in \omega_h, $$ в которой $$ a_i = \left( \frac{1}{h} \int_{x_{i-1}}^{x_{i}} \frac{dx}{k(x)} \right)^{-1}, $$ $$ \varphi_i = \frac{1}{h^2} \left( a_{i+1} \int_{x_i}^{x_{i+1}} \frac{1}{k(x)} \left(\int_{x_i}^{x} f(s) ds\right) dx + a_{i} \int_{x_{i-1}}^{x_{i}} \frac{1}{k(x)} \left(\int_{x}^{x_i} f(s) ds\right) dx \right) $$ для приближенного решения задачи $$ - \frac{d}{dx} \left( k(x) \frac{du}{dx} \right) = f(x), \quad 0 < x < l, $$ $$ u(0) = \mu_1, \quad u(l) = \mu_2 $$ является точной. Написать программу разностного решения краевых задач с использованием точной разностной схемы при использовании квадратурных формул трапеций для вычисления коэффициентов. Применить эту программу для решения краевой задачи при \( l = 1 \) и \( \mu_1 = \mu_2 = 0 \) с разрывным коэффициентом $$ k(x) = \begin{cases} 1, & 0 < x < 0.5,\\ 100, & 0.5 < x < 1 \end{cases} $$ и правой частью \( f(x) = 0 \). Сравнить приближенное решение на последовательности трех вложенных друг в сетках \( h = h_0, h_0/2, h_0/4 \).