Неявная разностная схема

В неявной схеме для получения решения на последующем слое необходимо решать систему алгебраических уравнений специального вида. Преимущество использования неявных схем заключается в существенном ослаблении требований к шагам сетки для выполнения условия устойчивости.

Приведем (4) к виду, удобному для применения метода прогонки:

Граничное условие:

Таким образом, построена неявная разностная схема аппроксимирующая краевую задачу с погрешностью аппроксимации порядка .

Полученная система линейных алгебраических уравнений (5)-(8) описывается трехдиагональной матрицей:

и в общем случае имеет вид:

где

Для решения таких систем применяется метод прогонки.

Прогоночные коэффициенты:

a=1;

c= ;

b=1;

с0= - из аппроксимации граничного условия.

Алгоритм решения системы состоит из двух этапов: прямого и обратного хода прогонки.

- правые части уравнений системы.

Обозначим , тогда из первого уравнения системы следует:

Подставим во второе уравнение системы при i=1 и выразим :

Продолжая подстановку далее, получим на к-ом шаге уравнение

, к =1,2,…М-1

где и ;

причём

В результате прямого хода прогонки рекуррентно вычисляются прогоночные коэффициенты и . Далее по известному коэффициенту из (12) определяются и , а затем находятся остальные . Это обратный ход прогонки.

Устойчивость и корректность метода прогонки обеспечивается при условии выполнения следующей теоремы:

если коэффициенты системы уравнений метода прогонки удовлетворяют следующие условия

причем хотя бы одно из неравенств {1} или {2} является строгим, тогда для метода прогонки имеют место неравенства:

которые гарантируют корректность и устойчивость метода прогонки.

Реализация в Matlab

U_n=zeros(N+1,K+1);

U_n(:,1)=2./(9/4*X.^2+1);

U_n(N+1,:)=8./(13+9.*T);

a=1;

c=2+h^2/tau;

b=1;

Delta=zeros(1,N-1);

Lambda=zeros(1,N-1);

for j=1:K

f_progon(1:N-1)=(h^2)*f(2:N,j)+(h^2)/tau*U_n(2:N,j);

%для предпоследнего элемента

f_progon(N-1)=f_progon(N-1)+U_n(N+1,j+1);

%для 2 элемента

znam=1+h^2/(2*tau);

f_progon(1)=f_progon(1)+(h^2/(2*tau)*U_n(1,j)+h^2/2*f(1,j+1))/znam;

c0=c-1/znam;

%%%%%%%%%

Delta(1)=b/c0;

Lambda(1)=f_progon(1)/c0;

for i=2:N-1

Delta(i)=b/(c-a*Delta(i-1));

Lambda(i)=(f_progon(i)+a*Lambda(i-1))/(c-a*Delta(i-1));

end

U_n(N,j+1)=Lambda(N-1);

for i=N-1:-1:2

U_n(i,j+1)=Delta(i-1)*U_n(i+1,j+1)+Lambda(i-1);

end

U_n(1,j+1)=(U_n(2,j+1)+h^2/(2*tau)*U_n(1,j)+h^2/2*f(1,j+1))/znam;

end;

Результат работы: