Метод Кранка – Николсона - Crank–Nicolson method

В числовой анализ, то Метод Кранка – Николсона это метод конечных разностей используется для численного решения уравнение теплопроводности и подобные уравнения в частных производных.^[1] Это второго порядка метод во времени. это скрытый во времени и может быть записан как неявный метод Рунге – Кутты, и это численно стабильный. Метод был разработан Джон Крэнк и Филлис Николсон в середине 20 века.^[2]

За уравнения диффузии (и многие другие уравнения), можно показать, что метод Кранка – Николсона безусловно стабильный.^[3] Однако приближенные решения могут все же содержать (затухающие) паразитные колебания, если отношение шага по времени Δт раз температуропроводность к квадрату пространственного шага, ΔИкс², большой (обычно больше 1/2 на Анализ устойчивости фон Неймана ). По этой причине, когда необходимы большие временные шаги или высокое пространственное разрешение, менее точный обратный метод Эйлера часто используется, который является стабильным и невосприимчивым к колебаниям.^{[нужна цитата ]}

Метод

Шаблон Кранка – Николсона для одномерной задачи.

Метод Кранка – Николсона основан на трапеция, дающий сходимость второго порядка по времени. Для линейных уравнений правило трапеций эквивалентно неявный метод средней точки^{[нужна цитата ]} - простейший пример Гаусс-Лежандр неявный метод Рунге-Кутта, который также имеет свойство быть геометрический интегратор. Например, в одном измерении предположим, что уравнение в частных производных является

${\displaystyle {\frac {\partial u}{\partial t}}=F\left(u,\,x,\,t,\,{\frac {\partial u}{\partial x}},\,{\frac {\partial ^{2}u}{\partial x^{2}}}\right)}$

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

${\displaystyle {\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)\qquad {\mbox{(forward Euler)}}}$

${\displaystyle {\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)\qquad {\mbox{(backward Euler)}}}$

${\displaystyle {\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]\qquad {\mbox{(Crank-Nicolson)}}.}$

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

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

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

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

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

${\displaystyle {\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)}$

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

${\displaystyle -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}\,}$

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

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

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

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

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

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

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

${\displaystyle {\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 соответствуют предыдущий и следующий канал.

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

${\displaystyle {\frac {\partial C}{\partial t}}\Rightarrow {\frac {C_{i}^{j+1}-C_{i}^{j}}{\Delta t}}}$

(2)

${\displaystyle {\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)

${\displaystyle {\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)

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

(5)

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

(6)

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

(7)

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

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

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

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

и заменить (2), (3), (4), (5), (6), (7), α, β и λ в (1). Затем мы помещаем Новое время условия слева (j + 1) и настоящее время условия справа (j) получить:

${\displaystyle -\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}=\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), поэтому выражение упрощается до:

${\displaystyle -(\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}=+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}+\beta C_{Mi}^{j}.}$

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

${\displaystyle -\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}=\beta C_{Ni}^{j}+(\lambda +\alpha )C_{i-1}^{j}+(1-2\lambda -\beta )C_{i}^{j}+(\lambda -\alpha )C_{i+1}^{j}.}$

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

${\displaystyle C_{0}^{j}}$: начальное условие для канала на текущем временном шаге
${\displaystyle C_{0}^{j+1}}$: начальное условие для канала на следующем временном шаге
${\displaystyle C_{N0}^{j}}$: начальное условие для предыдущего канала к проанализированному в текущий момент времени.
${\displaystyle C_{M0}^{j}}$: начальное условие для следующего канала после анализируемого на текущем временном шаге.

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

${\displaystyle {\frac {\partial C}{\partial x}}_{x=z}={\frac {(C_{i+1}-C_{i-1})}{2\Delta x}}=0.}$

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

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

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

${\displaystyle {\begin{bmatrix}AA\end{bmatrix}}{\begin{bmatrix}C^{j+1}\end{bmatrix}}=[BB][C^{j}]+[d]}$

куда

${\displaystyle \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}}}$ и ${\displaystyle \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 должны быть массивами, состоящими из четырех различных подмассивов (помните, что в этом примере рассматриваются только три канала, но он охватывает основную часть, описанную выше).

${\displaystyle \mathbf {AA} ={\begin{bmatrix}AA1&AA3&0\\AA3&AA2&AA3\\0&AA3&AA1\end{bmatrix}}}$ и

${\displaystyle \mathbf {BB} ={\begin{bmatrix}BB1&-AA3&0\\-AA3&BB2&-AA3\\0&-AA3&BB1\end{bmatrix}}}$  

где упомянутые выше элементы соответствуют следующим массивам и дополнительному 4x4, заполненному нулями. Обратите внимание, что размеры AA и BB составляют 12x12:

${\displaystyle \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}}}$   , 

${\displaystyle \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}}}$   , 

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

${\displaystyle \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}}}$   &  

${\displaystyle \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 вектор здесь используется для выполнения граничных условий. В этом примере это вектор 12x1:

${\displaystyle \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}}.}$

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

${\displaystyle {\begin{bmatrix}C^{j+1}\end{bmatrix}}={\begin{bmatrix}AA^{-1}\end{bmatrix}}([BB][C^{j}]+[d]).}$

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

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

${\displaystyle {\partial u \over \partial t}=a\nabla ^{2}u}$

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

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

${\displaystyle {\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})\\&\qquad {}+(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}}}$

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

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

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

${\displaystyle {\begin{aligned}&(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)\\&\quad =(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).\end{aligned}}}$

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

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

Дальнейшая информация: Конечно-разностные методы ценообразования опционов

Потому что ряд других явлений может быть смоделированный с уравнение теплопроводности (часто называемое уравнением диффузии в финансовая математика ), метод Кранка – Николсона также был применен к этим областям.^[5] В частности, Блэк – Скоулз модель ценообразования опционов дифференциальное уравнение можно преобразовать в уравнение теплопроводности, и, таким образом, численные решения за опционная цена может быть получен методом Кранка – Николсона.

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

Смотрите также

• Финансовая математика
• Трапецеидальная линейка

Рекомендации

1. Tuncer Cebeci (2002). Конвективная теплопередача. Springer. ISBN  0-9668461-4-1.
2. Crank, J .; Николсон, П. (1947). «Практический метод численной оценки решений дифференциальных уравнений в частных производных типа теплопроводности». Proc. Camb. Фил. Soc. 43 (1): 50–67. Дои:10.1017 / S0305004100023197..
3. Томас, Дж. У. (1995). Численные уравнения с частными производными: конечно-разностные методы. Тексты по прикладной математике. 22. Берлин, Нью-Йорк: Springer-Verlag. ISBN  978-0-387-97999-1.. Пример 3.3.2 показывает, что Кранк – Николсон безусловно устойчив в применении к ${\displaystyle u_{t}=au_{xx}}$.
4. «Многомерные параболические задачи» (PDF). Кафедра компьютерных наук. RPI. Получено 29 мая 2016.
5. Wilmott, P .; Howison, S .; Дьюинн, Дж. (1995). Математика финансовых деривативов: введение для студентов. Cambridge Univ. Нажмите. ISBN  0-521-49789-2. Математика производных финансовых инструментов Уилмотт.

внешняя ссылка

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