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

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

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

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

Содержание
  • 1 Метод
  • 2 Пример: 1D диффузия
  • 3 Пример: 1D диффузия с адвекцией для устойчивого потока, с несколькими соединениями каналов
  • 4 Пример: 2D диффузия
  • 5 Применение в финансовой математике
  • 6 См. Также
  • 7 Ссылки
  • 8 Внешние ссылки
Метод
Трафарет Крэнка – Николсона для одномерной задачи.

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

∂ u ∂ t = F (u, x, t, ∂ u ∂ x, ∂ 2 u ∂ x 2) {\ displaystyle {\ frac {\ partial u} {\ partial t}} = F \ left (u, \, x, \, t, \, {\ frac {\ partial u} {\ partial x}}, \, {\ гидроразрыва {\ partial ^ {2} u} {\ partial x ^ {2}}} \ right)}{\ frac {\ partial u} {\ partial t}} = F \ left (u, \, x, \, t, \, {\ frac {\ partial u} {\ partial x}}, \, {\ frac {\ partial ^ {2} u} {\ partial x ^ { 2}}} \ right)

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

uin + 1 - uin Δ t = F in ( u, x, t, ∂ u ∂ x, ∂ 2 u ∂ x 2) (прямой Эйлер) {\ displ aystyle {\ 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 { \ t_dv {(прямой Эйлер)}}}{\ 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) \ q quad {\ t_dv {(вперед Эйлер)}}
uin + 1 - uin Δ t = F in + 1 (u, x, t, ∂ u ∂ x, ∂ 2 u ∂ x 2) (обратный Эйлер) {\ 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 {\ t_dv {(обратный Эйлер)}}}{\ 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} {\ частичный x ^ {2}}} \ right) \ qquad {\ t_dv {(обратный Эйлер)}}
uin + 1 - uin Δ t = 1 2 [F in + 1 (u, x, t, ∂ u ∂ x, ∂ 2 u ∂ x 2) + F in (u, x, t, ∂ u ∂ x, ∂ 2 u ∂ x 2)] (Кранк-Николсон). {\ 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 {\ t_dv {(Crank-Nicolson)}}.}{\ 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 {\ t_dv {(Crank-Nicolson)}}.

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

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

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

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

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

uin + 1 - uin Δ t знак равно a 2 (Δ Икс) 2 ((ui + 1 n + 1-2 uin + 1 + ui - 1 n + 1) + (ui + 1 n - 2 uin + ui - 1 n)) {\ 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)}{\ 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 = a Δ t 2 (Δ x) 2 {\ displaystyle r = {\ frac {a \ Delta t} {2 (\ Delta x) ^ {2}}}}{\ displaystyle r = {\ frac {a \ Delta t} {2 (\ Delta x) ^ {2}}}} :

- rui + 1 n + 1 + (1 + 2 r) uin + 1 - rui - 1 n + 1 знак равно rui + 1 n + (1-2 r) uin + rui - 1 n {\ 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} \,}-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}} \,

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

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

∂ u ∂ t = a (u) ∂ 2 U ∂ Икс 2 {\ Displaystyle {\ frac {\ partial u} {\ partial t}} = a (u) {\ frac {\ partial ^ {2} u} {\ partial x ^ {2}} }}{\ frac {\ partial u} {\ partial t}} = a (u) {\ frac {\ partial ^ {2} u} {\ partial x ^ {2}}}

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

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

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

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

∂ C ∂ T знак равно D x ∂ 2 C ∂ Икс 2 - U x ∂ C ∂ x - k (C - CN) - k (C - CM) {\ 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})}{\ 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) преобразует каждый компонент PDE в следующее:

∂ C ∂ t ⇒ C ij + 1 - C ij Δ t {\ displaystyle {\ frac {\ partial C} {\ partial t}} \ Rightarrow {\ frac {C_ {i} ^ {j + 1} -C_ {i} ^ {j}} {\ Delta t}}}{\ frac {\ partial C} {\ partial t}} \ Rightarrow {\ frac {C _ {{i}} ^ {{j + 1}} - C _ {{i}} ^ {{j}} } {\ Delta t}}

