Метод Кранка – Николсона

В численном анализе метод Кранка-Николсона представляет собой метод конечных разностей, используемый для численного решения уравнения теплопроводности и аналогичных уравнений в частных производных . Это метод второго порядка по времени. Он неявен во времени, может быть записан как неявный метод Рунге-Кутты и численно устойчив . Метод был разработан Джоном Крэнком и Филлис Николсон в 1940-х годах.

Для уравнений диффузии (и многих других уравнений) можно показать, что метод Кранка – Николсона безоговорочно устойчив . Однако приближенные решения все же могут содержать (затухающие) паразитные колебания, если отношение шага по времени $\Delta t$ умноженное на коэффициент температуропроводности на квадрат пространственного шага, $\Delta x^{2}$, велико (обычно больше 1/2 согласно анализу устойчивости фон Неймана ). По этой причине, когда необходимы большие временные шаги или высокое пространственное разрешение, менее точный обратный метод Эйлера , который одновременно стабилен и невосприимчив к колебаниям. часто используется

Принцип

Кранка – Николсона Трафарет для одномерной задачи.

Метод Кранка-Николсона основан на правиле трапеций , дающем сходимость второго порядка по времени. Для линейных уравнений правило трапеций эквивалентно неявному методу средней точки. - простейший пример неявного метода Рунге-Кутты Гаусса-Лежандра , который также обладает свойством быть геометрическим интегратором . Например, предположим, что в одном измерении уравнение в частных производных имеет вид $${\frac {\partial u}{\partial t}}=F\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right).$$

Сдача в аренду $u(i\Delta x,n\Delta t)=u_{i}^{n}$ и $F_{i}^{n}=F$ оценено за $i,n$ и $u_{i}^{n}$уравнение для метода Кранка – Николсона представляет собой комбинацию прямого метода Эйлера при $n$ и обратный метод Эйлера при $n+1$ (однако обратите внимание, что сам метод не является просто средним значением этих двух методов, поскольку обратное уравнение Эйлера имеет неявную зависимость от решения):

${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)$	вперед Эйлер
${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}=F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)$	обратный Эйлер
${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {1}{2}}\left[F_{i}^{n+1}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)+F_{i}^{n}\left(u,x,t,{\frac {\partial u}{\partial x}},{\frac {\partial ^{2}u}{\partial x^{2}}}\right)\right]$	Кренк – Николсон

Обратите внимание, что это неявный метод : чтобы получить «следующее» значение $u$ со временем необходимо решить систему алгебраических уравнений. Если уравнение в частных производных нелинейно, дискретизация также будет нелинейной, так что продвижение во времени потребует решения системы нелинейных алгебраических уравнений, хотя линеаризация возможна. Во многих задачах, особенно линейной диффузии, алгебраическая задача является трехдиагональной и может быть эффективно решена с помощью алгоритма трехдиагональной матрицы , который дает быстрое решение. ${\mathcal {O}}(N)$ прямое решение, в отличие от обычного ${\mathcal {O}}(N^{3})$ для полной матрицы, в которой $N$ указывает размер матрицы.

Пример: 1D-диффузия

Метод Кранка-Николсона часто применяется к задачам диффузии . Например, для линейной диффузии:

${\frac {\partial u}{\partial t}}=a{\frac {\partial ^{2}u}{\partial x^{2}}},$

применяя пространственную дискретизацию с конечной разностью для правой части, тогда дискретизация Кранка – Николсона будет равна

${\frac {u_{i}^{n+1}-u_{i}^{n}}{\Delta t}}={\frac {a}{2(\Delta x)^{2}}}\left((u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1})+(u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n})\right)$

или, позволяя $r={\frac {a\Delta t}{2(\Delta x)^{2}}}$,

$-ru_{i+1}^{n+1}+(1+2r)u_{i}^{n+1}-ru_{i-1}^{n+1}=ru_{i+1}^{n}+(1-2r)u_{i}^{n}+ru_{i-1}^{n}.$

Учитывая, что члены в правой части уравнения известны, это трехдиагональная задача, так что $u_{i}^{n+1}$ может быть эффективно решена с использованием алгоритма трехдиагональной матрицы вместо гораздо более дорогостоящего обращения матрицы .

Квазилинейное уравнение, например (это минималистичный пример, а не общий)

