Одномерное уравнение теплопроводности. Численное решение.

[rudolfninja]

Всем доброго времени суток.
Есть задание: "Решить одномерное нестационарное уравнение (см. формулу ниже) теплопроводности методом сеток."

$\frac{\partial u}{\partial t}=\frac{d}{dx}\left( g\frac{\partial u}{\partial x}\right) + f$

При граничных условиях:

$q^0 \frac{\partial u}{\partial x} + \left.\alpha^0u\right|_{x=0}=\beta^0$

$q^1 \frac{\partial u}{\partial x} + \left.\alpha^1u\right|_{x=1}=\beta^1$

$u(x,0)=u^0(x)$

Область интегрирования $\Omega=\left\lbrace0\leqslant x\leqslant b, 0\leqslant t\leqslant T\right\rbrace$ покрываем равномерной сеткой $\omega_{h,r}=\left\lbrace(i-1)h, k\tau, i=1...N, k=0...K\right\rbrace$, $\tau - \text {шаг по времени}; h - \text {шаг по }x; k - \text {номер временного слоя}; i - \text {номер точки}$. Для нахождения таблицы решения $\overline{u}_h = \left\lbrace \overline{u}^k_i \cdot u(x_it^k) \right\rbrace$ воспользуемся следующей конечно разностной схемой:

$\frac{\overline{u}^{k+1}_i - \overline{u}^k_i}{\tau}=\frac{\omega_r}{h^2} \Bigl[g_{i-1/2}\overline{u}^{k+1}_{i-1} - (g_{i-1/2} + g_{i+1/2}) \overline{u}^{k+1}_i + g_{i+1/2}\overline{u}^{k+1}_{i+1}$ \Bigl] + \\ + \frac{1 - \omega_r}{h^2} \Bigl[g_{i-1/2}\overline{u}^k_{i-1} - (g_{i-1/2} + g_{i+1/2}) \overline{u}^k_i + g_{i+1/2}\overline{u}^k_{i+1} \Bigl] + \overline{f}_i; \omega_r = 0.5$

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

$q^0\frac{u_2 - u_1}{h} + \alpha^0 u_1 = \beta^0$,

$q^1\frac{u_{n+1} - u_n}{h} + \alpha^1 u_{n+1} = \beta^1$

Явная схема ($\omega_r=0$):

$u_1 = (\beta^0 - q^0/h \cdot u_2)/(\alpha^0 - q^0/h) $
$u_{n+1} = (\beta^1 + q^1/h \cdot u_n)/(\alpha^1 + q^1/h) $

Неявная схема ($\omega_r=0$):

$b_1u_1 + c_1u_2 = d_1; \Rightarrow d_1 = \beta^0; b_1=\alpha^0 - q^0/h; c_1=q^0/h;$

$a_{n+1}u_n + b_{n+1}u_{n+1} = d_{n+1}; \Rightarrow d_{n+1} = \beta^1; b_{n+1}=\alpha^1- q^1/h; \alpha_{n+1}=-q^1/h.$

Расчет выполнять до достижения сходимости, т.е. пока

$\max\limits_{i}\left\lvert \overline{u}^k_i - \overline{u}^{k+1}_i \right\rvert \leqslant\varepsilon$

Выдать графики изменения $u(x,t^k), t^k=0,m\tau,2m\tau,3m\tau,\ldots u(x,tk), tk=0, m$

Тест: при $f=1$ и $g=1$ точное решение в установившемся режиме

$u=- \frac{x^2}{2} + c_1x + c_2$

Подберите $c_1$ и $c_2$ такими, чтобы выполнялись ваши граничные условия и протестируйте программу для этого случая. Подберите значения $N, k, m$ такими, чтобы в установившемся режиме у вас получилось точное решение. После этого подставляйте нужные значения $f$ и $g$.

Затем идет пример решения для следущих условий:
$\omega_r={0,1}$
$g=1+x^2$
$f=1-x^2$
$q^0=0$
$\alpha^0=1$
$\beta^0=0$
$q^1=0$
$\alpha^1=1$
$\beta^1=0$
Пример программы ($\omega_r=0.5$ или $\omega_r=1$)

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

function V1_16_1(N,K,m,wr)
%решение диф ур по неявной схеме
%нестационарное уравнение типа теплопроводности
h=1/N; h2=h^2;
tau=h/2;
q0=0; al0=1; be0=0; q1=0; al1=1; be1=0;
for i=1:N+1
    u(i)=0; u1(i)=0; x(i)=(i-1)*h; x1(i)=(i-0.5)*h;
    f(i)=1-x(i)^2;
    g1(i)=1+x1(i)^2;
end;
 t=0;
 plot(x,u);
 %disp('t='); disp(t);
for kt=1:K
    t=t+tau;
for i=2:N
    a(i)=wr*tau/h2*g1(i-1);
    c(i)=wr*tau/h2*g1(i);
    b(i)=-a(i)-c(i)-1;
    ad=(1-wr)*tau/h2*g1(i-1);
    cd=(1-wr)*tau/h2*g1(i);
    bd=-ad-cd+1;
    d(i)=-ad*u(i-1)-bd*u(i)-cd*u(i+1)-tau*f(i);
