Теплопроводность. Обобщённая неявная разностная схема

[boh]

Сама задача описана в этой теме. Здесь прошу помочь с её численным решением. По заданию нужно пользоваться обобщённой неявной схемой.

Цитата:

Рассмотрим первую краевую задачу для уравнения теплопроводности с
постоянными коэффициентами

$$ \begin{cases} \frac{\partial u(x,t)}{\partial t}=a^{2}\frac{\partial^{2}u(x,t)}{\partial x^{2}}+F(x,t), & 0\leqslant t\leqslant T,\;0\leqslant x\leqslant l;\\ u(x,0)=\psi(x) & 0\leqslant x\leqslant l;\\ \frac{\partial u(0,t)}{\partial x}=\frac{\partial u(l,t)}{\partial x}=0 & 0<t\leqslant T.\end{cases}$$

Определим равномерную сетку, как множество узлов $(x_{i},t_{k})$,
где

$$ x_{i}=ih_{x},\; i=\overline{0,I},\; h_{x}=\frac{l}{I}$$

$$ t_{k}=kh_{t},\; k=\overline{0,K},\; h_{t}=\frac{T}{K}$$

Запишем приближенные выражения для производных, входящих в задачу:

$$ \frac{\partial^{2}u(x_{i},t_{k})}{\partial x^{2}}\thickapprox\frac{u(x_{i+1},t_{k})-2u(x_{i},t_{k})+u(x_{i-1},t_{k})}{h_{x}^{2}},$$

$$ \frac{\partial u(x_{i},t_{k})}{\partial t}\thickapprox\frac{u(x_{i},t_{k+1})-u(x_{i},t_{k})}{h_{t}}$$

$$ \frac{\partial u(x_{i},t_{k})}{\partial t}\thickapprox\frac{u(x_{i},t_{k})-u(x_{i},t_{k-1})}{h_{t}}$$

$$ \frac{\partial u(x_{i},t_{k})}{\partial x}\thickapprox\frac{u(x_{i},t_{k})-u(x_{i-1},t_{k})}{h_{t}}$$

Проводя дискретизацию задачи (...), заменим все функции непрерывных
аргументов $x$, $t$ их сеточными аналогами. Рассмотрим наряду с
основными узлами шаблона вспомогательный узел $(x_{i},t_{k}+\theta h_{t})$,
где $0\leqslant\theta<1$. Используя этот узел, запишем уравнения
простейших явной и неявной схем для задачи (...):

$$ \frac{\tilde{u}_{i}-u_{i}^{k}}{\theta h_{t}}=a^{2}\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{h_{x}^{2}}+F_{i}^{k};$$

$$ \frac{u_{i}^{k}-\tilde{u}_{i}}{\theta h_{t}}=a^{2}\frac{u_{i+1}^{k+1}-2u_{i}^{k+1}+u_{i-1}^{k+1}}{h_{x}^{2}}+F_{i}^{k+1},$$

где $\widetilde{u}_{i}$ - значение решения во вспомогательном узле.
Если умножить каждое из уравнений на $\frac{1}{\theta}$ и сложить,
а также дописать начальные и граничные условия, то получим обобщённую
неявную схему:

$$ \begin{cases} \frac{u_{i}^{k+1}-u_{i}^{k}}{h_{t}}=(1-\theta)(\Lambda_{x}^{h}u_{i}^{k+1}+F_{i}^{k+1})+\theta(\Lambda_{x}^{h}u_{i}^{k}+F_{i}^{k}), & i=\overline{1,I-1},\; k=\overline{0,K-1}\\ u_{i}^{0}=\psi & i=\overline{0,I}\\ \frac{u_{1}^{k}-u_{0}^{k}}{h_{x}}=0 & k=\overline{1,K}\\ \frac{u_{I}^{k}-u_{I-1}^{k}}{h_{x}}=0 & k=\overline{1,K}\end{cases}$$

где $\Lambda_{x}^{h}u_{i}^{k}=a^{2}\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{h_{x}^{2}}$.


Как это запрограммировать

Первое уравнение системы представим в следующем виде: в левой части
все члены $u^{k+1}$, остальные слагаемые - в правой части. Таким
образом, получаем систему вида

$$ A_{i}u_{i+1}^{k+1}+B_{i}u_{i}^{k+1}+C_{i}u_{i-1}^{k+1}=D_{i}^{k}$$

Система из $i$ уравнений с $i$ неизвестными решается методом прогонки.

Необходимо учесть, что заданы краевые условия $u_{0}^{k}=u_{1}^{k}$
и $u_{I-1}^{k}=u_{I}^{k}$. Уравнение примет следующий вид при $i=1$:

$$ A_{i}u_{2}^{k+1}+B_{i}u_{1}^{k+1}+C_{i}u_{0}^{k+1}=D_{i}^{k}\Longleftrightarrow A_{i}u_{2}^{k+1}+(B_{i}+C_{i})u_{1}^{k+1}=D_{i}^{k}$$

при $i=I-1$:

$$ A_{i}u_{I}^{k+1}+B_{i}u_{I-1}^{k+1}+C_{i}u_{I-2}^{k+1}=D_{i}^{k}\Longleftrightarrow C_{i}u_{I-2}^{k+1}+(A_{i}+B_{i})u_{I-1}^{k+1}=D_{i}^{k}.$$

Таким образом, трёхдиагональная матрица $M$ для метода прогонки будет
иметь элементы $B_{i}$ на главной диагонали, $A_{i}$ - на диагонали
над главной и $C_{i}$ - на диагонали под главной. Кроме того, с учётом
вышеприведённых краевых условий, исключениями будут

$$ M_{1\,1}=B_{1}+C_{1}$$

$$ M_{I-1\, I-1}=A_{I-1}+B_{I-1}$$

Система будет содержать $I-2$ уравнение с $I-2$ переменной. Ещё
2 значения получатся из краевых условий.

Итак, преобразования первого уравнения системы:

$$ \frac{u_{i}^{k+1}-u_{i}^{k}}{h_{t}}=(1-\theta)(a^{2}\frac{u_{i+1}^{k+1}-2u_{i}^{k+1}+u_{i-1}^{k+1}}{h_{x}^{2}}+F_{i}^{k+1})+\theta(a^{2}\frac{u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}}{h_{x}^{2}}+F_{i}^{k})$$

$$ \frac{1}{h_{t}}u_{i}^{k+1}-\frac{1}{h_{t}}u_{i}^{k}=\frac{a^{2}(1-\theta)}{h_{x}^{2}}\left(u_{i+1}^{k+1}-2u_{i}^{k+1}+u_{i-1}^{k+1}\right)+(1-\theta)F_{i}^{k+1}+\frac{a^{2}\theta}{h_{x}^{2}}\left(u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}\right)+\theta F_{i}^{k}$$

$$ u_{i+1}^{k+1}\left[\frac{a^{2}(1-\theta)}{h_{x}^{2}}\right]-u_{i}^{k+1}\left[\frac{1}{h_{t}}+2\frac{a^{2}\theta}{h_{x}^{2}}\right]+u_{i-1}^{k+1}\left[\frac{a^{2}(1-\theta)}{h_{x}^{2}}\right]=$$

$$ =-\frac{1}{h_{t}}u_{i}^{k}-(1-\theta)F_{i}^{k+1}-\frac{a^{2}\theta}{h_{x}^{2}}\left(u_{i+1}^{k}-2u_{i}^{k}+u_{i-1}^{k}\right)-\theta F_{i}^{k};$$

$$ M_{1\,1}=M_{I-1\, I-1}=-\frac{1}{h_{t}}-2\frac{a^{2}\theta}{h_{x}^{2}}+\frac{a^{2}(1-\theta)}{h_{x}^{2}}=-\frac{1}{h_{t}}-\frac{a^{2}\theta}{h_{x}^{2}}$$

Проблема в том, что этот метод даёт какие-то совершенно неправдоподобные результаты. В чём может быть проблема?

[ol_mer]

Ну, много чего может быть. Самое простое и легкопроверяемое - неправильное решение СЛАУ прогонкой.
Вы напишите, что и как вы проверяли.(ну а если не проверяли, так..редкие программы с первого раза работают)

[boh]

Мой код можно посмотреть здесь. Язык - Java, по ссылке проект для NetBeans.
В функции debugSweepButtonActionPerformed находится тестовая матрица, которая решается методом прогонки (функция gidsAbstractSweepSolve). Решение самой задачи (с описанными изменениями в матрице) происходит в функции gidsSolve.
Примеры работы программы:

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

[boh]

Ещё вопрос: можно ли применять эту схему к исходной задаче?

$$ \begin{cases} \bar{u}_{t}=a^{2}\bar{u}_{xx}+b(\bar{u}-u_{c})+f(x)\\ \bar{u}(x,0)=\psi(x) & 0\leqslant x\leqslant l\\ \bar{u}_{x}(0,t)=0 & 0<t\leqslant T\\ \bar{u}_{x}(l,t)=0\end{cases}$$