${\frac {\partial u}{\partial t}}=a(u){\frac {\partial ^{2}u}{\partial x^{2}}},$

привело бы к нелинейной системе алгебраических уравнений, которую невозможно было бы легко решить, как указано выше; однако в некоторых случаях можно линеаризовать задачу, используя старое значение для $a$, то есть, $a_{i}^{n}(u)$ вместо $a_{i}^{n+1}(u)$. В других случаях можно оценить $a_{i}^{n+1}(u)$ используя явный метод и сохраняя стабильность.

Пример: 1D-диффузия с адвекцией для устойчивого потока с несколькими соединениями каналов.

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

Здесь мы моделируем концентрацию растворенного загрязняющего вещества в воде. Эта задача состоит из трех частей: известного уравнения диффузии ( $D_{x}$ выбрана постоянной), адвективная составляющая (что означает, что система развивается в пространстве за счет поля скоростей), которую мы выбираем постоянной $U_{x}$и латеральное взаимодействие продольных каналов ( $k$):

${\frac {\partial C}{\partial t}}=D_{x}{\frac {\partial ^{2}C}{\partial x^{2}}}-U_{x}{\frac {\partial C}{\partial x}}-k(C-C_{N})-k(C-C_{M}),$ ( 1 )

где $C$ — концентрация загрязняющего вещества, а нижние индексы $N$ и $M$ соответствуют предыдущему и следующему каналу.

Метод Кранка–Николсона (где $i$ представляет позицию, и $j$ time) преобразует каждый компонент PDE в следующее:

${\frac {\partial C}{\partial t}}\Rightarrow {\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}},$ ( 2 )

${\frac {\partial ^{2}C}{\partial x^{2}}}\Rightarrow {\frac {1}{2(\Delta x)^{2}}}\left((C_{i+1}^{j+1}-2C_{i}^{j+1}+C_{i-1}^{j+1})+(C_{i+1}^{j}-2C_{i}^{j}+C_{i-1}^{j})\right),$ ( 3 )

${\frac {\partial C}{\partial x}}\Rightarrow {\frac {1}{2}}\left({\frac {(C_{i+1}^{j+1}-C_{i-1}^{j+1})}{2(\Delta x)}}+{\frac {(C_{i+1}^{j}-C_{i-1}^{j})}{2(\Delta x)}}\right),$ ( 4 )

$C\Rightarrow {\frac {1}{2}}(C_{i}^{j+1}+C_{i}^{j}),$ ( 5 )

$C_{N}\Rightarrow {\frac {1}{2}}(C_{Ni}^{j+1}+C_{Ni}^{j}),$ ( 6 )

$C_{M}\Rightarrow {\frac {1}{2}}(C_{Mi}^{j+1}+C_{Mi}^{j}).$ ( 7 )

Теперь мы создадим следующие константы для упрощения алгебры:

$\lambda ={\frac {D_{x}\,\Delta t}{2\,\Delta x^{2}}},$

$\alpha ={\frac {U_{x}\,\Delta t}{4\,\Delta x}},$

$\beta ={\frac {k\,\Delta t}{2}},$

и подставим ( 2 ), ( 3 ), ( 4 ), ( 5 ), ( 6 ), ( 7 ), $\alpha$, $\beta$ и $\lambda$ в ( 1 ). Затем мы помещаем новые временные условия слева ( $j+1$) и условия настоящего времени справа ( $j$) получить

$-\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +2\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}={}$

$\qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -2\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.$

Чтобы смоделировать первый канал, мы понимаем, что он может контактировать только со следующим каналом ( $M$), поэтому выражение упрощается до

$-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}-\beta C_{Mi}^{j+1}={}$

$\qquad {}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.$

Точно так же, моделируя последний канал, мы понимаем, что он может контактировать только с предыдущим каналом ( $N$), поэтому выражение упрощается до

$-\beta C_{Ni}^{j+1}-(\lambda +\alpha )C_{i-1}^{j+1}+(1+2\lambda +\beta )C_{i}^{j+1}-(\lambda -\alpha )C_{i+1}^{j+1}={}$

$\qquad \beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}.$

Чтобы решить эту линейную систему уравнений, мы должны теперь увидеть, что сначала необходимо задать граничные условия в начале каналов:

$C_{0}^{j}$: начальное состояние канала на текущем временном шаге,