(2)

∂ 2 C ∂ x 2 ⇒ 1 2 (Δ x) 2 ((C i + 1 j + 1-2 C ij + 1 + C i - 1 j + 1) + (C i + 1 j - 2 C ij + C i - 1 j)) {\ displaystyle {\ frac {\ partial ^ {2} C} {\ partial x ^ {2}}} \ Rightarrow {\ frac {1} {2 (\ Дельта 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)}{\ 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)

∂ C ∂ x ⇒ 1 2 ((C я + 1 j + 1 - C я - 1 j + 1) 2 (Δ x) + (C i + 1 j - C я - 1 j) 2 (Δ x)) {\ 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)}{\ 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 ⇒ 1 2 (C ij + 1 + C ij) {\ displaystyle C \ Rightarrow {\ frac {1} {2}} ( C_ {i} ^ {j + 1} + C_ {i} ^ { j})}C \ Rightarrow {\ frac {1} {2}} (C _ {{i}} ^ { {j + 1}} + C _ {{i}} ^ {{j}})

(5)

CN ⇒ 1 2 (CN ij + 1 + CN ij) {\ displaystyle C_ {N} \ Rightarrow {\ frac {1} {2}} (C_ {Ni} ^ {j + 1} + C_ {Ni} ^ {j})}C_ {N} \ Ri ghtarrow {\ frac {1} {2}} (C _ {{Ni}} ^ {{j + 1}} + C _ {{Ni}} ^ {{j}})

(6)

CM ⇒ 1 2 (CM ij + 1 + CM ij). {\ displaystyle C_ {M} \ Rightarrow {\ frac {1} {2}} (C_ {Mi} ^ {j + 1} + C_ {Mi} ^ {j}).}C_ {M} \ Rightarrow {\ frac {1} {2}} (C _ {{Mi}} ^ {{j + 1}} + C _ {{Mi} } ^ {{j}}).

(7)

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

λ = D x Δ t 2 Δ x 2 {\ displaystyle \ lambda = {\ frac {D_ {x} \ Delta t} {2 \ Delta x ^ { 2}}}}\ lambda = {\ frac {D_ {x} \ Delta t} {2 \ Дельта x ^ {2}}}
α = U x Δ t 4 Δ x {\ displaystyle \ alpha = {\ frac {U_ {x} \ Delta t} {4 \ Delta x}}}\ alpha = {\ frac {U_ {x} \ Delta t} {4 \ Delta x}}
β = k Δ t 2 {\ displaystyle \ beta = {\ frac {k \ Delta t} {2}}}\ beta = {\ frac {k \ Delta t} {2}}

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

- β CN ij + 1 - (λ + α) C i - 1 j + 1 + (1 + 2 λ + 2 β) C ij + 1 - (λ - α) C i + 1 j + 1 - β CM ij + 1 = β CN ij + (λ + α) C i - 1 j + (1-2 λ - 2 β) C ij + (λ - α) C i + 1 j + β CM ij. {\ 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}.}- \ 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), поэтому выражение упрощается до:

- (λ + α) C i - 1 j + 1 + (1 + 2 λ + β) C ij + 1 - (λ - α) C i + 1 j + 1 - β CM ij + 1 = + (λ + α) C i - 1 j + (1-2 λ - β) C ij + (λ - α) C i + 1 j + β CM ij. {\ displaystyle - (\ lambda + \ alpha) C_ {i-1} ^ {j + 1} + (1 + 2 \ lambda + \ beta) C_ {i} ^ {j + 1} - (\ lambda - \ альфа) 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}.}- (\ 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), поэтому выражение упрощается до:

