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

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

В численном анализе метод Кранка –Николсона представляет собой метод конечных разностей, используемый для численного решения уравнения теплопроводности и подобных уравнений в частных производных . ^[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}$

Пример: одномерная диффузия

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

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

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

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

Здесь мы моделируем концентрацию растворенного загрязняющего вещества в воде. Эта задача состоит из трех частей: известного уравнения диффузии ( выбранного как константа), адвективного компонента (который означает, что система развивается в пространстве из-за поля скорости), который мы выбираем как константу , и бокового взаимодействия между продольными каналами ( ): ${\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, заполненные нулями. Обратите внимание, что размеры AA и BB составляют 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}}}.}$

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

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

Решение такой линейной системы является дорогостоящим. Следовательно, для решения численного PDE может быть реализован неявный метод с чередующимся направлением , при котором одно измерение обрабатывается неявно, а другое измерение явно для половины назначенного временного шага и наоборот для оставшейся половины временного шага. Преимущество этой стратегии заключается в том, что неявному решателю требуется только алгоритм трехдиагональной матрицы для решения. Разница между истинным решением Крэнка–Николсона и приближенным решением 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. Crank, J.; Nicolson, P. (1947). "Практический метод численной оценки решений уравнений с частными производными типа теплопроводности". Proc. Camb. Phil. Soc . 43 (1): 50–67. Bibcode :1947PCPS...43...50C. doi :10.1017/S0305004100023197. S2CID  16676040.
3. Томас, Дж. В. (1995). Численные уравнения в частных производных: методы конечных разностей . Тексты по прикладной математике. Том 22. Берлин, Нью-Йорк: Springer-Verlag . ISBN 978-0-387-97999-1.. Пример 3.3.2 показывает, что уравнение Крэнка–Николсона безусловно устойчиво при применении к . ${\displaystyle u_{t}=au_{xx}}$
4. "Многомерные параболические задачи" (PDF) . Computer Science Department . RPI . Получено 29 мая 2016 г. .
5. Уилмотт, П.; Хоуисон, С.; Дьюинн, Дж. (1995). Математика финансовых производных: введение для студентов . Cambridge Univ. Press. ISBN 0-521-49789-2. Математика финансовых производных инструментов Уилмотт.

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

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