ДУ в частных производных гиперболического типа, неявная КРС

[evg_ever]

Условие задачи:

$\frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial^2 x}+2\frac{\partial^2 u}{\partial x^2}-2u$ (1)

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

$u(0,t)=\cos(2t)$ (2)

$u(1.57,t)=0$ (3)

Начальные условия:

$u(x,0)=e^{-x}\cos(x)$ (4)

$\frac{\partial u}{\partial t}(x,0)=0$ (5)

$0 \leqslant x \leqslant 1.57$

$0 \leqslant t \leqslant t_{\max}$

Нужно решить методом конечных разностей. $t_{\max}$, а также шаги по пространственной и временной координатам взять произвольные. Также надо сравнить на графиках полученное решение с $u(x,t)=e^{-x}\cos(x)\cos(2t)$, эта функция является решением уравнения.

Я вроде как решил, но графики совсем разные получаются, ошибку не могу никак найти. Решал с помощью неявной конечно-разностной схемы.

Шаг по каждой из переменных взял 0.157, $t_{\max}=1.57$. Ввёл сетку:

$x_j=0.157j$

$t^k=0.157k$

$0 \leqslant j,k \leqslant 10$

Обозначил за $u_j^k$ значение $u(x_j,t^k)$. Значения $u_j^0$ вычисляются из (4), $u_0^0,u_1^0,...,y_{10}^0$ равны:

1, 0.84419, 0.6948, 0.55639, 0.43184, 0.32265, 0.2293, 0.15144, 0.088178, 0.03825, 0.0001

$u_j^1=u_j^0$ из (5). Чтобы найти значение функции на следующих временных слоях, надо аппроксимировать значения производных в уравнении (1), здесь h_x и h_t - это шаги по переменным, равные 0.157:

$\frac{u^j_k-2u^j_{k-1}+u^j_{k-2}}{{h_t}^2}=\frac{u_{j+1}^k-2u_j^k+u_{j-1}^k}{{h_x}^2}+2\frac{u_{j+1}^k-u_j^k}{h_x}-2u_j^k$ (6)

Уравнения (2)-(3) подставляются значениями сетки очевидно, тогда $u_0^2, u_0^3, ..., u_0^{10}: 0.8092, 0.58817, 0.30962, 0.0008, -0.30811, -0.58688, -0.80827, -0.95061, -1. $
$u_{10}^k=0$.

Тогда относительно переменных $u_j^k$ на каждом временном слое, начиная с k=2, получаем трёхдиагональную СЛАУ, которую можно решить методом прогонки. Коэффициенты уравнений для временного слоя k:

$a_1=c_1=a_{10}=c_{10}=0; b_1=b_{10}=1$ - эти коэффициенты очевидно вытекают из (2)-(3).

Остальные коэффициенты берём из (6):

$a_j=\frac1{{h_x}^2}$

$b_j=-\frac2{{h_x}^2}-\frac2{h_x}-2-\frac1{{h_t}^2}$

$c_j=\frac1{{h_x}^2}+\frac2{h_x}$

$d_j=\frac{-2u_j^{k-1}-u_j^{k-2}}{{h_t}^2}$

Итоговые результаты вычисления ($u_j^k$ находится в j-той строке и k-том столбце):



В то время как решение уравнения такое:



Как видно, совсем не то, что надо И ещё видно, что моё решение сильно возрастает с увеличением t.

Код GNU/Octave (свободный аналог Matlab) для решения уравнения:

nx=11; % количество шагов по Ox
nt=11; % количество шагов по Ot
hx=0.157; % шаг по Ox
ht=0.157; % шаг по Ot
disp('начальное условие');
for x=1:nx;
    u(x,1)=u(x,2)=exp(-(x-1)*hx)*cos((x-1)*hx); % начальное условие
    disp('x');
    disp(x);
    disp('u(x,1)');
    disp(u(x,1));
end;
disp('граничные условия');
for t=3:nt;
    fl(t)=cos(2*(t-1)*ht); % левое граничное условие
    fr(t)=0; % правое граничное условие
