Промерзание грунта

[RAEman]

Здравствуйте ув.Форумчане.
Помогите разобраться в задаче Стефана:
$ c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2} , 0<x<S(t)\\ c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2} , S(t)<x<l\\ \lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}\\ \\ l = 50 \text{ метров}\\ c_1 \rho_1 = 2 \cdot10^6, \lambda_1 = 2\\ c_2 \rho_2 = 2 \cdot 10^6, \lambda_2 = 0,5\\ \\ x = S(t)\\ T_- = T_{\text{Ф}}=T_+\\ T(x,0) = T^0(x), t=0, 0 \leq x \leq l\\ T(0,t) = T_1(t) = \sin(\pi \cdot a \cdot x)\\ \begin{cases} T(l,t) = T_H &\\ &\text{ $t>0,x=l$;}\\ \frac{\partial T}{\partial x} = q & \end{cases}
Есть набросок программы которую я не совсем понимаю

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

#include "stdafx.h"
#include <math.h>
 
#define Pi 3.14159265358979
#define q_w 3.34*1000000
#define ro_w 1000
#define n 10
#define N 10
#define lambda2 0.58
#define lambda1 0.5
#define delta 0.01
#define h 0.1
#define c_s 700
#define ro_s 4000
#define c_w 4190
#define ro_w 1000
#define c_ice 2000
#define ro_ice 900
#define lambda_s 900
#define lambda_w 900
#define lambda_ice 900
#define m 0.4//пористость
#define Tfaza 0
#define Tnach -3
 
int c_ro(double S_w /*Водонасыщенность*/, double M /*Пористость*/){
    double temp;
    temp = (1-M)*c_s*ro_s+c_w*ro_w*S_w*M+c_ice*ro_ice*(1-S_w)*M;
    return temp;
}
int _lambda(double S_w /*Водонасыщенность*/, double M /*Пористость*/){
    double temp;
    temp = (1 - M)*lambda_s + lambda_w*S_w*M + lambda_ice*(1 - S_w)*M;
    return temp;
}
int Tgran(double a,double x){
    return sin(Pi*a*x);
}
int main(){
    int k;
    double tau = 1, t = 1, maxy = 1, alpha[n+1], beta[n+1], lambda, cro, y[n+1], yj[n+1],a , b, c, d;
    for (int j = 0; j <= N; j++){ 
        y[j] = Tnach; 
        yj[j] = y[j]; 
    }
    for (int j = 1; j < N; j++){
        tau = tau*1.1;
        if (tau>1e6) tau = 1e6;
        t = t + tau;
        k = 0;
        while ((maxy > 1e-4) && (k < 20)){
            k = k + 1;
            alpha[0] = 0; beta[0] = Tgran(1,1)/*Синус пока от балды написал*/;
            for (int i = 1; i < n; i++ ){
                if (y[i] < Tfaza) lambda = _lambda(1.0,m); 
                else lambda = _lambda(0.0, m);
                if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m); 
                else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m); 
                else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);
                a = lambda / (h * h);
                b = a;
                c = cro / tau;
                d = c*yj[i];
                c = a + b + c;
                alpha[i] = b / (c - a*alpha[i - 1]);
                beta[i] = (d + a*beta[i - 1]) / (c - a*alpha[i - 1]);
            }
            y[n] = Tnach;
            y[n] = beta[n - 1] / (1 - alpha[n - 1]);
            maxy = y[n];
            for (int i = n - 1; i >= 0; i--){
                a = y[i];
                y[i] = y[i + 1] * alpha[i] + beta[i];
                if (maxy < abs(y[i] - a)) maxy = abs(y[i] - a);
            }
            for (int i = 1; i <= n; i++){
                printf("%f", y[i]);
            }
            printf("\n");
            yj[j] = y[j];
            
        }
        getchar();
    }
}
 

Я понимаю, что из постановку задачи мы аппроксимируем, получаем сетку, (и потом вроде как сглаживание(? не знаю что это) проводим) и ее прогонкой решаем.
1)А что эта сетка означает и как она мне помогает построить график(скажем в той же проге) я не понимаю.
2)Зачем было введено maxy.
3)И как мне по времени разбить результат

[svv]

1. Понятен ли Вам физический смысл всех уравнений, которые Вы записали?
2. Понятно ли, зачем в численных методах решения ДУЧП вообще вводятся сетки?
3. Смогли бы Вы написать расчётные формулы, по которым должна работать программа (не та, что Вы привели, а Ваша, ещё не написанная)?