$C_{0}^{j+1}$: начальное состояние канала на следующем временном шаге,

$C_{N0}^{j}$: исходное состояние для предыдущего канала по отношению к каналу, анализируемому на текущем временном шаге,

$C_{M0}^{j}$: начальное состояние для канала, следующего за каналом, анализируемым на текущем временном шаге.

Для последней ячейки каналов ( $z$), наиболее удобным условием становится адиабатическое, поэтому

$\left.{\frac {\partial C}{\partial x}}\right|_{x=z}={\frac {C_{i+1}-C_{i-1}}{2\,\Delta x}}=0.$

Это условие выполняется тогда и только тогда, когда (независимо от нулевого значения)

$C_{i+1}^{j+1}=C_{i-1}^{j+1}.$

Решим эту задачу (в матричной форме) для случая 3 каналов и 5 узлов (включая начальное граничное условие). Мы выражаем это как задачу линейной системы:

$\mathbf {AA} \,\mathbf {C^{j+1}} =\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ,$

где

$\mathbf {C^{j+1}} ={\begin{bmatrix}C_{11}^{j+1}\\C_{12}^{j+1}\\C_{13}^{j+1}\\C_{14}^{j+1}\\C_{21}^{j+1}\\C_{22}^{j+1}\\C_{23}^{j+1}\\C_{24}^{j+1}\\C_{31}^{j+1}\\C_{32}^{j+1}\\C_{33}^{j+1}\\C_{34}^{j+1}\end{bmatrix}},\quad \mathbf {C^{j}} ={\begin{bmatrix}C_{11}^{j}\\C_{12}^{j}\\C_{13}^{j}\\C_{14}^{j}\\C_{21}^{j}\\C_{22}^{j}\\C_{23}^{j}\\C_{24}^{j}\\C_{31}^{j}\\C_{32}^{j}\\C_{33}^{j}\\C_{34}^{j}\end{bmatrix}}.$

Теперь мы должны понимать, что AA и BB должны быть массивами, состоящими из четырех разных подмассивов (помните, что в этом примере рассматриваются только три канала, но он охватывает основную часть, рассмотренную выше):

$\mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}},\quad \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}},$

где упомянутые выше элементы соответствуют следующим массивам и дополнительным 4×4, полным нулей. Обратите внимание, что размеры АА и ББ — 12×12:

$\mathbf {AA1} ={\begin{bmatrix}(1+2\lambda +\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +\beta )\end{bmatrix}},$

$\mathbf {AA2} ={\begin{bmatrix}(1+2\lambda +2\beta )&-(\lambda -\alpha )&0&0\\-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )&0\\0&-(\lambda +\alpha )&(1+2\lambda +2\beta )&-(\lambda -\alpha )\\0&0&-2\lambda &(1+2\lambda +2\beta )\end{bmatrix}},$

$\mathbf {AA3} ={\begin{bmatrix}-\beta &0&0&0\\0&-\beta &0&0\\0&0&-\beta &0\\0&0&0&-\beta \end{bmatrix}},$

$\mathbf {BB1} ={\begin{bmatrix}(1-2\lambda -\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -\beta )\end{bmatrix}},$

$\mathbf {BB2} ={\begin{bmatrix}(1-2\lambda -2\beta )&(\lambda -\alpha )&0&0\\(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )&0\\0&(\lambda +\alpha )&(1-2\lambda -2\beta )&(\lambda -\alpha )\\0&0&2\lambda &(1-2\lambda -2\beta )\end{bmatrix}}.$

Вектор d здесь используется для хранения граничных условий. В этом примере это вектор 12×1:

$\mathbf {d} ={\begin{bmatrix}(\lambda +\alpha )(C_{10}^{j+1}+C_{10}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{20}^{j+1}+C_{20}^{j})\\0\\0\\0\\(\lambda +\alpha )(C_{30}^{j+1}+C_{30}^{j})\\0\\0\\0\end{bmatrix}}.$

Чтобы найти концентрацию в любой момент времени, необходимо выполнить следующее уравнение:

$\mathbf {C^{j+1}} =\mathbf {AA} ^{-1}(\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ).$

Пример: 2D-диффузия

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

${\frac {\partial u}{\partial t}}=a\,\nabla ^{2}u,$

${\frac {\partial u}{\partial t}}=a\left({\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}\right)$

можно решить с помощью дискретизации Кранка – Николсона

