Конечно-разностный метод решения для уравнений параболического типа (стр. 3 из 4)

, (20)

т. е. условие (20) является необходимым условием устойчивости.

1.5 Схема Кранка-Николсона

параболическое дифференциальное уравнение конечная разность

Явная конечно разностная схема, записанная в форме

(21)

обладает тем достоинством, что решение на верхнем временном слое t_k+l получается сразу (без решения СЛАУ) по значениям сеточной функции на нижнем временном слое t^k, где решение известно (при k = 0 значения сеточной функции формируются из начального условия). Но эта же схема обладает существенным недостатком, поскольку она является условно устойчивой. С другой стороны, неявная конечно-разностная схема, записанная форме

(22)

приводит к необходимости решать СЛАУ, но зато эта схема абсолютно устойчива.

Проанализируем схемы (21) и (22). Пусть точное решение, которое неизвестно, возрастает по времени, т.е.

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

Для неявной схемы (22) на возрастающем решении, наоборот, решение завышено по сравнению с точным, поскольку оно определяется по значениям сеточной функции на верхнем временном слое.

На убывающем решении картина изменяется противоположным образом: явная конечно-разностная схема завышает решения, а неявная — занижает (Рисунок 4).

На основе этого анализа возникла идея о построении более точной неявно-явной конечно-разностной схемы с весами при пространственных конечно-разностных операторах, причем при измельчении шагов тик точное (неизвестное) решение может быть взято в «вилку» сколь угодно узкую, так как если явная и неявная схемы аппроксимируют дифференциальную задачу и эти схемы устойчивы, то при стремлении сеточных характеристик τ и h к нулю решения по явной и неявной схемам стремятся к точному решению с разных сторон.

Рисунок 4 – Двусторонний метод аппроксимации

Проведенный анализ дал блестящий пример так называемых двусторонних методов, исследованных В. К. Саульевым

Рассмотрим неявно-явную схему с весами для простейшего уравнения теплопроводности:

(23)

где θ – вес неявной части конечно-разностной схемы,

θ-1 – вес для явной части

Причем

. При θ=1 имеем полностью неявную схему, при θ=0 – полностью явную схему, а при θ=1/2 – схему Кранка-Николсона.

В соответствии с гармоническим анализом для схемы (23) получаем неравенство

,

откуда

(24)

причем правое неравенство выполнено всегда.

Левое неравенство имеет место для любых значений σ, если

. Если же вес θ лежит в пределах , то между σ и θ из левого неравенства устанавливается связь

(25)

являющаяся условием устойчивости неявно-явной схемы с весами (23), когда вес находится в пределах

.

Таким образом, неявно-явная схема с весами абсолютно устойчива при

и условно устойчива с условием (25) при .

Рассмотрим порядок аппроксимации неявно-явной схемы с весами, для чего разложим в ряд Тейлора в окрестности узла (x^j,t_k) на точном решении значения сеточных функций

по переменной t, , по переменной х и полученные разложения подставим в (23):

В этом выражении дифференциальный оператор

от квадратной скобки в соответствии с дифференциальным уравнением равен дифференциальному оператору , в соответствии с чем вышеприведенное равенство приобретает вид

После упрощения получаем

,

откуда видно, что для схемы Кранка-Николсона (θ = 1/2) порядок аппроксимации схемы (23) составляет

, т.е. на один порядок по времени выше, чем для обычных явных или неявных схем. Таким образом, схема Кранка-Николсона при θ = 1/2 абсолютно устойчива и имеет второй порядок аппроксимации по времени и пространственной переменной х.

Используем в уравнение (23) подстановку r=a^2k/h². Но в то же время его нужно решить для трех "еще не вычисленных" значений

, , и . Это возможно, если все значения перенести в левую часть уравнения. Затем упорядочим члены уравнения (23) и в результате получим неявную разностную формулу

(26)

для i=2,3,…,n-1. Члены в правой части формулы (26) известны. Таким образом, формула (26) имеет вид линейной трехдиагональной системы АХ=В. Шесть точек, используемых в формуле Кранка-Николсона (26), вместе с промежуточной точкой решетки, на которой основаны численные приближения, показаны на рисунке 5.

Рисунок 5 – Шаблон (схема) метода Кранка-Николсона

Иногда в формуле (26) используется значение r=1. В этом случае приращение по оси t равно

, формула (26) упрощается и принимает вид

, (27)

для i=2,3,…,n-1. Граничные условия используются в первом и последнем уравнениях (т. е. в

и соответственно).

Уравнения (27) особенно привлекательны при записи в форме трехдиагональной матрицы АХ = В.

Если метод Кранка-Николсона реализуется на компьютере, то линейную систему АХ = В можно решить либо прямым методом, либо итерационным.

2. Практическая часть

2.1 Постановка задачи

Используем метод Кранка-Николсона, чтобы решить уравнение

,

где xϵ(0;1),

tϵ(0;0.1),

с начальным условием

,

где t=0,

xϵ(0;1),

и граничными условиями

u(0,t) = g₁(t) ≡ 0,

u(1,t) = g₂(t) ≡ 0.

2.2 Решение в ППП MatLab

Решение будем искать в ППП MatLab 7. Создадим четыре выполняемых m-фала: crnich.m – файл-функция с реализацией метода Кранка-Николсона; trisys.m – файл-функция метода прогонки; f.m – файл-функция задающая начальное условие задачи; fе.m – файл-функция задающая функцию определяющую точное решение задачи(найдена аналитическим путем). Листинги программ представлены в приложении А.

Для простоты возьмем шаг Δх = h = 0,1 и Δt = к = 0,01. Таким образом, соотношение r =1. Пусть решетка имеет n=11 столбцов в ширину и m=11 рядов в высоту.

2.3 Анализ результатов

Решения для данных параметров отразим в таблице 1. Трехмерное изображение данных из таблицы покажем на рисунке 5.

Таблица 1 – Значения u(х_i, t_i), полученные методом Кранка-Николсона

xi	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9	1
ti
0	0	1.1180	1.5388	1.1180	0.3633	0	0.3633	1.1180	1.5388	1.1180	0
0.01	0	0.6169	0.9288	0.8621	0.6177	0.4905	0.6177	0.8621	0.9288	0.6169	0
0.02	0	0.3942	0.6480	0.7186	0.6800	0.6488	0.6800	0.7186	0.6480	0.3942	0
0.03	0	0.2887	0.5067	0.6253	0.6665	0.6733	0.6665	0.6253	0.5067	0.2887	0
0.04	0	0.2331	0.4258	0.5560	0.6251	0.6458	0.6251	0.5560	0.4258	0.2331	0
0.05	0	0.1995	0.3720	0.4996	0.5754	0.6002	0.5754	0.4996	0.3720	0.1995	0
0.06	0	0.1759	0.3315	0.4511	0.5253	0.5504	0.5253	0.4511	0.3315	0.1759	0
0.07	0	0.1574	0.2981	0.4082	0.4778	0.5015	0.4778	0.4082	0.2981	0.1574	0
0.08	0	0.1419	0.2693	0.3698	0.4338	0.4558	0.4338	0.3698	0.2697	0.1419	0
0.09	0	0.183	0.2437	0.3351	0.3936	0.4137	0.3936	0.3351	0.2437	0.1283	0
0.1	0	0.1161	0.2208	0.3038	0.3570	0.3753	0.3570	0.3038	0.2208	0.1161	0

Величины, полученные методом Кранка-Николсона, достаточно близки к