[RAEman]

svv

Цитата:

1. Понятен ли Вам физический смысл всех уравнений, которые Вы записали?

Нет, что я и хочу понять....

Цитата:

2. Понятно ли, зачем в численных методах решения ДУЧП вообще вводятся сетки?

Не очень, видимо для численной реализации...

Цитата:

3. Смогли бы Вы написать расчётные формулы, по которым должна работать программа (не та, что Вы привели, а Ваша, ещё не написанная)?

Аппроксимацию, по которой и строится сетка, я вроде бы усвоил и подогнать для программы могу

[svv]

Вообще я в этой теме не специалист, но смысл уравнений, вроде, понятен.
Здесь, похоже, используется такая модель. Есть две среды, «одномерные», с одной координатой $x$ — глубиной (от других координат никакие величины не зависят и потому они в задаче не фигурируют). Среда 1 — промёрзший грунт, среда 2 — непромёрзший. Вся «фишка» задачи в том, что глубина границы раздела $S(t)$, то есть глубина промерзания, зависит от времени. Грубо говоря, если становится теплее, граница поднимается, а если холоднее — опускается. Каждая среда характеризуется набором известных постоянных параметров:
$c$ — удельная теплоёмкость (массовая), $\frac{\text{Дж}}{\text{кг}\cdot\text{K}}$;
$\rho$ — плотность, $\frac{\text{кг}}{\text{м}^3}$;
$\lambda$ — удельная теплопроводность, $\frac{\text{Вт}}{\text{м}\cdot\text{K}}$, она же коэффициент теплопроводности.
Кроме того, для описания перехода воды из твердого в жидкое состояние и обратно нужна ещё плотность воды $\rho_w$ и её удельная теплота плавления $q_w$, $\frac{\text{Дж}}{\text{кг}}$. Похоже, в программе ошибка, потому что по таблицам последняя величина равна $3.3477\cdot 10^5 \frac{\text{Дж}}{\text{кг}}$, а в программе q_w=3.34*1000000, то есть в 10 раз больше.
В задаче рассматриваются глубины не до ядра Земли , а только до некоторой глубины $\ell=50\text{м}$. Вероятно, ниже этой глубины температура практически постоянна, независимо от времени года.

Помимо констант, в задаче имеется неизвестная функция — температура $T(x, t)$, зависящая не только от координаты $x$, но и от времени $t$. Собственно, в задаче и требуется найти зависимость $T(x, t)$. Как известно, изменение распределения температуры в среде со временем подчиняется уравнению теплопроводности. В Ваших обозначениях оно имеет вид $c\rho\frac{\partial T}{\partial t} = \lambda \frac{\partial^2 T}{\partial x^2}$.
Уравнение, приведённое в Википедии, отличается от Вашего тем, что там 1) константы собраны в одну, 2) имеются дополнительные внутренние источники тепла и 3) рассматривается трехмерный случай.

$c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2}, \quad 0<x<S(t)$
Это уравнение теплопроводности для среды 1 (промёрзший грунт). Среда 1 имеет координаты от $0$ (поверхность земли) до $S(t)$ (глубина промерзания в момент $t$).

$c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2} , \quad S(t)<x<\ell$
Это уравнение теплопроводности для среды 2 (непромёрзший грунт). Среда 2 имеет координаты от $S(t)$ до $\ell$ (дальше не рассматриваем).

$\lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}$
Это уравнение описывает смещение границы раздела сред со временем, то есть процесс промерзания и оттаивания грунта. В левой части стоит разность потоков тепла в обеих средах у границы раздела. Эта разница идёт на замораживание или размораживание грунта на границе раздела, что и приводит к смещению границы. По сути, это уравнение баланса энергии на границе. Тепловая энергия, необходимая для перевода 1 килограмма льда в жидкое состояние, известна — это $q_w$, а $\rho_w$ позволяет пересчитать массу в объём.

$\ell = 50 \text{м}$
Выбрали такую предельную глубину, до которой рассматриваются процессы теплопередачи.

$c_1 \rho_1 = 2 \cdot10^6, \lambda_1 = 2$
$c_2 \rho_2 = 2 \cdot 10^6, \lambda_2 = 0,5$
Для каждой среды можно было вообще все константы перенести в правую часть и заменить их одним коэффициентом температуропроводности: $\alpha=\frac{\lambda}{c\rho}$ (как сделано в Википедии), тогда бы уравнения теплопроводности приняли вид $\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}$, где $\alpha$ в каждой из сред своё.

