Странное поведение неявных разностных схем

[robot80]

Здравствуйте, уважаемые участики форума!
В рамках предмета уравнения математической физики разбираю численную реализацию решения базовых УРЧП (волновое и теплопередачи). Рассматриваю явные и неяные разностные схемы. С явными схемами проблем нет в том смысле, что практическая их реализация (кодирование), расчёты, поведение решения согласуется с теорией.
А вот с неявными схемами что-то не ладится. Для уравнения теплопередачи я рассматриваю схему Крэнка-Николсон, а для волнового уравнения трёхслойную схему с весами (вес беру равным 0,5, что приводит к схеме похожей на схему Крэнка-Николсон - см. Калиткин Н.Н.). Краевые условия самые простые: $u(0,t) = \mu_1(t), u(1,t) = \mu_2(t), u(x,0) = \mu_3(x) , \frac {\partial u (x, 0)} {\partial t} = \mu_4(x) $ Я приведу волновое уравнение, разностную схему и код на Scilab. Я не могу понять, где я ошибаюсь? При реализации появляются жуткие осцилляции, которых вроде бы не должно быть.

$ \frac {\partial^2 u } {\partial t^2} =a^2 \frac {\partial^2 u} {\partial x^2} + f(t,x) $ - уравнение.
$\frac {u_j^{i-1} -2u_j^i +u_j^{i+1}}{h_t^2} = \frac{a^2}{2h_x^2} \left( {u_{j-1}^{i-1} -2u_j^{i-1} +u_{j+1}^{i-1}} +{u_{j-1}^{i+1} -2u_j^{i+1} +u_{j+1}^{i+1}} \right) +f_j^i$ - разностная схема
ДОБАВЛЕНИЕ. $a = 0.5, h_x = 0.01, h_t = 0.02, X = 1.0, T = 2.0, f(t, x) =0 $
$ \gamma u_{j-1}^{i+1} - (1 + 2 \gamma)u_j^{i+1} +\gamma u_{j+1}^{i+1} = u_j^{i-1} -2 u_j^i - \gamma \left(u_{j-1}^{i-1} -2u_j^{i-1} +u_{j+1}^{i-1} \right) - h_t^2 f_j^i $ - разностная схема приведена к виду системы уравнений из nx-1 уравнения, при условии, что пространственная сетка состоит из nx+1 узла. Система решается стандартной прогонкой.
Код в синтаксисе scilab

Используется синтаксис Matlab M

function [x,t,u] = WOLNA (T, ht, X, hx, a)
    nt = round(T./ht); nx = round(X./hx); 
    u = zeros(nt+1,nx+1); gam = zeros(nx+1);
    
           gam = 0.5.*(ht.*a./hx).^2.0;

        x = zeros(nx+1); u = zeros (nt+1, nx +1); t = zeros(nt+1);

    for j = 1 : 1: nx+1
        x(j) =  (j - 1).*hx;
        u(1,j) = mu3(x(j));
    end
    
    t(1) = 0.0;
    for i = 2 : 1 : nt+1
        t(i) = (i - 1).*ht;
        u(i,1) = mu1(t(i));
        u(i, nx+1) = mu2(t(i)); 
    end  
    
    for j = 2 : 1 : nx
      u(2,j) = u(1, j) + 0.5.*gam.*(u(1,j-1) - 2.0.*u(1,j) + u(1,j+1)) + 0.5.*ht.^2.0.*func(t(1), x(j))+ ht.*mu4(x(j));
    end
    
    for i = 2 : 1 : nt 
        A = zeros(nx-1); B = zeros(nx-1); C = zeros(nx-1); F = zeros(nx-1); xm = zeros(nx-1);
        F(1) = u(i - 1, 2) - 2.0.*u(i, 2) - ht.*ht.*func(t(i), x(2)) - ...
        gam.*(u(i+1, 1) + u(i-1, 1) -2.0.*u(i-1,2) + u(i-1, 3));
        A(1) = 0.0; C(1) = -(1.0 + 2.0.*gam); B(1) = gam;
        
        F(nx - 1) = u(i - 1, nx) - 2.0.*u(i, nx) - ht.*ht.*func(t(i), x(nx)) -...
        gam.*(u(i + 1, nx + 1) + u(i-1, nx-1) -2.0.*u(i-1,nx) + u(i-1, nx+1)); 
        
        A(nx - 1) = gam; C(nx-1) = -(1.0 + 2.0.*gam); B(nx-1) = 0.0;
        
        for j = 2 : 1 : nx-2
            F(j) = (1.0 +2.0.*gam).*u(i - 1, j) - 2.0.*u(i, j) - ht.*ht.*func(t(i), x(j))- ...
            gam.*(u(i - 1, j - 1) + u(i - 1, j + 1)); 
            A(j) = gam; C(j) = -(1.0 +2.0.*gam); B(j) = gam;
        end
        
        xm = progonka(A,B,C,F,nx-1);
        for j = 2 : 1 : nx
            u(i+1, j) = xm(j-1); 
        end
    end;  
