Требуется написать программу, численно решающую неявное уравнение теплопроводности неявной схемой.
Методичка гласит, что
Цитата:
Неявную схему рассчитываем итерациями. Простые итерации сходятся медленно, запас устойчивости для их сходимости невелик, поэтому чаще всего используется итерационная процедура Ньютона.
Суть этой процедуры:
Присваиваем неизвестной величине

индекс итерации

.
Нелинейное относительно неизвестных

слагаемые в правой части
уравнения:
![$$h_j\frac{T^{n+1}_j-T^n_j}{\tau}=\sigma\left[\frac{T^{n+1}_{j+1}-T^{n+1}_j}{h_{j+\frac{1}{2}}}\kappa^{n+1}_{j+\frac{1}{2}}-\frac{T^{n+1}_j-T^{n+1}_{j-1}}{h}\kappa^{n+1}_{j-\frac{1}{2}}\right]$$ $$h_j\frac{T^{n+1}_j-T^n_j}{\tau}=\sigma\left[\frac{T^{n+1}_{j+1}-T^{n+1}_j}{h_{j+\frac{1}{2}}}\kappa^{n+1}_{j+\frac{1}{2}}-\frac{T^{n+1}_j-T^{n+1}_{j-1}}{h}\kappa^{n+1}_{j-\frac{1}{2}}\right]$$](https://dxdy-03.korotkov.co.uk/f/6/6/7/667a3b258ab829183c7a2a9aba54353482.png)
раскладываем в ряд по малому параметру:


Ограничимся линейными членами разложения. Имеем:
![$$\begin{gathered}h\frac{T^{n+1}_j+\Delta{T}-T^n_j}{\tau}=\sigma\left[\frac{T^{n+1,k}_{j+1}-T^{n+1,k}_j}{h}\kappa^{n+1,k}_{j+\frac{1}{2}}-\frac{T^{n+1,k}_j-T^{n+1,k}_{j-1}}{h}\kappa^{n+1,k}_{j-\frac{1}{2}}\right]+\\
+\frac{\sigma\Delta{T}_{j+1}}{h_{j+\frac{1}{2}}}\left[\kappa^{n+1,k}_{j+\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j+\frac{1}{2}}(T^{n+1,k}_{j+1}-T^{n+1,k}_j)\right]+\\
+\frac{\sigma\Delta{T}_{j}}{h_{j+\frac{1}{2}}}\left[-\kappa^{n+1,k}_{j+\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j+\frac{1}{2}}(T^{n+1,k}_{j+1}-T^{n+1,k}_j)\right]-\\
-\frac{\sigma\Delta{T}_{j}}{h_{j-\frac{1}{2}}}\left[\kappa^{n+1,k}_{j-\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j-\frac{1}{2}}(T^{n+1,k}_{j}-T^{n+1,k}_{j-1})\right]+\\
+\frac{\sigma\Delta{T}_{j-1}}{h_{j-\frac{1}{2}}}\left[-\kappa^{n+1,k}_{j-\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j-\frac{1}{2}}(T^{n+1,k}_{j}-T^{n+1,k}_{j-1})\right]\end{gathered}$$ $$\begin{gathered}h\frac{T^{n+1}_j+\Delta{T}-T^n_j}{\tau}=\sigma\left[\frac{T^{n+1,k}_{j+1}-T^{n+1,k}_j}{h}\kappa^{n+1,k}_{j+\frac{1}{2}}-\frac{T^{n+1,k}_j-T^{n+1,k}_{j-1}}{h}\kappa^{n+1,k}_{j-\frac{1}{2}}\right]+\\
+\frac{\sigma\Delta{T}_{j+1}}{h_{j+\frac{1}{2}}}\left[\kappa^{n+1,k}_{j+\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j+\frac{1}{2}}(T^{n+1,k}_{j+1}-T^{n+1,k}_j)\right]+\\
+\frac{\sigma\Delta{T}_{j}}{h_{j+\frac{1}{2}}}\left[-\kappa^{n+1,k}_{j+\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j+\frac{1}{2}}(T^{n+1,k}_{j+1}-T^{n+1,k}_j)\right]-\\
-\frac{\sigma\Delta{T}_{j}}{h_{j-\frac{1}{2}}}\left[\kappa^{n+1,k}_{j-\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j-\frac{1}{2}}(T^{n+1,k}_{j}-T^{n+1,k}_{j-1})\right]+\\
+\frac{\sigma\Delta{T}_{j-1}}{h_{j-\frac{1}{2}}}\left[-\kappa^{n+1,k}_{j-\frac{1}{2}}+(\frac{1}{2}\frac{\partial{f}}{\partial{T}})^{n+1,k}_{j-\frac{1}{2}}(T^{n+1,k}_{j}-T^{n+1,k}_{j-1})\right]\end{gathered}$$](https://dxdy-01.korotkov.co.uk/f/c/3/2/c32ecfc6f8b7c644dd4f0eae2f73f66382.png)
Система трехдиагональна относительно неизвестных

и вычисляется методом прогонки. В качестве

обычно используют

. Для получения
хорошей точности (4-6 знаков) достаточно 3-5 итераций.
Соотвественно, для этой задачи я переписал уравнение относительно

, составил формулу для вычисления прогоночных коэффициентов и написал программу, которая должна всё это вычислять:
void Plot::setImplicit()
{
double *A = new double [N], *B = new double [N],
*dT = new double [N+1];
for (int k=0; k<5; k++)
{
double *T;
if (k==0) T = Tb;
else T = Td;
A[0]=0;
B[0]=0;
for (int j=1; j<N; j++)
{
//qDebug()<<T[j];
double k1 = pow((T[j-1]+T[j])*0.5,25);
double k2 = pow((T[j]+T[j+1])*0.5,25);
double dk1 = 25*pow((T[j-1]+T[j])*0.5,24);
double dk2 = 25*pow((T[j]+T[j+1])*0.5,24);
//if (k==0) qDebug()<<k1<<" "<<k2;
A[j] = (-1/h*(k2+0.5*dk2*(T[j+1]-T[j])))/
(1/h*(-k1*(A[j-1]+1)+0.5*dk1*(A[j-1]-1)*(T[j]-T[j-1])-k2+0.5*dk2*(T[j+1]-T[j]))-h/tau);
B[j] = ((h/tau*(T[j]-Tb[j])-k2/h*(T[j+1]-T[j])+k1/h*(T[j]-T[j-1]))-1/h*(-k1+0.5*dk1*(T[j]-T[j-1]))*B[j-1])/
(1/h*(-k1*(A[j-1]+1)+0.5*dk1*(A[j-1]-1)*(T[j]-T[j-1])-k2+0.5*dk2*(T[j+1]-T[j]))-h/tau);
qDebug()<<A[j]<<" "<<B[j];
}
dT[N]=B[N-1]/(1-A[N-1]);
Td[N]=T[N]+dT[N];
for (int j=N-1; j>=0; j--)
{
dT[j]=A[j]*dT[j+1]+B[j];
Td[j]=T[j]+dT[j];
//qDebug()<<dT[j]<<" "<<T[j]<<" "<<Td[j];
}
}
for (int j=0; j<=N; j++)
{
//qDebug()<<Td[j];
Tb[j]=Td[j];
}
}
Результаты, мягко говоря, неадекватные. Где у меня ошибка?