$T_- = T_{\text{Ф}}=T_+, \quad\quad x=S(t)$
На границе раздела сред $x = S(t)$ температура сверху от границы $T_-$ равна температуре снизу от границы $T_+$, и обе совпадают с температурой фазового перехода лёд-вода $T_{\text{Ф}}$.

В задаче заданы начальные и граничные условия.

$T(x,0) = T^0(x), \quad\quad 0\leqslant x \leqslant l,\;\;t=0$
Начальная температура в момент времени $t=0$ в каждой точке из рассматриваемого диапазона координат равна заданной функции $T^0(x)$.

$T(0,t) = T_1(t) = \sin(\pi a t),\quad\quad x=0,\;\;t>0$
Температура в точке $x=0$, т.е. на поверхности земли, в любой момент времени $t$ равна заданной функции $T_1(t)=\sin(\pi at)$. Эта функция, вероятно, описывает годовые колебания температуры: один период синусоиды — один год. У Вас в функции было $\pi a x$, я исправил на $\pi a t$, здесь нужна зависимость от времени.

$T(\ell,t) = T_H, \;\;\frac{\partial T}{\partial x}(\ell,t) = q, \quad\quad x=\ell,\;\;t>0$
Температура и её производная по координате на нижней границе среды 2 в любой момент времени постоянны и заданы (на такой глубине погода уже не сказывается ).

[RAEman]

svv
Спасибо за столь подробное разьяснение постановки
Вы можете мне сказать зачем вообще было введено ограничение (maxy>1e-4) и что оно означает?

-- 20.05.2015, 16:50 --

Я забыл указать в постановке еще одно уравнение
$(\tilde{c \rho} + \delta(T_\text{Ф} - T) \cdot q_w \cdot \rho_w)\frac{\partial T}{\partial t} = \tilde{\lambda}\frac{\partial T}{\partial x}$
И мне кажется видимо первые два (или три уравнения) судя по программе

Цитата:

$c_1 \rho_1\frac{\partial T}{\partial t} = \lambda_1 \frac{\partial^2 T}{\partial x^2} , 0<x<S(t)\\ c_2 \rho_2\frac{\partial T}{\partial t} = \lambda_2 \frac{\partial^2 T}{\partial x^2} , S(t)<x<l\\ \lambda_1 \frac{\partial T}{\partial x} - \lambda_2 \frac{\partial T}{\partial x} = q_w \rho_w \frac{\partial S}{\partial t}\\$

были заменены им для выставления условий:

if (y[i] < Tfaza) lambda = _lambda(1.0,m);
else lambda = _lambda(0.0, m);
if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m);
else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m);
else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);

-- 20.05.2015, 16:52 --

либо я что-то путаю

[svv]

Программу я ещё не смотрел, чуть позже посмотрю и постараюсь ответить.

По-хорошему, между дифференциальными уравнениями и программой должен быть ещё один этап — расчётные формулы, или запись алгоритма «почти как в будущей программе», но в математических обозначениях. Чтобы было понятно — вот пример. (Даже всю ту тему бегло просмотрите.)

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

[RAEman]

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

По-хорошему, между дифференциальными уравнениями и программой должен быть ещё один этап — расчётные формулы

Эти формулы у меня уже есть, программа на них вроде нормально работает, но с какими-то ошибками(вероятнее всего в вычислении maxy)
Используется неявная схема(она мне больше нравится ибо там не надо проверять на устойчивость)
$c \rho \frac{y_i}{\tilde{y_i}}{\tau} = \lambda \frac{y_{i+1}+2y_{i}+y_{i-1}}{h^2}\\ y_i(\frac{\tilde{c \rho}}{\tau} + \frac{2 \tilde{\lambda}}{h^2}) = (\frac{c \rho \tilde{y_i}}{\tau} + y_{i+1}\frac{\lambda}{h^2}+y_i-1\frac{\lambda}{h^2})
Используется метод прогонки, т.е. $a_i {y__{i-1}} - c_i y_i + b_i y_{i+1} +d_i$
Вы про этот этап имели ввиду?

[svv]

Да, про этот.
Пожалуйста, подождите день-два, пока я посмотрю Вашу программу.

[RAEman]

svv
Хорошо, Спасибо.

[svv]

Предварительно: я нашел в программе несколько ошибок, из-за которых она вряд ли будет работать. Вы можете аргументированно возражать . Для удобства я их пронумерую.

