Квазилинейное уравнение теплопроводности(разн. схема)

[dasalam]

Хочу я решить вот такое квазиленейное одномерное уравнение теплопроводности:
$\\ \frac{\partial U}{\partial t} = \frac{\partial }{\partial x}(U^{k} \frac{\partial U}{\partial x})\\ 0 \le t \le T; 0 \le x \le X; \\ U(x,0) = 0; U(0,t) = Ct^{1/k}; U(X, t) = 0 $
Для этого использую неявную разностную схему

, где
$\\a_{m+1} = 0.5(U^{n}_{m}+U^{n}_{m+1})\\ a_{m-1} = 0.5(U^{n}_{m}+U^{n}_{m-1}) $
Точное решение известно (можно посмотреть тут http://www.intuit.ru/department/calcula ... q/2/3.html). Но при таких НУ и ГУ получаем, что
$U^{0}_{m} = 0; U^{1}_{m}=0 ... \forall m$
Что тут неверно?

[dasalam]

Может кто подскажет какую-нибудь другую разностную схему?

[svv]

Вопросы:
1) Чтобы посмотреть точное решение на intuit.ru, там надо зарегистрироваться?
2) $U^k$ в уравнении -- это $k$-я степень? Тогда см. 5)
3) Верхний индекс $n, n+1$ соответствует $t$, а нижний $m, m+1, m-1$ соответствует $x$, так?
4) $U$ и $u$ -- это одно и то же? В чем разница?
5) Если в $a_{m+1} = 0.5(U^{n}_{m}+U^{n}_{m+1})$ индекс $n$ -- это индекс, то где степень $k$?
6) Откуда взялось $f_m^n$, если уравнение однородное?
7) Для чего в разностном уравнении в правой части верхние индексы везде $n+1$, а не $n$? Ведь эти значения на очередном шаге по $n$ еще неизвестны.

[dasalam]

svv в сообщении #496978 писал(а):

Вопросы:
1) Чтобы посмотреть точное решение на intuit.ru, там надо зарегистрироваться?
2) $U^k$ в уравнении -- это $k$-я степень? Тогда см. 5)
3) Верхний индекс $n, n+1$ соответствует $t$, а нижний $m, m+1, m-1$ соответствует $x$, так?
4) $U$ и $u$ -- это одно и то же? В чем разница?
5) Если в $a_{m+1} = 0.5(U^{n}_{m}+U^{n}_{m+1})$ индекс $n$ -- это индекс, то где степень $k$?
6) Откуда взялось $f_m^n$, если уравнение однородное?
7) Для чего в разностном уравнении в правой части верхние индексы везде $n+1$, а не $n$? Ведь эти значения на очередном шаге по $n$ еще неизвестны.

1) Да надо, но вот оно



2) Да, это степепень.
3) Да, верхний индекс - время, нижний - координата.
4) $U$ и $u$ одно и то же, просто не хотелось набирать, то что уже есть.
5) Извиняюсь, $a_{m+1} = 0.5((U^{n}_{m})^k+(U^{n}_{m+1})^k)$
6) $f_m^n = 0$
Извиняюсь за неточность
7) Да неизвестны, это неявная схема. Их надо найти как раз (я предполагаю сделать это методом прогонки).

[svv]

$u^0_m=0$ -- просто по начальному условию $u(x, 0)=0$.
Ваша физическая ситуация такова: начальный момент -- везде ледяной холод. Правая граница -- всегда ледяной холод. Левая граница медленно нагревается, и от нее разогревается вся область. Но состояние левой границы на некотором шаге $n$ сказывается на соседних пространственных точках только на следующих шагах. А так как в момент $n=0$ левая граница ещё холодная как лед ($Ct^{1/k}=0}$) и не греет, то и в момент $n=1$ вся область (кроме левой границы) ещё сохраняет исходную температуру: $u^1_m=0$ -- её ещё ничто не разогрело.

Измените в правой части верхние индексы на $n$, левую не трогайте. Тогда Ваше разностное уравнение будет давать $u^{n+1}_m$ из $u^n_{m-1}$, $u^n_m$, $u^n_{m+1}$. Схема станет явной.

