Всем доброго времени суток.
Есть задание: "Решить одномерное нестационарное уравнение (см. формулу ниже) теплопроводности методом сеток."
![$\frac{\partial u}{\partial t}=\frac{d}{dx}\left( g\frac{\partial u}{\partial x}\right) + f$ $\frac{\partial u}{\partial t}=\frac{d}{dx}\left( g\frac{\partial u}{\partial x}\right) + f$](https://dxdy-04.korotkov.co.uk/f/7/8/1/7817af48f9e22c6416470f2258a93ac282.png)
При граничных условиях:
![$q^0 \frac{\partial u}{\partial x} + \left.\alpha^0u\right|_{x=0}=\beta^0$ $q^0 \frac{\partial u}{\partial x} + \left.\alpha^0u\right|_{x=0}=\beta^0$](https://dxdy-04.korotkov.co.uk/f/f/7/8/f786fa98f2d81160f03b5d957a22639c82.png)
![$q^1 \frac{\partial u}{\partial x} + \left.\alpha^1u\right|_{x=1}=\beta^1$ $q^1 \frac{\partial u}{\partial x} + \left.\alpha^1u\right|_{x=1}=\beta^1$](https://dxdy-04.korotkov.co.uk/f/7/e/b/7ebcabe233f4cda087abf2a1e70c90f582.png)
![$u(x,0)=u^0(x)$ $u(x,0)=u^0(x)$](https://dxdy-03.korotkov.co.uk/f/6/9/6/69602b6daff37f5b9b77df6bff36452c82.png)
Область интегрирования
![$\Omega=\left\lbrace0\leqslant x\leqslant b, 0\leqslant t\leqslant T\right\rbrace$ $\Omega=\left\lbrace0\leqslant x\leqslant b, 0\leqslant t\leqslant T\right\rbrace$](https://dxdy-03.korotkov.co.uk/f/a/5/c/a5c32af69a5a9b3ab23800e32a05308d82.png)
покрываем равномерной сеткой
![$\omega_{h,r}=\left\lbrace(i-1)h, k\tau, i=1...N, k=0...K\right\rbrace$ $\omega_{h,r}=\left\lbrace(i-1)h, k\tau, i=1...N, k=0...K\right\rbrace$](https://dxdy-02.korotkov.co.uk/f/9/0/c/90cd563db7ae4fdf43a96fa7cd23c47d82.png)
,
![$\tau - \text {шаг по времени}; h - \text {шаг по }x; k - \text {номер временного слоя}; i - \text {номер точки}$ $\tau - \text {шаг по времени}; h - \text {шаг по }x; k - \text {номер временного слоя}; i - \text {номер точки}$](https://dxdy-03.korotkov.co.uk/f/6/5/2/6521ef85783d0f87f71f0960baa735a182.png)
. Для нахождения таблицы решения
![$\overline{u}_h = \left\lbrace \overline{u}^k_i \cdot u(x_it^k) \right\rbrace$ $\overline{u}_h = \left\lbrace \overline{u}^k_i \cdot u(x_it^k) \right\rbrace$](https://dxdy-03.korotkov.co.uk/f/e/f/c/efcefe4fe2d1916c3871656437e6a60182.png)
воспользуемся следующей конечно разностной схемой:
![$\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$ $\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$](https://dxdy-02.korotkov.co.uk/f/5/4/e/54e07d02f9e7b1f8f2ec732b0ce3136682.png)
Граничные условия аппроксимировать следующим образом:
![$q^0\frac{u_2 - u_1}{h} + \alpha^0 u_1 = \beta^0$ $q^0\frac{u_2 - u_1}{h} + \alpha^0 u_1 = \beta^0$](https://dxdy-02.korotkov.co.uk/f/1/d/8/1d8a1e6723047988a23a3b508586f79982.png)
,
![$q^1\frac{u_{n+1} - u_n}{h} + \alpha^1 u_{n+1} = \beta^1$ $q^1\frac{u_{n+1} - u_n}{h} + \alpha^1 u_{n+1} = \beta^1$](https://dxdy-02.korotkov.co.uk/f/1/3/5/1351aae5490f6a74d6d4fd933c2119f182.png)
Явная схема (
![$\omega_r=0$ $\omega_r=0$](https://dxdy-03.korotkov.co.uk/f/2/0/3/2030a5d54ffeb6b4395d02a973cb182482.png)
):
![$u_1 = (\beta^0 - q^0/h \cdot u_2)/(\alpha^0 - q^0/h) $ $u_1 = (\beta^0 - q^0/h \cdot u_2)/(\alpha^0 - q^0/h) $](https://dxdy-03.korotkov.co.uk/f/a/8/2/a829acc21e247c232122183d725d94ed82.png)
![$u_{n+1} = (\beta^1 + q^1/h \cdot u_n)/(\alpha^1 + q^1/h) $ $u_{n+1} = (\beta^1 + q^1/h \cdot u_n)/(\alpha^1 + q^1/h) $](https://dxdy-01.korotkov.co.uk/f/8/6/4/86429b4d55e0774cd3a475e23818c83b82.png)
Неявная схема (
![$\omega_r=0$ $\omega_r=0$](https://dxdy-03.korotkov.co.uk/f/2/0/3/2030a5d54ffeb6b4395d02a973cb182482.png)
):
![$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;$ $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;$](https://dxdy-03.korotkov.co.uk/f/a/f/8/af83c1e48a253beeeeca4b61f28152ea82.png)
![$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.$ $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.$](https://dxdy-03.korotkov.co.uk/f/e/a/6/ea68d5ada64a60e0660bb186fdfe315682.png)
Расчет выполнять до достижения сходимости, т.е. пока
![$\max\limits_{i}\left\lvert \overline{u}^k_i - \overline{u}^{k+1}_i \right\rvert \leqslant\varepsilon$ $\max\limits_{i}\left\lvert \overline{u}^k_i - \overline{u}^{k+1}_i \right\rvert \leqslant\varepsilon$](https://dxdy-04.korotkov.co.uk/f/b/a/9/ba9e11a7ae48f66dc7682c1e9c1b46a682.png)
Выдать графики изменения
![$u(x,t^k), t^k=0,m\tau,2m\tau,3m\tau,\ldots u(x,tk), tk=0, m$ $u(x,t^k), t^k=0,m\tau,2m\tau,3m\tau,\ldots u(x,tk), tk=0, m$](https://dxdy-04.korotkov.co.uk/f/f/6/0/f600c27b4bce282906cdf847ad87198882.png)
Тест: при
![$f=1$ $f=1$](https://dxdy-01.korotkov.co.uk/f/0/5/6/05677363f3635f772075a606e086503582.png)
и
![$g=1$ $g=1$](https://dxdy-02.korotkov.co.uk/f/1/8/8/188fdf6c12e3fe5e9e4904c24af8957f82.png)
точное решение в установившемся режиме
![$u=- \frac{x^2}{2} + c_1x + c_2$ $u=- \frac{x^2}{2} + c_1x + c_2$](https://dxdy-02.korotkov.co.uk/f/d/d/2/dd2ff997035ad54f34de18c370f8b1be82.png)
Подберите
![$c_1$ $c_1$](https://dxdy-02.korotkov.co.uk/f/9/8/8/988584bba6844388f07ea45b7132f61c82.png)
и
![$c_2$ $c_2$](https://dxdy-03.korotkov.co.uk/f/e/3/5/e355414b8774603011922d600510b1df82.png)
такими, чтобы выполнялись ваши граничные условия и протестируйте программу для этого случая. Подберите значения
![$N, k, m$ $N, k, m$](https://dxdy-02.korotkov.co.uk/f/1/0/b/10b233dec1bd970abe0bf59dee159c8082.png)
такими, чтобы в установившемся режиме у вас получилось точное решение. После этого подставляйте нужные значения
![$f$ $f$](https://dxdy-02.korotkov.co.uk/f/1/9/0/190083ef7a1625fbc75f243cffb9c96d82.png)
и
![$g$ $g$](https://dxdy-04.korotkov.co.uk/f/3/c/f/3cf4fbd05970446973fc3d9fa3fe3c4182.png)
.
Затем идет пример решения для следущих условий:
![$\omega_r={0,1}$ $\omega_r={0,1}$](https://dxdy-02.korotkov.co.uk/f/5/a/a/5aaf0799c8983b25b7198e0bf46dae2e82.png)
![$g=1+x^2$ $g=1+x^2$](https://dxdy-03.korotkov.co.uk/f/a/3/2/a3220c01687f7de8f4c70aae37dcb69a82.png)
![$f=1-x^2$ $f=1-x^2$](https://dxdy-04.korotkov.co.uk/f/f/5/5/f5518a412c2477e88f521d36536f119e82.png)
![$q^0=0$ $q^0=0$](https://dxdy-01.korotkov.co.uk/f/c/8/8/c88ba9d9a5b030702cba5f63d720914e82.png)
![$\alpha^0=1$ $\alpha^0=1$](https://dxdy-03.korotkov.co.uk/f/2/e/8/2e8012b09957023ef1a9333c3efa528882.png)
![$\beta^0=0$ $\beta^0=0$](https://dxdy-01.korotkov.co.uk/f/c/2/1/c213107dc5c3c651385147f82fea7b1d82.png)
![$q^1=0$ $q^1=0$](https://dxdy-02.korotkov.co.uk/f/5/c/d/5cd1581e7781cefc9a2673582b70e3dc82.png)
![$\alpha^1=1$ $\alpha^1=1$](https://dxdy-02.korotkov.co.uk/f/1/b/f/1bfe8aa0a19ea34e18e2796007d5364e82.png)
![$\beta^1=0$ $\beta^1=0$](https://dxdy-01.korotkov.co.uk/f/8/b/c/8bc9433dc36d0f8c5fd2244c69fa124d82.png)
Пример программы (
![$\omega_r=0.5$ $\omega_r=0.5$](https://dxdy-04.korotkov.co.uk/f/3/9/e/39ee9ed42e0af198e2b914f2cdfe0d6282.png)
или
![$\omega_r=1$ $\omega_r=1$](https://dxdy-02.korotkov.co.uk/f/5/e/f/5ef73c503284acc4c0a730a68cdf757d82.png)
)
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$ $\omega_r=0$](https://dxdy-03.korotkov.co.uk/f/2/0/3/2030a5d54ffeb6b4395d02a973cb182482.png)
)
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$ $\omega_r$](https://dxdy-01.korotkov.co.uk/f/0/0/e/00e132f3f2840c12fa39ff487aac86fc82.png)
, откуда берутся такие значения шага, откуда берутся такие значения векторов
![$u$ $u$](https://dxdy-03.korotkov.co.uk/f/6/d/b/6dbb78540bd76da3f1625782d42d6d1682.png)
и
![$u1$ $u1$](https://dxdy-01.korotkov.co.uk/f/4/8/6/486bd388d14cc06ed0b885f3a714d01182.png)
. В общем, пожалуйста, прокомментируйте код, в соответсвии с тем, что написано в прилогаемом материале. А то как я не пытался понять почему код написан именно так, у меня не получается.
Заранее спасибо!