1. Согласно таблицам, для воды (льда) удельная теплота плавления воды равна $3.35\cdot 10^5 \frac{\text{Дж}}{\text{кг}}$, а в программе в 10 раз больше.

2. Дальше обратите внимание на формулы для $c\rho$ и $\lambda$. Пусть индекс $s$ означает грунт сам по себе, $w$ — вода в жидком состоянии, $i$ — лёд. Пусть $M$ — пористость (т.е. доля воды в любом состоянии в общем количестве вещества), $S_w$ — водонасыщенность (доля воды в жидком состоянии в общем количестве воды). Тогда, действительно, считая, что каждая фракция вносит вклад пропорционально своей доле, получаем формулы вроде такой:
$\lambda=(1 - M)\lambda_s +M S_w \lambda_w + M (1-S_w) \lambda_i$
Именно так в программе. Здесь всё правильно и логично. Первое слагаемое — доля грунта в веществе, т.е. $1-M$. Грунт не делится на твердый и жидкий. Доля воды $M$. Доля жидкой воды в веществе $MS_w$ (второе слагаемое), а доля твердой воды $M(1-S_w)$ (третье слагаемое).
Чтобы лучше почувствовать написанную формулу для $\lambda$, посмотрите, что она даёт в трёх случаях:
$\bullet$ $M=0,\; S_w$ любое
$\bullet$ $M=1,\; S_w=1$
$\bullet$ $M=1,\; S_w=0$

Понятно, что $M$ от температуры не зависит: если в грунте есть вода, изменение температуры может только перевести её в другое агрегатное состояние, но не может изменить долю воды. А вот $S_w$ от температуры зависит: когда очень тепло, вся вода жидкая ($S_w=1$), когда очень холодно, вся вода твёрдая ($S_w=0$). Так должно быть.

Но теперь посмотрите, как в программе вычисляется «суммарная» $\lambda$ грунта:
if (y[i] < Tfaza) lambda = _lambda(1.0,m);
else lambda = _lambda(0.0, m);
и аналогично вычисляется «суммарное» $c\rho$. Всё наоборот. Когда температура ниже фазового перехода, используется водонасыщенность $1$, иначе $0$.

3. Функции _lambda и c_ro вычисляют число типа double (с плавающей точкой, двойной точности), это правильно. Но поскольку тип возвращаемого значения у них обеих int, заботливо вычисленная дробная часть отбрасывается и остаётся только целая часть.

4. Посмотрите, как вычисляется $c\rho$:
if (y[i] < (Tfaza - delta)) cro = c_ro(1.0,m);
else if (y[i] < (Tfaza + delta)) cro = c_ro(0.0,m);
else cro = c_ro(1.0,m) + m*q_w*ro_w / (2 * delta);
Первая альтернатива должна сработать, когда температура определённо ниже фазового перехода ($T<T_{\text{Ф}}-\delta$). Про неправильную водонасыщенность я уже сказал в п.2, но знак < правильный.
Вторая альтернатива должна сработать, когда температура определённо выше фазового перехода ($T>T_{\text{Ф}}+\delta$). Во второй строке кода знак < неправильный.
Третья альтернатива должна сработать в остальных случаях, когда температура с точностью до небольшого зазора $\pm\delta$ равна температуре фазового перехода $T_{\text{Ф}}$. Здесь теплоёмкость записана с дополнительным слагаемым, которое представляет собой удельную теплоту плавления, «размазанную» по зазору $T_{\text{Ф}}-\delta<T<T_{\text{Ф}}+\delta$. Это я могу понять.

Исправьте все эти ошибки и покажите, что получилось. Скоро дойдём и до maxy.

[RAEman]

svv
Я не знаю чему изначально будет равна $\tau$ и $t$ и поэтому приравнял от балды, как и $T_{gran}$
Вот исправленная программа

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