Пусть $M$ -- максимальное значение нижнего индекса.
Начало. $u^0_m=0$ для всех $m$.
Шаг итерации. Пусть Вы уже вычислили $u^n_m$ для данного $n$ и всех $m$. Тогда:
$u^{n+1}_0=Ct^{1/k}$, где $t$ соответствует значению временнОго индекса $n+1$.
$u^{n+1}_M=0$
$u^{n+1}_m$ для $m=1..M-1$ находятся из разностного уравнения.

[svv]

Подумал, и вот к каким выводам пришел. Вот эта логика:

svv писал(а):

А так как в момент $n=0$ левая граница ещё холодная как лед ($Ct^{1/k}=0}$) и не греет, то и в момент $n=1$ вся область (кроме левой границы) ещё сохраняет исходную температуру: $u^1_m=0$ -- её ещё ничто не разогрело.

-- годилась бы для явной схемы.

Но у Вас неявная схема, и дело не в этом. Виноваты коэффициенты $a$. Они равны нулю для $n=0$ (всё дело в том, что вычисляются они для $n$, а не для $n+1$). Но тогда вся правая часть равна $0$. И Ваше разностное уравнение для $n=0$ можно записать просто как "левая часть равна нулю". Тогда, конечно, и $u^1_m=0$.

[dasalam]

svv в сообщении #496989 писал(а):

Подумал, и вот к каким выводам пришел. Вот эта логика:

svv писал(а):

А так как в момент $n=0$ левая граница ещё холодная как лед ($Ct^{1/k}=0}$) и не греет, то и в момент $n=1$ вся область (кроме левой границы) ещё сохраняет исходную температуру: $u^1_m=0$ -- её ещё ничто не разогрело.

-- годилась бы для явной схемы.

Но у Вас неявная схема, и дело не в этом. Виноваты коэффициенты $a$. Они равны нулю для $n=0$ (всё дело в том, что вычисляются они для $n$, а не для $n+1$). Но тогда вся правая часть равна $0$. И Ваше разностное уравнение для $n=0$ можно записать просто как "левая часть равна нулю". Тогда, конечно, и $u^1_m=0$.

Спасибо, но тогда будет и $u^2_m=0$ и $u^3_m=0$ и т.д. Но точное решение не нулевое

[svv]

dasalam писал(а):

Спасибо, но тогда будет и $u^2_m=0$

Давайте проверим.
Для определенности добавим еще один индекс в коэффициенты $a$:
$a^n_{m+1} = \frac{(u^{n}_{m})^k+(u^{n}_{m+1})^k}2$, что эквивалентно $a^n_m = \frac{(u^n_{m-1})^k+(u^n_m)^k}2$.
Ваше разностное уравнение я переписал так:
$\frac{u^{n+1}_m-u^n_m}\tau=\frac 1 h \left[a^n_{m+1} \frac{u^{n+1}_{m+1}-u^{n+1}_m} h-a^n_{m} \frac{u^{n+1}_m-u^{n+1}_{m-1}} h\right]$
Здесь еще два изменения в индексах коэффициентов $a$. Cогласитесь ли Вы с обоими?
1. Те индексы, что были -- пространственные, там должна быть буква $m$, а не $n$.
2. Не $a_{m+1}$ и $a_{m-1}$, а $a_{m+1}$ и $a_{m}$, чтобы у них было единое определение, приведенное выше. Оба слагаемых в правой части должны строиться по одному принципу, они отличаются только сдвигом на $1$ по $m$.

Дальше, в любой схеме -- явной или неявной -- должно быть $u^n_0=Ct^{1/k}=C(n\tau)^{1/k}$. Вы закладываете это в алгоритм? Если нет, тогда нуль будет тождественный, везде и всюду, "ненулю" просто неоткуда взяться.

Тезис. При $n=2$ величина $u^n_m$ отлична от нуля хотя бы при одном $m>0$.
Доказательство. Рассмотрим уравнение при $n=1$, $m=1$:
$\frac{u^2_1-u^1_1}\tau=\frac 1 h \left[a^1_2} \frac{u^2_2-u^2_1} h-a^1_1 \frac{u^2_1-u^2_0} h\right]$
Для $n=1$ только $u^1_0$ отлично от нуля. Допустим, для $n=2$ также только $u^2_0$ отлично от нуля. Тогда
$0=\frac 1 {h^2} a^1_1 u^2_0 \right]$
Однако $u^2_0$ и $a^1_1=\frac{(u^1_0)^k+(u^1_1)^k}2=\frac{(u^1_0)^k} 2$ оба отличны от нуля.
Значит, в случае, когда равны нулю все $u^2_m$, кроме $u^2_0$, нарушается разностное уравнение.

[dasalam]

2svv, очередное спасибо. Уже сам понял, что ошибался, вы подтвердили мои предположения.

[dasalam]

2svv, очередное спасибо. Уже сам понял, что ошибался, вы подтвердили мои предположения.
А может еще подскажите, вот есть такая разностная схема
$(1+\xi)*\frac{(U^{n+1}_{m}-U^{n}_{m})}t-\xi*\frac{(U^{n}_{m}-U^{n-1}_{m})}t = \frac{a_{m+1/2}(U^{n+1}_{m+1}-U^{n+1}_m)-a_{m-1/2}(U^{n+1}_{m}-U^{n+1}_{m-1})}{h^2}$
Здесь $\xi = const$. Как лучше такую систему решать?

[svv]

Я не специалист в этих вопросах -- отвечаю из общих соображений. Мне кажется, что примерно так же, просто теперь набор значений $U^{n+1}$ (для всех $m$ от $1$ до $M-1$) определяется не только набором значений $U^n$, но и набором значений $U^{n-1}$. Как рекуррентная формула $z_{n+1}=f(z_n, z_{n-1})$ -- не только по предыдущему, но и по пред-предыдущему элементу (только здесь "элемент" -- это значения $U$ для данного $n$ и всех $m$).

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

Скажите, а вот эта добавка
$\xi \left[ \frac{U^{n+1}_{m}-U^{n}_{m}}\tau-\frac{U^{n}_{m}-U^{n-1}_{m}}\tau \right]$,
соответствующая $\frac {\partial^2 U}{\partial t^2}$, она из усложненной физической модели происходит, или это нечто для улучшения сходимости?

[dasalam]

А как мы узнаем значения $U^n$? Если например для$ n=1$ из НУ и ГУ получаем только, что $U^0_{m} = 0$. Если я не ошибаюсь, то при $\xi = 1/2$ порядок апроксимации у этой схемы $O(t^2, h^2)$

[svv]

Я вот как раз и спрашивал по вторую производную по времени: если она есть в дифуравнении в частных производных, то начальные условия должны содержать помимо $U(x, t)\vline_{t=0}$ еще $\frac{\partial U(x, t)}{\partial t}\vline_{t=0}$.

Так эта добавка в левой части не "физическая", она для того, чтобы повысить порядок аппроксимации по $t$ с первого до второго?

[dasalam]

Интересно, а можно ли посчитать $U^1$ по схеме, котораю была предложена ранее, а потом подставить в "новую" схему. Я точно не знаю смысл этой добавки, но да она повышает порядок

[svv]

Да, я как раз Вам хотел это посоветовать. $U^1$ считаем по старому варианту (и получаем почти все нули), а начиная с $U^2$ -- уже по новому.

Все, понял. Не вторая производная, а просто линейная комбинация правой и левой конечной разности. Чтобы стала центральной разностью, надо взять $\xi=-1/2$, тогда
$(1+\xi) \frac{U^{n+1}_{m}-U^{n}_{m}}\tau-\xi \frac{U^{n}_{m}-U^{n-1}_{m}}\tau =\frac 1 2 \frac{U^{n+1}_{m}-U^{n}_{m}}\tau + \frac 1 2 \frac{U^{n}_{m}-U^{n-1}_{m}}\tau = \frac{U^{n+1}_{m}-U^{n-1}_{m}}{2\tau}$,
а центральная разность имеет второй порядок $O(\tau^2)$, в отличие от правой и левой $O(\tau)$.

Странный этот случай, $\xi=-1/2$. В уравнение не входят $U^n$, только $U^{n-1}$ и $U^{n+1}$. Выходит, $U$ с четными $n$ и $U$ с нечетными $n$ никак не "сцепляются" друг с другом. Я бы поостерегся использовать такую схему.