- β CN ij + 1 - (λ + α) C i - 1 j + 1 + (1 + 2 λ + β) C ij + 1 - (λ - α) C i + 1 j + 1 = β CN ij + (λ + α) C i - 1 j + (1 - 2 λ - β) C ij + (λ - α) C i + 1 j. {\ Displaystyle - \ бета 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}.}- \ 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}}.

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

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

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

∂ C ∂ xx = z = (C i + 1 - C i - 1) 2 ∆ x Знак равно 0. {\ displaystyle {\ frac {\ partial C} {\ partial x}} _ {x = z} = {\ frac {(C_ {i + 1} -C_ {i-1})} {2 \ Дельта x}} = 0.}{\ frac {\ partial C} {\ partial x}} _ {{x = z}} = {\ frac {(C _ {{i + 1}} -C _ {{i-1}})} {2 \ Delta x}} = 0.

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

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

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

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

где

C j + 1 = [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 С 34 j + 1] {\ 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}}}{\ 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}} и C j = [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]. {\ 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}}.}{\ 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 должны быть массивами, состоящими из четырех разных подмассивов (помните, что в этом примере рассматриваются только три канала, но он охватывает основную часть, описанную выше).

AA = [AA 1 AA 3 0 AA 3 AA 2 AA 3 0 AA 3 AA 1] {\ displaystyle \ mathbf {AA} = {\ begin {bmatrix} AA1 AA3 0 \\ AA3 AA2 AA3 \\ 0 AA3 AA1 \ end {bmatrix} }}{\ mathbf {AA}} = {\ begin {bmatrix} AA1 AA3 0 \\ AA3 AA2 AA3 \\ 0 AA3 AA1 \ end {bmatrix}} и
BB = [BB 1 - AA 3 0 - AA 3 BB 2 - AA 3 0 - AA 3 BB 1] {\ displaystyle \ mathbf {BB} = {\ begin {bmatrix } BB1 -AA3 0 \\ - AA3 BB2 -AA3 \\ 0 -AA3 BB1 \ end {bmatrix}}}{\ mathbf {BB}} = {\ begin {bmatrix} BB1 -AA3 0 \\ - AA3 BB2 -AA3 \\ 0 -AA3 BB1 \ end {bmatrix}}

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

AA 1 = [(1 + 2 λ + β) - (λ - α) 0 0 - (λ + α) (1 + 2 λ + β) - (λ - α) 0 0 - (λ + α) (1 + 2 λ + β) - (λ - α) 0 0-2 λ (1 + 2 λ + β)] {\ 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}}}{\ 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}} ,
AA 2 = [(1 + 2 λ + 2 β) - (λ - α) 0 0 - (λ + α) (1 + 2 λ + 2 β) - (λ - α) 0 0 - (λ + α) (1 + 2 λ + 2 β) - (λ - α) 0 0 - 2 λ (1 + 2 λ + 2 β)] {\ 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}}}{\ 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}} ,
AA 3 = [- β 0 0 0 0 - β 0 0 0 0 - β 0 0 0 0 - β] {\ Displaystyle \ math bf {AA3} = {\ begin {bmatrix} - \ beta 0 0 0 \\ 0 - \ beta 0 0 \\ 0 0 - \ beta 0 \\ 0 0 0 - \ 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}} ,
BB 1 = [( 1 - 2 λ - β) (λ - α) 0 0 (λ + α) (1-2 λ - β) (λ - α) 0 0 (λ + α) (1-2 λ - β) (λ - α) 0 0 2 λ (1-2 λ - β)] {\ 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}}}{\ mathbf {BB1}} = {\ begin {bmatrix} (1-2 \ lambda - \ beta) (\ lambda - \ alpha) 0 0 \\ (\ лямбда + \ альфа) (1-2 \ лямбда - \ бета) (\ лямбда - \ альфа) 0 \\ 0 (\ лямбда + \ альфа) (1-2 \ лямбда - \ бета) (\ лямбда - \ alpha) \\ 0 0 2 \ lambda (1-2 \ lambda - \ beta) \ end {bmatrix}}
BB 2 = [(1-2 λ - 2 β) (λ - α) 0 0 (λ + α) (1-2 λ - 2 β) (λ - α) 0 0 (λ + α) (1-2 λ - 2 β) (λ - α) 0 0 2 λ (1 - 2 λ - 2 β)]. {\ displaystyle \ mathbf {BB2} = {\ begin {bmatrix} (1-2 \ lambda -2 \ beta) (\ lambda - \ alpha) 0 0 \\ (\ lambda + \ alpha) (1-2 \ лямбда -2 \ бета) (\ лямбда - \ альфа) 0 \\ 0 (\ лямбда + \ альфа) (1-2 \ лямбда -2 \ бета) (\ лямбда - \ альфа) \\ 0 0 2 \ лямбда (1-2 \ lambda -2 \ beta) \ end {bmatrix}}.}{\ mathbf {BB2}} = {\ begin {bmatrix} (1-2 \ lambda -2 \ бета) (\ lambda - \ alpha) 0 0 \\ (\ lambda + \ alpha) (1-2 \ l ambda -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:

d = [(λ + α) (C 10 j + 1 + C 10 j) 0 0 0 (λ + α) (C 20 j + 1 + C 20 j) 0 0 0 (λ + α) (C 30 j + 1 + C 30 j) 0 0 0]. {\ 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}}.}{\ 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 }}.

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

