Нелинейное ур-е теплопроводности, неявная схема - реализация

[byulent]

Требуется написать программу, численно решающую неявное уравнение теплопроводности неявной схемой.
Методичка гласит, что

Цитата:

Неявную схему рассчитываем итерациями. Простые итерации сходятся медленно, запас устойчивости для их сходимости невелик, поэтому чаще всего используется итерационная процедура Ньютона.
Суть этой процедуры:
Присваиваем неизвестной величине $T^{n+1}$ индекс итерации $k$.
Нелинейное относительно неизвестных $T^{n+1}$ слагаемые в правой части
уравнения:

$$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]$$

раскладываем в ряд по малому параметру:

$$\Delta T_j \equiv T^{n+1,k+1}_j-T^{n+1,k}_j,$$

$$f(T^{n+1,k+1})=f(T^{n+1,k})+(\frac{\partial{f}}{\partial{T}})_{n+1,k}\Delta{T}+...$$

Ограничимся линейными членами разложения. Имеем:

$$\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}$$

Система трехдиагональна относительно неизвестных $\Delta T_j$ и вычисляется методом прогонки. В качестве $T^{n+1, 0}$ обычно используют $T^n$. Для получения
хорошей точности (4-6 знаков) достаточно 3-5 итераций.

Соотвественно, для этой задачи я переписал уравнение относительно $\Delta T$, составил формулу для вычисления прогоночных коэффициентов и написал программу, которая должна всё это вычислять:

Используется синтаксис C++

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];
    }
}
 

Результаты, мягко говоря, неадекватные. Где у меня ошибка?

[Pphantom]

byulent в сообщении #1179327 писал(а):

double k1 = pow((T[j-1]+T[j])*0.5,25);

и последующее в том же духе...

Вы всерьез считаете, что $T^n_j$ - это $T_j$ в степени $n$? Откуда там такие бешеные степени?

[byulent]

Pphantom в сообщении #1179336 писал(а):

Вы всерьез считаете, что $T^n_j$ - это $T_j$ в степени $n$?

Нет, а с чего вы взяли?

Используется синтаксис C++

double *T;
        if (k==0) T = Tb;
        else T = Td; 

Это значит, что на 0-й итерации у меня берутся значения на предыдущем временном шаге ($T^{n+1,0}=T^n$), а на остальных итерациях - на недосчитанном текущем.

Pphantom в сообщении #1179336 писал(а):

Откуда там такие бешеные степени?

Из задания. Явная схема с этим вполне себе работала.

[Pphantom]

byulent в сообщении #1179338 писал(а):

Из задания. Явная схема с этим вполне себе работала.

А где они в задании?

[byulent]

Pphantom в сообщении #1179339 писал(а):

А где они в задании?

Ну, собственно, вот задание: $\kappa^n_j = \kappa_0 (T^n_j)^{25}, \kappa_0=1$.

[Pphantom]

Да уж, "вовремя" Вы признались... Ну что ж, тогда возражение выше действительно снимается, но задача становится довольно веселенькой. Обычно в таких случаях используют не просто прогонку, но ее модификацию - т.н. метод потоковой прогонки. У Вас то вроде бы не так. Почему?

[byulent]

Pphantom в сообщении #1179343 писал(а):

Обычно в таких случаях используют не просто прогонку, но ее модификацию - т.н. метод потоковой прогонки. У Вас то вроде бы не так. Почему?

"Я и словей-то таких не знаю!" (с).
Ну, если в методичке написано, что эта задача может решаться обычной прогонкой, значит, она может ей решаться.
Вы считаете, тут проблема именно в степени?

[Pphantom]

byulent в сообщении #1179397 писал(а):

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

Ну, совсем уж ручаться не буду, но скорее всего. А Вы можете написать полностью условие задачи? Все-таки восстанавливать его из промежуточных результатов несколько неудобно...

[Munin]

Ничего себе степень! Это ядерная реакция в центре звезды?