#include <math.h>
#include <stdio.h>
#include <stdlib.h> 
#define Pi 3.14159265358979
#define q_w 3.34*100000
#define ro_w 1000
#define n 100
#define N 100
#define lambda2 0.58
#define lambda1 0.5
#define delta 0.01
#define h 0.1
#define c_s 700
#define ro_s 4000
#define c_w 4190
#define ro_w 1000
#define c_ice 2000
#define ro_ice 900
#define lambda_s 900
#define lambda_w 900
#define lambda_ice 900
#define m 0.4
#define Tfaza 0
#define Tnach -3
double c_ro(double S_w, double M)
{
        double temp;
        temp = (1-M)*c_s*ro_s+c_w*ro_w*S_w*M+c_ice*ro_ice*(1-S_w)*M;
        return temp;
}
double _lambda(double S_w , double M)
{
        double temp;
        temp = (1 - M)*lambda_s + lambda_w*S_w*M + lambda_ice*(1 - S_w)*M;
        return temp;
}
double Tgran(double a,double t)
{
        return sin(Pi*a*t);
}
int main()
{
        FILE *f = fopen("1.txt", "w");
        int k,i,j;
        double tau = 1, t = 0, maxy = 1, alpha[n+1], beta[n+1], lambda, cro, y[n+1], yj[n+1],a , b, c, d;
        for (j = 0; j <= N; j++)
        { 
                y[j] = Tnach;
                yj[j] = y[j]; 
        }
        for (j = 1; j < 100; j++)
        {
                tau = tau*1.1;
                if (tau>1e6) tau = 1e6;
                t = t + tau;
                k = 0;
                while ((maxy > 1e-4) && (k < 20))
                {
                        k = k + 1;
                        alpha[0] = 0; beta[0] = Tgran(1,t);
                        for (i = 1; i < n; i++)
                        {
                                if (y[i] < Tfaza) lambda = _lambda(0.0, m);
                                else lambda = _lambda(1.0, m);
                                if (y[i] < (Tfaza - delta)) cro = c_ro(0.0, m);
                                else if (y[i] > (Tfaza + delta)) cro = c_ro(1.0, m);
                                else cro = c_ro(0.0, m) + m*q_w*ro_w / (2 * delta);
                                a = lambda / (h * h);
                                b = a;
                                c = cro / tau;
                                d = c*yj[i];
                                c = a + b + c;
                                alpha[i] = b / (c - a*alpha[i - 1]);
                                beta[i] = (d + a*beta[i - 1]) / (c - a*alpha[i - 1]);
                        }
//                      y[n] = Tnach; //1 rod
                        y[n] = beta[n - 1] / (1 - alpha[n - 1]); //2 rod
                        maxy = y[n];
                        for (i = n - 1; i >= 0; i--)
                        {
                                a = y[i];
                                y[i] = y[i + 1] * alpha[i] + beta[i];
                                if (maxy < abs(y[i] - a)) maxy = abs(y[i] - a);
                        }
//                      printf("%d::",k);
                        for (int i = 1; i <= n; i++) printf("%.5f |",y[i]);
                        printf("\n\n");
                        yj[j] = y[j];
                }
        }
}

и результат практически тот же самый
https://yadi.sk/i/m6XXXzOOgpHj6

[RAEman]

svv
И еще одно замечание по вашим исправлениям:

Цитата:

Понятно, что $M$ от температуры не зависит: если в грунте есть вода, изменение температуры может только перевести её в другое агрегатное состояние, но не может изменить долю воды.

Пористость у меня по условию равна 0,4.

[svv]

В своих попытках понять логику работы программы (ох, легче свою с нуля написать!) я привёл программу к такому виду:

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

const int
   n=100,
   N=100;
const double
   Pi=M_PI,
   q_w=3.34e5,
   delta=0.01,
   h=0.1,
   c_s=700,
   ro_s=4000,
   c_w=4190,
   ro_w=1000,
   c_ice=2000,
   ro_ice=900,
   lambda_s=900,
   lambda_w=900,
   lambda_ice=900,
   m=0.4,
   Tfaza=0,
   Tnach=-3;

double get_lambda(double T)
{
   double S_w = T<Tfaza? 0: 1;
   return (1-m)*lambda_s + lambda_w*S_w*m + lambda_ice*(1-S_w)*m;
}

double get_cro(double T)
{
   double
      S_w = (T <= Tfaza+delta)? 0: 1,
      c_ro = (1-m)*c_s*ro_s + c_w*ro_w*S_w*m + c_ice*ro_ice*(1-S_w)*m;
   if (fabs(T - Tfaza)<=delta) c_ro+=m*q_w*ro_w / (2 * delta);
   return c_ro;
}

double Tgran(double a, double t)
{
   return sin(Pi*a*t);
}

double max(double a, double b)
{
   return a>b? a: b;
}

double min(double a, double b)
{
   return a<b? a: b;
}