${\begin{aligned}u_{i,j}^{n+1}={}&u_{i,j}^{n}+{\frac {1}{2}}{\frac {a\Delta t}{(\Delta x)^{2}}}{\big [}(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}-4u_{i,j}^{n+1})\\&+(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}-4u_{i,j}^{n}){\big ]},\end{aligned}}$

предполагая, что используется квадратная сетка, так что $\Delta x=\Delta y$. Это уравнение можно несколько упростить, переставив члены и используя число CFL.

$\mu ={\frac {a\,\Delta t}{(\Delta x)^{2}}}.$

Для численной схемы Кранка – Николсона низкое число КЛЛ не требуется для стабильности, однако оно необходимо для численной точности. Теперь мы можем записать схему как

$(1+2\mu )u_{i,j}^{n+1}-{\frac {\mu }{2}}\left(u_{i+1,j}^{n+1}+u_{i-1,j}^{n+1}+u_{i,j+1}^{n+1}+u_{i,j-1}^{n+1}\right)$

$\qquad =(1-2\mu )u_{i,j}^{n}+{\frac {\mu }{2}}\left(u_{i+1,j}^{n}+u_{i-1,j}^{n}+u_{i,j+1}^{n}+u_{i,j-1}^{n}\right).$

Решение такой линейной системы является дорогостоящим. Следовательно, для решения численного УЧП может быть реализован неявный метод чередования направлений , при котором одно измерение обрабатывается неявно, а другое измерение явно для половины назначенного временного шага и, наоборот, для оставшейся половины временного шага. Преимущество этой стратегии заключается в том, что неявного решателя требуется только алгоритм трехдиагональной матрицы для решения . Разница между истинным решением Кранка – Николсона и приближенным решением ADI имеет порядок точности $O(\Delta t^{4})$ и, следовательно, им можно пренебречь с достаточно малым шагом по времени.

Кранка – Николсона для нелинейных задач

Поскольку метод Кранка-Николсона является неявным , точное решение обычно невозможно. Вместо этого для сходимости к решению следует использовать итерационный метод. Один из вариантов — использовать метод Ньютона для сходимости предсказания, но для этого требуется вычисление якобиана . Для многомерных систем, подобных тем, которые используются в вычислительной гидродинамике или численной теории относительности , вычислить этот якобиан может оказаться невозможным.

Альтернативой без якобиана является итерация с фиксированной точкой . Если $f$ - скорость системы, то предсказание Кранка-Николсона будет фиксированной точкой на карте. $\Phi (x)=x_{0}+{\frac {h}{2}}\left[f(x_{0})+f(x)\right].$ Если итерация карты $x^{(i+1)}=\Phi (x^{(i)})$ не сходится, параметризованное отображение $\Theta (x,\alpha )=\alpha x+(1-\alpha )\Phi (x)$, с $\alpha \in (0,1)$, может вести себя лучше. В развернутом виде формула обновления имеет вид

$x^{i+1}=\alpha x^{i}+(1-\alpha )\left[x_{0}+{\frac {h}{2}}\left(f(x_{0})+f(x^{i})\right)\right],$

где $x^{i}$ это текущее предположение и $x_{i-1}$ это предыдущий временной шаг.

Даже для многомерных систем итерация этого отображения может сходиться на удивление быстро.

Duration: 11 seconds.0:11

Численное решение уравнений Навье–Стокса в форме завихренности. В этом случае $\alpha =7/8$ было необходимо для сходимости итерации Кранка – Николсона с фиксированной точкой.

Применение в финансовой математике

Дополнительная информация: Методы конечных разностей для ценообразования опционов.

ряд других явлений ) можно смоделировать Поскольку с помощью уравнения теплопроводности (часто называемого уравнением диффузии в финансовой математике , к этим областям также был применен метод Кранка-Николсона. В частности, Блэка-Шоулза модели ценообразования опционов дифференциальное уравнение можно преобразовать в уравнение теплопроводности, и, таким образом, численные решения для ценообразования опционов можно получить с помощью метода Кранка-Николсона.

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

См. также

• Финансовая математика
• Правило трапеции

Внешние ссылки

• Численные методы PDE для ученых и инженеров , лекции и коды для числовых PDE в открытом доступе
• Пример применения и реализации метода Кранка – Николсона для уравнения адвекции.