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

Метод конечных разностей для численного решения параболических дифференциальных уравнений

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

Для уравнений диффузии (и многих других уравнений) можно показать, что метод Кранка – Николсона безоговорочно устойчив . ^[3] Однако приближенные решения могут по-прежнему содержать (затухающие) паразитные колебания, если отношение шага по времени, умноженного на температуропроводность, к квадрату пространственного шага, велико (обычно больше, чем 1/2 на анализ устойчивости фон Неймана ). . По этой причине, когда необходимы большие временные шаги или высокое пространственное разрешение, часто используется менее точный обратный метод Эйлера , который одновременно стабилен и невосприимчив к колебаниям. ^{[ нужна цитата ]} ${\displaystyle \Delta t}$ ${\displaystyle \Delta x^{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}$ ${\displaystyle n+1}$

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

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

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

${\displaystyle {\frac {\partial u}{\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}}$ ${\displaystyle U_{x}}$ ${\displaystyle k}$

где – концентрация примеси, а индексы и соответствуют предыдущему и следующему каналу. ${\displaystyle C}$ ${\displaystyle N}$ ${\displaystyle M}$

Метод Кранка-Николсона (где представляет положение и время) преобразует каждый компонент УЧП в следующее: ${\displaystyle i}$ ${\displaystyle j}$

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

${\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 ). Затем мы помещаем новые временные условия слева ( ) и текущие временные условия справа ( ), чтобы получить ${\displaystyle \alpha }$ ${\displaystyle \beta }$ ${\displaystyle \lambda }$ ${\displaystyle j+1}$ ${\displaystyle 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}={}}$

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

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

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

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

${\displaystyle \qquad \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}}$: начальное состояние для канала, следующего за каналом, анализируемым на текущем временном шаге.

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

${\displaystyle \left.{\frac {\partial C}{\partial x}}\right|_{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 \mathbf {AA} \,\mathbf {C^{j+1}} =\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {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}},\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 должны быть массивами, состоящими из четырех разных подмассивов (помните, что в этом примере рассматриваются только три канала, но он охватывает основную часть, рассмотренную выше):

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

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

${\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 \mathbf {C^{j+1}} =\mathbf {AA} ^{-1}(\mathbf {BB} \,\mathbf {C^{j}} +\mathbf {d} ).}$

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

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

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

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

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

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

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

${\displaystyle \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 имеет порядок точности и, следовательно, может быть проигнорирована с достаточно малым шагом по времени. ^[4] ${\displaystyle O(\Delta t^{4})}$

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

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

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

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

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

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

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

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

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

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

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

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

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

1. Тунцер Чебечи (2002). Конвективный теплообмен. Спрингер. ISBN 0-9668461-4-1.
2. Кранк, Дж.; Николсон, П. (1947). «Практический метод численного определения решений уравнений в частных производных типа теплопроводности». Учеб. Кэмб. Фил. Соц . 43 (1): 50–67. Бибкод : 1947PCPS...43...50C. дои : 10.1017/S0305004100023197. S2CID  16676040.
3. Томас, JW (1995). Численные уравнения в частных производных: конечно-разностные методы . Тексты по прикладной математике. Том. 22. Берлин, Нью-Йорк: Springer-Verlag . ISBN 978-0-387-97999-1.. Пример 3.3.2 показывает, что Крэнк–Николсон безусловно устойчив применительно к . ${\displaystyle u_{t}=au_{xx}}$
4. «Многомерные параболические задачи» (PDF) . Кафедра компьютерных наук . РПИ . Проверено 29 мая 2016 г.
5. Уилмотт, П.; Хауисон, С.; Дьюинн, Дж. (1995). Математика производных финансовых инструментов: введение для студентов . Кембриджский университет. Нажимать. ISBN 0-521-49789-2. Математика финансовых производных Уилмотт.

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

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