end
c(1)=q0/h; b(1)=al0-q0/h; d(1)=be0;
a(N+1)=-q1/h; b(N+1)=al1+q1/h; d(N+1)=be1; %c(N+1)=0;
%u=mprog(a,b,c,d,N);
ks(1)=-c(1)/b(1);
et(1)=d(1)/b(1);
for i=2:N
    z=b(i)+a(i)*ks(i-1);
    ks(i)=-c(i)/z;
    et(i)=(d(i)-a(i)*et(i-1))/z;
end
u1(N+1)=(d(N+1)-a(N+1)*et(N))/(b(N+1)+a(N+1)*ks(N));
for i=N:-1:1
    u1(i)=ks(i)*u1(i+1)+et(i);
end
em=0;
for i=1:N+1
        e=abs(u1(i)-u(i));
            u(i)=u1(i);
     if e>em
     em=e;
    end
end
if mod(kt,m)==0
   hold on
   pause(0.1);
   plot(x,u);
  %t
   %disp(u);
end
%if em<0.001
%    break
%end
end
hold off
return

Пример программы ($\omega_r=0$)

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

function V1_16_0(N,K,m)
%решение диф ур по явной схеме
%нестационарное уравнение типа теплопроводности
h=1/N; h2=h^2;
tau=h2/4;
q0=0; al0=1; be0=0; q1=0; al1=1; be1=0;
for i=1:N+1
    u(i)=0; u1(i)=0; x(i)=(i-1)*h; x1(i)=(i-0.5)*h;
    f(i)=1-x(i)^2;
    g1(i)=1+x1(i)^2;
end;
u(1)=(be0-q0/h*u(2))/(al0-q0/h);
u(N+1)=(be1+q1/h*u(N))/(al1+q1/h);
plot(x,u);
 t=0;
 %disp('t='); disp(t);
for kt=1:K
    t=t+tau;
    em=0;
for i=2:N
    u1t=u(i)+tau/h^2*(g1(i-1)*u(i-1)-(g1(i-1) + g1(i)) * u(i) + g1(i) * u(i+1)) + tau*f(i);
    e=abs(u1t-u1(i));
if  e>em;
    em=e;
end
u1(i)=u1t;
end
for i=2:N
    u(i)=u1(i);
end
u(1)=(be0-q0/h*u(2))/(al0-q0/h);
u(N+1)=(be1+q1/h*u(N))/(al1+q1/h);
if mod(kt ,m)==0
    hold on
    pause(0.5);
    plot(x,u)
    disp(u);
   % if em<0.001
   %     break
   %end
end
end
hold off
 
return
 

Вопрос, собсветнно в чем: помогите, пожалуйста, разобраться что за переменная $\omega_r$, откуда берутся такие значения шага, откуда берутся такие значения векторов $u$ и $u1$. В общем, пожалуйста, прокомментируйте код, в соответсвии с тем, что написано в прилогаемом материале. А то как я не пытался понять почему код написан именно так, у меня не получается.

Заранее спасибо!

[Lia]

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

rudolfninja
Просьба изложить основную задачу полностью здесь, с помощью $\TeX$. Таблицу тоже. Остальные формулы тоже.
Сопроводительные материалы можно оставить на картинке.

- неправильно набраны формулы (краткие инструкции: «Краткий FAQ по тегу [math]» и видеоролик Как записывать формулы);


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

[Lia]

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

[мат-ламер]

rudolfninja в сообщении #1070180 писал(а):

Вопрос, собсветнно в чем: помогите, пожалуйста, разобраться что за переменная $\omega_r$

ИМХО, она обеспечивает плавный переход от неявной к схеме к явной.

[cool.phenon]

Когда значение параметра равно нулю, имеем чисто явную схему, она условно устойчива, поэтому шаг нужно подбирать исходя из условия Куранта-Фридрихса-Леви. Когда $w_r=\frac{1}{2}$, имеем схему Кранка-Николсона, она безусловно устойчива. То же самое можно сказать о $w_r=1$. Если взять любое другое значение параметра, то решение может развалиться. Параметр показывает, какую из комбинаций неявного и явного метода будем использовать.

-- 11.11.2015, 02:07 --

З.Ы. Вектор $u_1$ выражен исходя из краевых условий и того, что в неявном методе используется алгоритм Томаса, который у нас называют прогонкой.

[cool.phenon]

Утром прочитал, немного неточно ночью написал. Параметр можно брать от нуля до единицы, нельзя брать отрицательных значений или значений больше единицы

[rudolfninja]

Спасибо за разъяснения.
А можете еще рассказать про параметры, передаваемые в функцию и в соответствии с чем их выбирать (в частности N, K, m)?

[TOTAL]

rudolfninja в сообщении #1070180 писал(а):

Расчет выполнять до достижения сходимости, т.е. пока
$\max\limits_{i}\left\lvert \overline{u}^k_i - \overline{u}^{k+1}_i \right\rvert \leqslant\varepsilon$

Берите очень маленький шаг по времени. Тогда "сходимость" (согласно этому великолепному критерию) случится за одну итерацию.

[rudolfninja]

TOTAL, спасибо за совет, но меня пока что больше интересуют параметры, которыю передаются в функцию.