Дискретизация Задачи Коши. Общая теория дискретизации ОДУ.#
Пусть заданы дифференциальная задача Коши (ОДУ):
Здесь \(\hat{L}\) - дифферециальный оператор. К примеру, \(\hat{L} = \frac{d}{dx}\) или \(\hat{L} = x \cdot \frac{d^2}{dx^2} -1 \).
Поставим ей в соответствие разностную задачу:
где \(y^{(h)}, F^{(h)}\) - сеточные функции, значения которых заданы в конечном наборе равноотстоящих (с шагом \(h\) ) точек \(\left\{x_{k}\right\}_{k=0}^{N}\). Заметим, что \(y^{(h)}\) мы заранее не знаем и определяем его из решения даной разностной задачи. Также отметим, что сеточная функция \(F^{(h)}\) не обязана быть \(F(x)\), спроецированной на сетку, - это необходимо только в пределе \(h \rightarrow 0\).
Соответственно, \(\hat{L}_{h}\) - оператор, действующий на множестве сеточных функций. Выбирается нами так, чтобы приблизить \(\hat{L}\). Для предыдущих примеров, соответственно, \(\hat{L}_h = \cfrac{y_{n+1}-y_n}{h}\) или \(\hat{L}_h = x_n \cdot \cfrac{y_{n+1} - 2 y_{n} + y_{n-1}}{h^2} - 1 \)
Введём проекцию истинного решения на сетку как
\(\left[y(x)\right]^{(h)}\) - тоже сеточная функция, полученная на основе значений истинного решения \(y(x)\) взятых точках области определения \(y^{(h)}\).
Примечание. Всё вышеперечисленное остаётся верным при трактовке \(y\) и \(f\) как векторов, т.е. рассмотрении систем дифференциальных уравнений.
Сходимость#
Введём невязку как сеточную функцию
Тогда будем говорить, что метод «сходится» при
При этом, если
то сходимость имеет порядок \(p\).
На практие её достаточно тяжело определить, т.к. нужно много знать про истинное решение \(y(x)\) - нужны некоторые хитрости. Одну из таких хитростей мы рассматривали на семинаре по численному интегрированию - оценка погрешности по Рунге.
Как оценивать порядок сходимости на практике#
Введём \(\hat{y}^{(h)}\) - сеточная функция, полученная на основе значений \(y^{(h)}\) взятых в каждой второй точке начиная с первой (размерности \(\hat{y}^{(h)}\) и \(y^{(2 h)}\) совпадают, и там в два раза меньше точек, чем в \(y^{(h)}\) - но вычисляются они по-разному!).
Тогда порядок сходимости можно явно оценить из следующего равенства:
Последнюю величину вполне себе можно считать в ходе расчёта и смотреть её ассимптотику от \(h\).
Доказательство:
Введём \(e^{(h)}\) и \(e^{(2 h)}\) - постоянные во всех точках сеточные функции, с единичной нормой. Пусть \(p\) - порядок сходимости. Тогда
Вычтем второе из третьего:
Заметим, что \(y^{(h)}-[y(x)]^{(h)} = C h^{p} e^{(h)}\).
А тогда, т.к. \(\left\|e^{(2 h)}\right\|=\left\|e^{(h)}\right\|=1\) получаем:
Аппроксимация#
Подействуем разностным оператором на проекцию истинного (!) решения \(\hat{L}_{h}[y(x)]^{(h)}\).
Посмотрим, решает ли эта проекция разностную задачу (если решает, то найденное в лоб решение разностной задачи является проекцией истинного решения в силу теоремы единственности).
Введём сеточную функцию
Если \(\psi^{(h)}\) не ноль, то \([y(x)]^{(h)}\) - не является решением разностной задачи \(\hat{L}_{h} y^{(h)}=F^{(h)}\). А она при любых \(h\) будет ненулевой :).
Поэтому говорят, что данная разностная задача аппроксимирует исходную задачу Коши, если
При этом, если \(\left\|\psi^{(h)}\right\|=O\left(h^p\right)\), то аппроксимация имеет порядок \(p\).
Примечание. Аппроксимация не эквивалентна сходимости! Нужен ещё один параметр.
Примечание. Для определения порядка аппроксимации как правило не нужно знать истинное решение. Пример ниже.
Устойчивость_#
Введём возмущение на свободный член \(\delta F^{(h)}\), и найдём решение возмущённой разностной задачи
Определение 1. Считаем, что задача устойчива к возмущению, при
Примечание. Норму можно брать любую на конечномерном векторном пространстве.
Такое определение устойчивости не очень удобно на практике, поэтому вводят более узкое определение.
Определение 2. Разностная схема \(\hat{L}_{h} y^{(h)}=F^{(h)}\) называется корректной, если
Заметим, что из корректности следует устойчивость.
Теорема Рябенького-Лакса#
Теорема. Решение разностной задачи сходится к решению задачи Коши тогда и только тогда, когда разностная задача устойчива и аппроксимирует задачу Коши. При этом порядок аппроксимации равен порядку сходимости.
Пример 1 исследования задачи#
Для решения задачи Коши
предложена разностная схема (заметим, что начальные условия также также дискретизовали! Пояснение ниже.)
Опеределить порядок аппроксимации и сравнить с порядком сходимости.
Примечание. Как дискретизовывать начальные (граничные) условия? В данном случае мы выбрали такую численную задачу, в которой для определения следующей точки \(y_{l+1}\) нам нужно две предыдущих: \(y_{l}, \; y_{l-1}\). Следовательно, чтобы «запустить» алгоритм, нам нужно уже иметь две точки: \(y_0\) и \(y_1\).
В «идеале» наше численное решение \(y^{(h)}\) должно быть равно проекции истинного решения \([y(x)]^{(h)}\). Давайте эти самые первые точки постулируем равными первым точкам сеточной функции \([y(x)]^{(h)}\):
Решение:
Аппроксимация (напоминаю про обозначение \([y]^{(h)} = [y(x)]^{(h)}\)):
Обратите внимание на первую и вторую точки. В ней чисто технически у нас другая численная схема, поэтому равенство именно такое.
Хотим понять ассимптотику этой сеточной функции при \(h \rightarrow 0\). Для этого воспользуемся свойствами проецирования и рядом Тейлора:
Подставляем всё в аппроксимацию:.
Упростим первые два слагаемых. Из начальных условий: \([y]_{0}=0,\left[y^{\prime}\right]_{0}=2 ;\) из уравнения: \(\left[y^{\prime \prime}\right]_{0}=-6\left[y^{\prime}\right]_{0}-5[y]_{0}=-12\).
В остальных слагаемых подставим выраженный из уравнения \(y'''=-6y''-5y'\).
Таким образом, \(\left\|\psi^{(h)}\right\|=\max _{l}\left|\psi_{l}^{(h)}\right| = O\left(h^{2}\right)\). Т.е. порядок аппроксимации равен двум. Заметим, что для анализа на аппроксимацию нам не понадобилось аналитически решать ДУ.
Более того, порядок аппроксимации можно было определить сходу - численные производные имели 2-й порядок аппроксимации, а начальные условия вообще 3-й.
Сходимость:
Аналитическое решение задачи Коши:
Решение разностной схемы также можно явно получить, т.к. тут простая реккурента:
Используя sympy, можно легко получить общее решение этой реккуренты и, соответственно, вычислить ассимптотику невязки:
from IPython.display import display, Math
from sympy import symbols, Function, rsolve, limit, Rational, exp
h = symbols('h')
l = symbols('l', integer=True)
y = Function('y')
recurrence = y(l+1) - ((2-5*h**2)/(1+3*h))*y(l) - ((3*h-1)/(1+3*h))*y(l-1) # Определяем реккуренту
solution = rsolve(recurrence, y(l), {y(0):0, y(1):2*h-6*h**2}) # Решаем реккуренту
display(Math("$y_l=$"))
display(solution)
delta = solution - Rational(1, 2)*(exp(-h*l) - exp(-5*h*l)) # Явным образом определяем невязку
display(Math("$\\delta_l = $"))
display(delta)
display(Math("$\\delta_l\\approx$"))
taylor_series = delta.series(h, 0) # Смотрим ассимптотику
display(taylor_series)
Видим первый ненулевой член невязки оказался ~ \(h^2\).
То есть, порядок сходимости оказался равен двум - такой же как и порядок аппроксимации, как и должно быть по теореме Рябенького-Лакса.
Отметим, что аналитически найти порядок сходимости можно примерно никогда, т.к.
Аналитическое решение может быть неизвестно
Практически никогда нельзя аналитически решить разностную задачу - для этого нам компьютор и нужон
Таким образом, если определить порядок аппроксимации можно всегда, то для определения сходимости нам нужно дополнительно научиться определять, устойчив ли метод.
Пример 2 - Аппроксимация начальных условий с произвольной точностью#
Решим ещё одну задачку.
Пусть имеем сложную задачу Коши:
Необходимо аппроксимировать начальные условия с произвольной заданной точностью.
Зачем? Нам может понадобится более чем две точки для старта алгоритма. Научимся сразу же приближать их с нужной точностью (порядком аппроксимации).
Сначала поймём, что любая производная в начальной точке \(y^{(k)}(1)\) нам известна - её легко можно выразить из уравнения, беря от него производные и подставляя предыдущие найденные значения. С помощью sympy вычислим первые 5 производных.
from IPython.display import display, Math, Latex
import sympy as sp
x = sp.symbols('x')
y = sp.Function('y')(x)
k = 5
# Define the differential equation
DE = y.diff(x, 2) + (1+x)/(1+x**2)*y.diff(x) + sp.cos(x)*y - 1/(1+x**2)
# Define the initial conditions
init_conds = {y.subs(x, 1): 3, y.diff(x).subs(x, 1): -1}
# Solve for the derivatives at x=1
for i in range(2, k):
DE_derivatived = DE.diff(x, i - 2).subs(x, 1)
deriv = y.diff(x, i).subs(x, 1)
sol = sp.solve(
DE_derivatived.subs(init_conds),
deriv,
list = True
)[0].simplify()
init_conds[deriv] = sol
# Print the results
print("Истинные начальные условия:")
for i in range(0, k):
deriv = y.diff(x, i)
_ = "$ y^{("+ str(i) + ")} (1)$ = $" + sp.pretty(sp.latex(init_conds[deriv.subs(x, 1)])) + "$"
display(Latex(_))
Истинные начальные условия:
Далее определим произвольное количество первых точек искомого решения \(y^{(h)}\) через Тейлора. Напомним, что для начальных условий в данном случае:
x = sp.symbols('x')
h = sp.symbols('h')
y = Function('y')
N = 6
print("Дискретизованные начальные условия с произвольной точностью:")
for i in range(len(init_conds.keys())):
y_curr = (y(1 + i*h).series(h, 0, 5).xreplace(init_conds)) # xreplace for reccurent expressions
_ = "$ y_{"+ str(i) + "}$ = $" + sp.pretty(sp.latex(y_curr), wrap_line = False) + "$"
display(Latex(_))
Дискретизованные начальные условия с произвольной точностью: