Решение уравнения теплопроводности в MatLab

[Kafari]

Всем доброго времени суток!
У меня такая задача - решить численным методом уравнение теплопроводности в MatLab и как-то его изобразить (графически). Даны граничные условия, начальная температура. То есть вот у меня есть стержень длины L при какой-то начальной температуре $T_0=T(x,0)$, в какой-то момент времени его начинают греть с двух концов с постоянными температурами $T_1=T(0,t)$ и $T_2=T(L,t)$. Собственно, нужно найти зависимость $T(x,t)$. Само уравнение такое:
$\frac{\partial T}{\partial t} = A \frac{\partial^2 T}{\partial x^2}$
Если заменить производные на разностные отношения (как нам говорили, так нужно делать), то получается следующее:

$$\frac{T(x_{i},t_{j+1})-T(x_i,t_j)}{\tau} = \frac A {h^2}(T(x_{i+1},t_{j+1})-2T(x_{i},t_{j+1})+T(x_{i-1},t_{j+1}))$$

.

, где h - интервал по х, $\tau$ - интервал по t.
И что с этим дальше делать? Как и что из этого нужно выразить? Подскажите, пожалуйста, хоть какой дальше план работы... С рисованием я как-нибудь справлюсь, но по решению чего-либо численными методами опыта пока нет.

[ewert]

Kafari в сообщении #363753 писал(а):

$$\frac{T(x_{i},t_{j+1})-T(x_i,t_j)}{\tau} = \frac A {h^2}(T(x_{i+1},t_{j+1})-2T(x_{i},t_{j+1})+T(x_{i-1},t_{j+1}))$$

.

, где h - интервал по х, $\tau$ - интервал по t.
И что с этим дальше делать? Как и что из этого нужно выразить?

Это -- "неявная" разностная схема. Выразите каждое узловое значение на старом временном слое ($t_{j}$) через комбинацию соотв. значений на новом временном слое ($t_{j+1}$). Получится некоторая система линейных уравнений для столбца новых узловых значений. Выпишите матрицу этой системы (она будет трёхдиагональной) -- и вперёд.

Это если не думать (а в учебной задаче, скорее всего, так и предполагается). Если же по-хорошему, то трёхдиагональную систему следовало бы решать не стандартными средствами Матлаба (умножением на обратную матрицу), а методом прогонки. Только соотв. процедура в Матлабе мне как-то не попадалась, а писать её самому -- хоть и небольшая, но морока (да и неизвестно, насколько она эффективно реализуется потом транслятором).

[panika_]

Доброго Вам времени суток Kafari, а не смогли бы Вы скинуть программу, которую написали))))))
Я была бы Вам ооооочень благодарна)))))

[Kafari]

panika_, программы еще пока нет)) Она в процессе. Но когда родится, выложу сюда на проверку

ewert в сообщении #364397 писал(а):

Выразите каждое узловое значение на старом временном слое ($t_{j}$) через комбинацию соотв. значений на новом временном слое ($t_{j+1}$). Получится некоторая система линейных уравнений для столбца новых узловых значений.

Видимо, я туплю, но у где взять другие уравнения для системы? У меня получается

$$T(x_i,t_j)=T(x_{i},t_{j+1})-\frac {A*\tau} {h^2} (T(x_{i+1},t_{j+1})-2T(x_{i},t_{j+1})+T(x_{i-1},t_{j+1}))$$

Но мне ведь, кажется, нужно выразить наоборот, температуру в следующий момент через предыдущий.
Поэтому появилась идея записать произодную по времени не как $\frac{T(x_{i},t_{j+1})-T(x_i,t_j)}{\tau}$, а как $\frac{T(x_{i},t_{j+2})-T(x_i,t_{j+1})}{\tau}$.
Тогда вроде получается выразить. Так нельзя делать?

[мат-ламер]

Цитата:

Но мне ведь, кажется, нужно выразить наоборот, температуру в следующий момент через предыдущий.

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

[panika_]

вот такая вот програмка у меня получилась, но там ошибки есть...исправлять надо....
задачу решаем методом прогонки....обратной и прямой....
прямая: находим прогоночные коэфф-ты...
обратная: зная Т=431 граничную(правая граница...Temp2=20*(N^2+N)), находим Т в предыдущей точке(в проге Temperature(L-1,s+1))....прогоночные коэфф-ты знаем....находим Т в L-2, и т.д......т.е. зная Т справа при L=10 находим T в L=9 и т.д......вот так вот.....


L = 10; %%Максимальное расстояние
T = 10; %Максимальное время
N=11; %Просто число...номер варианта=N
kappa=1;
%Устанавливаем шаг
h = L./(N-1); %по координате
tau = 0.01; %по времени

x = [0:h:L];
s = [0:tau:T];
% Инициализируем сеточную функцию и обнуляем её
for n=1:L+1,
for s=1:T+1,
Temperature(n,s)=0;
end;
end;
%Граничные условия
Temp1=N*20;
Temp2=20*(N^2+N);
%Начальные условия
Temperature=zeros(length(s), length(x));
Temperature(1, :) = (x>=0).*(x<L/2)*Temp1 + (x>=L/2).*(x<=L)*Temp2;
Temperature(:, 1) = 220;
Temperature(:, 10) = 341;

%Задаём коэффициенты для метода прогонки
A=kappa./h.^2;;
C=1./tau-2*kappa./h.^2;
B=kappa./h.^2;
% Устанавливаем начальные прогоночные коэффициенты
F = (1/tau)*Temperature(2,2);
alpha(3) = -kappa./h.^2;;
betta(3)=-Temp1*kappa./h.^2-(1./kappa)*Temperature(2,s);

%Запускаем обсчет по всем временам s
for s=1:length(s)
for n=4:L-1 %прямая прогонка
F = (1/tau)*Temperature(n,s);
%Прогоночные коэффициенты
alpha(n+1)=B/(C - alpha(n)*A);
betta(n+1)=(A*betta(n) + F)/(C - alpha(n)*A);
end;

Temperature(L,s+1) = Temp2; %Правая граница
Temperature(L-1,s+1)=-((betta(L-1).*A+A*Temp2+(Temperature(n+1,s)./tau)))./((A*alpha(L-1)+2*A-1./tau));
for n=L-1:-1:3, %обратная погонка
Temperature(n-2,s+1)=alpha(n-1)*Temperature(L-1,s+1)+betta(n-1)
end;
end;

plot(s, Temperature(length(s)),'-r') % Здесь надо T от времени нарисовать


ну просто аааагромная просьба....если Вы исправите эту прогу, выкиньте её сюда)))

-- Пт окт 22, 2010 22:17:48 --

Решение явной схемой у меня есть, и аналитическое решение тоже есть.....давайте напишем програмку методом прогонки)))))

[Kafari]

Как сложно...
Я не знаю метод прогонки, ни прямой, ни обратной, я просто выразила температуру в следующий момент времени через предыдущий. Разве так нельзя?
И у меня все вышло намного короче:

function [  ] = Temp( T1,T2,T0,A,L )

s=figure;                              \\ создаем окошко
step=0.05;                           \\ шаг по времени
t=0;                                     \\ начальное время
h=1;                                    \\ шаг по координате
T=T0*ones(L,1);                   \\ создаем вектор температуры от координаты
colormap(bone);                    \\ устанавливаем цветовую гамму
while ishandle(s)                    \\ цикл по времени:
    T(1)=T1;                          \\ первая стенка
    T(L)=T2;                          \\ вторая стенка
    for x=2:h:(L-1)                 \\ цикл по координате
        T(x)=(A*step/h^2)*(T(x-1)+T(x+1)-2*T(x)) - T(x-1);    \\собственно, уравнение.
    end
    image(T);                         \\ визуализируем
    t=t+step;
    pause(step);
end

end

Программа работает, но рисует что-то странное. В чем дело, пока не пойму, но никаких логических нарушений тут не вижу...

[ewert]

Kafari в сообщении #365160 писал(а):

    for x=2:h:(L-1)                 \\ цикл по координате
        T(x)=(A*step/h^2)*(T(x-1)+T(x+1)-2*T(x)) - T(x-1); 

Это и не явная схема, и не неявная -- это метод релаксации. И он тоже далеко не при всех step'ах устойчив.

(И, кстати, не имеет отношения к нестационарной задаче теплопроводности. И вообще это какая-то ерунда -- Ва зачем-то отождествили временные и пространственные индексы.)

[Kafari]

Не, я не отождествляю индексы. Просто у меня есть массив (т.е. вектор), в котором хранятся значения температур для каждой координаты. Координата - это номер элемента в массиве. Их я отождествляю, да. А в следующий момент времени просто пересчитываю этот массив, у меня получается новый, с новыми значениями температуры от координаты. Почему так нельзя?

[panika_]

Kafari, странный какой-то график....вообще должна рисоваться "ступенька", скругленная по углам....
в моей программе проблемы только с размерами массивов, я в Matlabе недавно программировать начала, поэтому найти ошибку не могу.....а вот алгоритм там правильный.....

[Kafari]

panika_, у меня программа не рисует график. Она просто визуализирует процесс, и по времени он не ограничен (то есть пока не закрыто окно.) Какая ступенька? У меня в идеале должна быть полоска, ну или прямоугльник, слева одного цвета, справа другого, и между ними установившийся градиент.
А метод прогонки я вообще не знаю, что это? Это для решения дифференциального уравнения? И неужели без него никак нельзя?

[panika_]

можно конечно...метод прогонки это один из численных методов, используемый для решения уравнения теплопроводности....и график получается ступенькой

[ewert]

Kafari в сообщении #365367 писал(а):

Не, я не отождествляю индексы.

Привет. А что это за "-T(x-1)" в самом конце того оператора? Что здесь-то означает "x"?... (и, кстати, откуда минус-то?...)

-- Вс окт 24, 2010 16:44:09 --

Kafari в сообщении #365607 писал(а):

А метод прогонки я вообще не знаю, что это?

Метод прогонки в данном случае -- это просто обычный метод Гаусса с обратным ходом. В случае, когда матрица системы трёхдиагональна (т.е. когда отличны от нуля элементы только на главной диагонали и на двух соседних), он реализуется очень быстро и с помощью вполне явных рекуррентных формул (общее количество операций -- порядка всего лишь $8n$).

[Kafari]

ewert в сообщении #365690 писал(а):

А что это за "-T(x-1)" в самом конце того оператора?

Да, там должен быть плюс. Это у меня был глюк. Тогда так:

T(x)=(A*step/h^2)*(T(x-1)+T(x+1)-2*T(x)) + T(x);

"х" означает координату.
Это тоже неправильно?
С другой стороны, мне сейчас показалось, что там может возникнуть одна ошибка... Когда на следующем шаге вычисляется температура в следующей координате, то в выражении используется "T(x-1)", для которого уже посчитано новое значение. А мне нужно старое. Из-за этого все плохо?

Да, и все же... Где взять другие уравнения для системы? Из начальных условий? А то я что-то в упор не вижу...

[ewert]

Kafari в сообщении #366644 писал(а):

Когда на следующем шаге вычисляется температура в следующей координате, то в выражении используется "T(x-1)", для которого уже посчитано новое значение. А мне нужно старое. Из-за этого все плохо?

Совершенно верно. Из-за этого и выходит метод релаксации, который отношения к нестационарной теплопроводности не имеет.

Или левая часть оператора присваивания, или последнее слагаемое справа должны относиться к другому массиву, нежели все остальные слагаемые (потом, после прогонки этого цикла, вновь пересчитанные значения можно переместить обратно в "старый" массив, если предыдущие хранить не требуется). В зависимости от варианта и получается или явная, или неявная схема для задачи теплопроводности.