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

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

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

${\ 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 {\ text {( обратный Эйлер)}},}$

${\ 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 {\ text {(Crank – Nicolson)}}.}$

Обратите внимание, что это неявный метод : чтобы получить «следующее» значение 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 г = {\ гидроразрыва {а \ 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 и_ {я} ^ {п + 1}}$может быть эффективно решена с помощью алгоритма трехдиагональной матрицы в пользу гораздо более дорогостоящего обращения матрицы .

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

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

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

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

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

Метод Кранка – Николсона (где i представляет положение, а j время) преобразует каждый компонент 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} = {}}$

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

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

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

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

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

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

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

Для последней ячейки каналов ( 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 \ 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 \ лямбда -2 \ бета) & (\ лямбда - \ альфа) & 0 \\ 0 & (\ лямбда + \ альфа) & (1-2 \ лямбда -2 \ бета) & (\ лямбда - \ альфа) \\ 0 & 0 & 2 \ лямбда & (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}). }$

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

${\ displaystyle {\ frac {\ partial u} {\ partial t}} = а \ набла ^ {2} и,}$

${\ displaystyle {\ frac {\ partial u} {\ partial t}} = a \ left ({\ frac {\ partial ^ {2} u} {\ partial x ^ {2}}} + {\ frac {\ частичный ^ {2} u} {\ partial y ^ {2}}} \ right)}$

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

${\ displaystyle {\ begin {align} u_ {i, j} ^ {n + 1} & = u_ {i, j} ^ {n} + {\ frac {1} {2}} {\ frac {a \ Дельта 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 {align}}}$

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

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

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

Одна альтернатива - определить карту, для которой предсказание Крэнка-Николсона является фиксированной точкой. Умные итерации этой карты могут привести к правильному прогнозу. Обратите внимание, что для нелинейных задач не гарантируется, что прогноз существует или даже уникален. Для достаточно малых размеров шага эти итерационные методы должны быть успешными. Позволять${\ Displaystyle {\ точка {х}} = е (х),}$ где ${\ Displaystyle х \ в \ mathbb {R} ^ {п}}$ а также ${\ displaystyle f: \ mathbb {R} ^ {n} \ rightarrow \ mathbb {R} ^ {n}}$. Прогноз Крэнка-Николсона будет фиксированной точкой карты.${\ displaystyle \ Phi (x) = x_ {0} + {\ frac {h} {2}} \ left [f (x_ {0}) + f (x) \ right].}$ Если итерация карты ${\ Displaystyle х ^ {(я + 1)} = \ Phi (х ^ {(я)})}$ не сходится, параметризованная карта ${\ Displaystyle \ Theta (х, \ альфа) = \ альфа х + (1- \ альфа) \ Phi (x)}$, с участием ${\ Displaystyle \ альфа \ в (0,1)}$, может вести себя лучше.