end;
disp('левое граничное условие');
disp(fl);
disp('правое граничное условие');
disp(fr);
for t=3:nt;
    % нахождение коэффициентов a, b, c, d для прогонки
    a(1)=c(1)=a(nx)=c(nx)=0;
    b(1)=b(nx)=1;
    d(1)=fl(t);
    d(nx)=fr(t);
    for i=2:nx-1;
        a(i)=1/hx^2;
        b(i)=-2/hx^2-2/hx-2-1/ht^2;
        c(i)=1/hx^2+2/hx;
        d(i)=(-2*u(i,t-1)-u(i,t-2))/ht^2;
    end;
    pa(1)=-c(1)/b(1); % коэффициент прогонки A
    pb(1)=d(1)/b(1); % коэффициент прогонки B
    % нахождение остальных коэффициентов прогонки
    disp('a');
    disp(a);
    disp('b');
    disp(b);
    disp('c');
    disp(c);
    disp('d');
    disp(d);
    for i=2:nx;
        pa(i)=-c(i)/(b(i)+a(i)*pa(i-1));
        pb(i)=(d(i)-a(i)*pb(i-1))/(b(i)+a(i)*pa(i-1));
    end;
    disp('pa');
    disp(pa);
    disp('pb');
    disp(pb);
    % обратный шаг: нахождение u(x,t)
    u(nx,t)=pb(nx);
    for i=1:nx-1;
        u(nx-i,t)=pa(nx-i)*u(nx-i+1,t)+pb(nx-i);
    end;
    prov(1)=d(1)-b(1)*u(1,t)-c(1)*u(2,t);
    prov(nx)=d(nx)-a(nx)*u(nx-1,t)-b(nx)*u(nx,t);
    for i=2:nx-1;
        prov(i)=d(i)-a(i)*u(i-1,t)-b(i)*u(i,t)-c(i)*u(i+1,t);
    end;
    disp('prov');
    disp(prov);
end;
for t=1:nt;
    for x=1:nx;
        g(x,t)=exp(-(x-1)*hx)*cos((x-1)*hx)*cos(2*(t-1)*ht);
    end;
end;
disp(g);

[a.k.]

Как я понял, уравнение должно быть таким: $\frac{\partial^2 u}{\partial t^2}=\frac{\partial^2 u}{\partial^2 x}+2\frac{\partial u}{\partial x}-2u$

[evg_ever]

Ну да, оно ведь и написано в (1) в моём первом посте в самом верху. :) Есть идеи?

[a.k.]

evg_ever в сообщении #780653 писал(а):

Ну да, оно ведь и написано в (1) в моём первом посте в самом верху. :) Есть идеи?

В (1) в данный момент написано что-то, противоречащее здравому смыслу. А ваше уравнение заменой $u=e^{-x}v$ сводится к телеграфному уравнению $v_{xx}=v_{tt}+3v$, которое решается в явном аналитическом виде, см. скажем В.А.Ильин, Е.И.Моисеев, ДАН, 2002, Т. 387, №5, с. 600-603. Не вижу смысла здесь писать разностную схему, а тем более выискивать ошибки в коде...

[evg_ever]

a.k. в сообщении #780972 писал(а):

В (1) в данный момент написано что-то, противоречащее здравому смыслу.

Ой, ошибся Да, вы правильно написали, что должно там быть в уравнении.

a.k. в сообщении #780972 писал(а):

А ваше уравнение заменой $u=e^{-x}v$ сводится к телеграфному уравнению $v_{xx}=v_{tt}+3v$, которое решается в явном аналитическом виде, см. скажем В.А.Ильин, Е.И.Моисеев, ДАН, 2002, Т. 387, №5, с. 600-603. Не вижу смысла здесь писать разностную схему, а тем более выискивать ошибки в коде...

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

Код я привёл на всякий случай, вдруг кто-то там найдёт ошибки. Но по идее там должно быть всё нормально, а вот с логикой ошибка где-то может быть, можете помочь?