void Calc()
{
   double t = 0, maxy = 1, y[n+1], yj[n+1];
   for (int j = 0; j <= N; j++)
   {
      y [j] = Tnach;
      yj[j] = Tnach;
   }
   for (int j = 1; j < 100; j++)
   {
      double tau = min(pow(1.1, j), 1e6);
      t = t + tau;
      for (int k = 1; k <= 20 && maxy > 1e-4; k++)
      {
         double alpha[n+1], beta[n+1];

         alpha[0] = 0;
         beta [0] = Tgran(1,t);
         for (int i = 1; i < n; i++)
         {
            double
               lambda = get_lambda(y[i]),
               cro = get_cro(y[i]),
               a = lambda/(h*h),
               b = a,
               c = cro / tau,
               d = c*yj[i],
               e = a + b + c;
            alpha[i] = b / (e - a*alpha[i - 1]);
            beta[i] = (d + a*beta[i - 1]) / (e - a*alpha[i - 1]);
         }

         y[n] = beta[n - 1] / (1 - alpha[n - 1]);
         double maxy = y[n];
         for (int i = n - 1; i >= 0; i--)
         {
            double a = y[i];
            y[i] = y[i + 1] * alpha[i] + beta[i];
            maxy = max(maxy, abs(y[i] - a));
         }
         yj[j] = y[j];
      }
   }
}
 

Вам не нужно делать так, как у меня. Это только для пояснения хода мыслей. Здесь язык C++, а у Вас — C или совместимое с ним подмножество C++. Я постарался сделать область видимости переменных как можно более узкой, чтобы легче было анализировать их поведение.

В строке
if (maxy < abs(y[i] - a)) maxy = abs(y[i] - a);
вместо abs точно надо поставить fabs, иначе будет преобразование к целому.

Много непонятного. Самый внешний цикл, с управляющей переменной j, очевидно, по дискретному времени. Почему с каждым следующим j шаг по времени tau растёт, вплоть до астрономической величины $10^6$?
Группа операторов, вычисляющих a, b, c, d и затем alpha[i], beta[i], малопонятна. Не стоит одну и ту же переменную (c) использовать в двух разных смыслах, это источник ошибок, лучше завести дополнительную (у меня e). Всю эту группу надо тщательно проверять.
Совершенно непонятен оператор yj[j] = y[j];. Почему надо копировать в yj только одно значение, почему именно это, и почему это надо делать при каждом k? Можно ведь сделать это один раз, после цикла по $k$, когда элемент y[j] примет окончательное для данного j, «устаканившееся» значение.

Что касается maxy. Это максимальное значение модуля производной от y[i] (т.е. температуры). Вот только пространственной или временной производной — Вам виднее. Когда этот максимальный модуль производной становится меньше $10^{-4}$, итерации (т.е. цикл по k) можно завершить, не дожидаясь k=20 — предельного количества итераций.

Похоже, дальше сможете продвинуться только Вы сами.

[RAEman]

svv
Спасибо большое, за столь подробный и полный ответ. Мне многое наконец понятно Благодаря Вам.

svv писал(а):

Много непонятного.
Совершенно непонятен оператор yj[j] = y[j];. Почему надо копировать в yj только одно значение, почему именно это, и почему это надо делать при каждом k? Можно ведь сделать это один раз, после цикла по $k$, когда элемент y[j] примет окончательное для данного j, «устаканившееся» значение.

Видимо это надо для того чтобы из каждого вычисленного слоя сделать "исходный" слой, от которого можно вычислить последующий.
Насчет fabc спасибо теперь у меня есть хоть какие-то нормальные результаты. Однако при изменении константы n у меня выводит ошибку сегментации :(
И Спасибо за разъяснение про maxy, но меня волнует не столь его значение, сколько вопрос "для чего он тут? и для чего сделан цикл с maxy и k?".

[svv]

По всей видимости, в этой программе уравнение решается методом последовательных приближений. Процесс последовательных приближений должен быстро сходиться к точному решению уравнения. Для этого в программе имеется цикл по k — это номер последовательного приближения (итерации). Этот цикл не должен быть бесконечным. Из него существует два выхода:
1) выход «по-хорошему» — достигли необходимой точности, а maxy — индикатор этой точности: когда он становится достаточно малым (1e-4), процесс можно прекращать; почему именно maxy является индикатором — другой вопрос;
2) выход «по-плохому», когда совершили уже 20 итераций, а процесс почему-то не сходится; в такой ситуации ясно, что «пора кончать», несмотря на неудачу — и для этого проверяется условие k<=20.