[C j + 1] = [AA - 1] ([BB] [C j] + [d]). {\ displaystyle {\ begin {bmatrix} C ^ {j + 1} \ end {bmatrix}} = {\ begin {bmatrix} AA ^ {- 1} \ end {bmatrix}} ([BB] [C ^ {j }] + [d]).}{\ begin {bmatrix} C ^ {{j +1}} \ end {bmatrix}} = {\ begin {bmatrix} AA ^ {{- 1}} \ end {bmatrix}} ([BB] [C ^ {{j}}] + [d]).
Пример: 2D-диффузия

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

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

.

∂ u ∂ t = a (∂ 2 u ∂ Икс 2 + ∂ 2 u ∂ Y 2) {\ Displaystyle {\ frac {\ partial u} {\ partial t}} = a \ left ({\ frac {\ partial ^ {2} u} { \ partial x ^ {2}}} + {\ frac {\ partial ^ {2} u} {\ partial y ^ {2}}} \ right)}{\ frac {\ partial u} {\ partial t}} = a \ left ({\ frac {\ partial ^ { 2} u} {\ partial x ^ {2}}} + {\ frac {\ partial ^ {2} u} {\ partial y ^ {2}}} \ right)

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

ui, jn + 1 = ui, jn + 1 2 a Δ t (Δ x) 2 [(ui + 1, jn + 1 + ui - 1, jn + 1 + ui, j + 1 n + 1 + ui, j - 1 n + 1 - 4 ui, jn + 1) + (ui + 1, jn + ui - 1, jn + ui, j + 1 n + ui, j - 1 n - 4 ui, jn)] { \ Displaystyle {\ begin {align} 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 {align}}}{\ begin {align} 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 {align}}

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

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

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

(1 + 2 μ) ui, jn + 1 - μ 2 (ui + 1, jn + 1 + ui - 1, jn + 1 + ui, j + 1 n + 1 + ui, j - 1 n + 1) = (1-2 μ) ui, jn + μ 2 (ui + 1, jn + ui - 1, jn + ui, j + 1 n + ui, j - 1 n). {\ displaystyle {\ begin {align} (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 {align}}}{\ begin {align} (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 {align}}

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

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

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

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

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

.

Внешние ссылки
Последняя правка сделана 2021-05-16 07:53:48
Содержание доступно по лицензии CC BY-SA 3.0 (если не указано иное).
Обратная связь: support@alphapedia.ru
Соглашение
О проекте