endfunction

function z = func (t,x)
  z = 0.0;
endfunction

function [ym] = mu1(t)
    ym = 0.0;    
endfunction

function [ym] = mu2(t)
    ym = 0.0;    
endfunction

function [ym] = mu3(x)
if (x <=1/3) then
    ym = 0.03.*x;
else
    ym = 0.015.*(1.0 - x);    
end
endfunction

function [ym] = mu4(x)
     ym = 0.0;
endfunction


function [xm] = progonka(A,B,C,F,n)
    for  i = 2: 1: n
        m = A(i)./C(i-1);
        C(i) = C(i) - m.*B(i-1);
        F(i) = F(i) - m.*F(i-1);
    end
    x(n) = F(n-1)./C(n-1);
    for i = n-1 :-1 : 1
        x(i) = (F(i) -B(i)*x(i+1))./C(i);
    end
    xm = x;
endfunction 

 

ДОБАВЛЕНИЕ 2. Если имеется в виду вызов функции, то вот он

Используется синтаксис Matlab M

[x,t,u] = WOLNA (2.0, 0.02, 1.0, 0.01, 0.5);

. Если что-то другое, то поясните пожалуйста, что именно имеется в виду.

[Pphantom]

Чему у вас равна скорость волны $a$ и шаги по координате и времени $h_x$ и $h_t$? Остальное (вид источникового члена и краевые условия) можно выцедить из кода, хотя и их было бы не лишним написать непосредственно.

[Pphantom]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Карантин»
по следующим причинам:

- см.выше;
- стоит дополнить код так, чтобы его можно было запустить.

Исправьте все Ваши ошибки и сообщите об этом в теме Сообщение в карантине исправлено.
Настоятельно рекомендуется ознакомиться с темами Что такое карантин и что нужно делать, чтобы там оказаться и Правила научного форума.

[Pphantom]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»

[robot80]

Спасибо за возвращение из карантина!
В общем, главный косяк я нашёл (пришлось уравнения ситемы полностью расписать и выявить зависимость между всеми индексами).
Изменения произошли во внутреннем цикле.

Используется синтаксис Matlab M

  for j = 2 : 1 : nx-2
            F(j) = (1.0 +2.0.*gam).*u(i - 1, j+1) - 2.0.*u(i, j+1) - ht.*ht.*func(t(i), x(j+1))- ...
            gam.*(u(i - 1, j ) + u(i - 1, j + 2)); 
            A(j) = gam; C(j) = -(1.0 +2.0.*gam); B(j) = gam;
        end
 

Но всё равно. Поведение решения, полученного с помощью неявной схемы отличается от аналогичного решения, полученного с помощью явной схемы (все параметры естесственно одинаковые). Чем больше временной интервал, тем сильнее проявляется волнообразность...

[Pphantom]

Я не уверен, что помню синтаксис правильно, но вроде во всех MatLAB'оподобных штуках функция zeros(n) возвращает нулевую матрицу размером $n\times n$. У вас же везде предполагается, что получится вектор. Как это повлияет на работоспособность кода - сходу сказать трудно, я не очень привык пользоваться этими инструментами, но, наверное, лучше бы поправить на всякий случай.

[robot80]

Спасибо, исправил. В данном случае это никак не повлияет на результат, но для исключения возможных неожиданностей надо писать как-то так

Используется синтаксис Matlab M

zeros(1, nx+1)

[Pphantom]

robot80 в сообщении #1432404 писал(а):

Чем больше временной интервал, тем сильнее проявляется волнообразность...

Еще бы понять, как именно узреть эту волнообразность...

[robot80]

Вот картинка для явной схемы "крест":
А вот картинка неявной схемы
На первой картинке видно, что поверхность гладкая и в конечный момент времени чётко отзеркаливает начальное состояние.
На второй картике видна рябь на поверхности и в конечный момент времени поверхность изогнутая получается.
Видно конечно не слишком хорошо. Но если изменить парметры вызова расчётной фуцнкции, то разница будет заметнее.

[robot80]

Собственно, вопрос и проблема здесь скорее библиографическая. Может быть кто-то встречал книги или статьи, где детально исследовалось бы поведение решений, полученных явной или неявной схемами. То, что нашёл я, ограничевается стандартной устойчивостью и точностью. Может быть есть ещё какие-то условия для неявной схемы, которая бы улучшила её поведение?