Численное решение уравнения теплопроводности

@Ivan789 · Регистрация: 02.04.2022

Студворк — интернет-сервис помощи студентам

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

Уравнение:
$https://www.cyberforum.ru/cgi-bin/latex.cgi?\small \begin{cases} & {u'}_{t} = {u''}_{xx}, 0\leq x\leq 1, t\geq 0\\ &\\ & u(0,x) = -2{x}^{3}-x+3, 0\leq x\leq 1\\ &\\ & {u'}_{x}(0, t) = -1, t\geq 0\\ &\\ & u(1, t) = 0, t\geq 0\end{cases}$

C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
#define _USE_MATH_DEFINES 
#include <cmath>
#include <locale>
#include <fstream>
#include <iostream>
#include <vector>
using namespace std;
 
 
int N = 40; //Число узлов
double h = 1.0 / double(N - 1); //Шаг
double tau = h * h / 4.0; //Шаг по времени
vector<vector<double>>Ui; //Значения численного решения на слоях
vector<double>R(N - 1); //Правая часть системы
vector<double>progonka(N); //Резы прогонки
vector<double>alfa(N - 1); //Прогоночный коэф-т
vector<double>beta(N - 1); //Прогоночный коэф-т
vector<double>Xi(N); //Узлы
vector<double>first(N); //Сначала значение f(x), потом рез-ы прогонки
double e = 1e-7; //Разницы между слоями численного решения для сходимости
 
 
double C = (1.0 / (12.0 * tau) - 1.0 / (2.0 * h * h)); 
double B = (5.0 / (6.0 * tau) + 1.0 / (h * h));
double D = (1.0 / (12.0 * tau) + 1.0 / (2.0 * h * h));
double F = (5.0 / (6.0 * tau) - 1.0 / (h * h));
double G = D, A = C;
 
void initXi() {
    Xi[0] = 0;
    for (int i = 1; i < Xi.size(); i++)
        Xi[i] = Xi[i - 1] + h;
}
double Sk(double x, double tau, int k) {
    double Ak = -48.0 * (4.0 * pow(-1, k) * M_PI * (2.0 * k - 1.0) + 8.0) / pow((M_PI) * (2.0 * k - 1.0), 4);
    double T = exp(-tau * pow((M_PI * (2.0 * k - 1.0) / 2.0), 2));
    double X = cos(M_PI * (2.0 * k - 1.0) * x / 2.0);
    return Ak * T * X;
}
double solution(double x, double tau) {
    int k = 1;
    double alfa = 1.0 - x;
    while (abs(Sk(x, tau, k)) > 0.00000000001) {
        alfa += Sk(x, tau, k);
        k++;
    }
    return alfa;
}
 
double f(double x) {
    return -2.0 * pow(x, 3) - x + 3.0;
}
 
//Проверка сходимости решения
int check(vector<double>a, vector<double>b) {
 
    for (int i = 0; i < a.size(); i++) {
        if (abs(a[i] - b[i]) > e)
            return 0;
    }
    return 1;
}
 
//Правая часть
vector<double>Ri(vector<double> k) {
    //Второй порядок
    R[0] = (-h * 2.0 + (k[0] * D + k[1] * F + k[2] * G) / C);
    
    //Первый порядок
    //R[0] = h; 
 
    for (int i = 1; i < R.size(); i++)
        R[i] = (k[i - 1] * D + k[i] * F + k[i + 1] * G);
 
    return R;
}
 
//Прогонка
vector<double> sweep(vector<double> rig) {
    //Первый порядок
    //alfa[0] = 1.0;
    //beta[0] = rig[0];
 
 
    //Второй порядок
    alfa[0] = (-(4.0 + B / C)) / (-3.0 + A / C);
    beta[0] = rig[0] / (-3.0 + A / C);
 
    for (int i = 1; i <= alfa.size() - 1; i++) {
        alfa[i] = -C / (B + alfa[i - 1] * A);
        beta[i] = (rig[i] - beta[i - 1] * A) / (B + alfa[i - 1] * A);
    }
 
    progonka[progonka.size() - 1] = 0;
 
    for (int i = progonka.size() - 2; i >= 0; i--) {
        progonka[i] = (alfa[i] * progonka[i + 1] + beta[i]);
    }
 
    return progonka;
}
 
 
 
int main()
{
    ofstream Sol("sol.txt"); //Для аналитического решения
    ofstream Progon("progon.txt"); //Для численного конечного
 
    //Для численного по слоям
    ofstream Progon1("progon1.txt");
    ofstream Progon2("progon2.txt");
    ofstream Progon3("progon3.txt");
    ofstream Progon4("progon4.txt");
    setlocale(LC_ALL, "rus");
 
    initXi();
 
    for (int i = 0; i < first.size(); i++) {
        first[i] = f(Xi[i]);
    }
    Ui.push_back(Xi);
    Ui.push_back(first);
    
 
    int i = 1;
    do {
        R = Ri(first);
        progonka = sweep(R);
        first = progonka;
        Ui.push_back(progonka);
        i++;
    } while (check(Ui[i], Ui[i - 1]) == 0);
 
    double t = tau * i; //Затраченное время для сходимости
    cout << "Шаг: h = " << h << endl;
    cout << "Время: T = " << t << endl << endl;
 
    //Максимальная ошибка 
    double tp = 0;
    double errt = -1;
    for (int j = 1; j < Ui.size(); j++) {
        double err = -1;
        for (int k = 0; k < N; k++) {
            if (abs(solution(Xi[k], tp) - Ui[j][k]) > err)
                err = abs(solution(Xi[k], tp) - Ui[j][k]);
        }
        if (err > errt)
            errt = err;
        tp += tau;
    }
    cout << "N = " << N << "\nОшибка: " << errt << endl;
    
    for (int j = 0; j < N; j++) {
        Progon1 << Xi[j] << ' ' << Ui[4][j] << endl;
        Progon2 << Xi[j] << ' ' << Ui[21][j] << endl;
        Progon3 << Xi[j] << ' ' << Ui[Ui.size()/25 - 1][j] << endl;
        Progon4 << Xi[j] << ' ' << Ui[Ui.size()/10][j] << endl;
        Progon << Xi[j] << ' ' << Ui[Ui.size() - 1][j] << endl;
        Sol << Xi[j] << ' ' << solution(Xi[j], t) << endl;
    }
}

Новые блоги и статьи Все статьи Все блоги /
Модель заражения группы наркоманов alhaos 17.04.2026 Условия задачи сформулированы тут Суть: - Группа наркоманов из 10 человек. - Только один инфицирован ВИЧ. - Колются одной иглой. - Колются раз в день. - Колются последовательно через. . .	Мысли в слух. Про "навсегда". kumehtar 16.04.2026 Подумалось тут, что наверное очень глупо использовать во всяких своих установках понятие "навсегда". Это очень сильное понятие, и я только начинаю понимать край его смысла, не смотря на то что давно. . .	My Business CRM MaGz GoLd 16.04.2026 Всем привет, недавно возникла потребность создать CRM, для личных нужд. Собственно программа предоставляет из себя базу данных клиентов, в которой можно фиксировать звонки, стадии сделки, а также. . .	Знаешь почему 90% людей редко бывают счастливыми? kumehtar 14.04.2026 Потому что они ждут. Ждут выходных, ждут отпуска, ждут удачного момента. . . а удачный момент так и не приходит.
Фиксация колонок в отчете СКД Maks 14.04.2026 Фиксация колонок в СКД отчета типа Таблица. Задача: зафиксировать три левых колонки в отчете. Процедура ПриКомпоновкеРезультата(ДокументРезультат, ДанныеРасшифровки, СтандартнаяОбработка) / / . . .	Настройки VS Code Loafer 13.04.2026 { "cmake. configureOnOpen": false, "diffEditor. ignoreTrimWhitespace": true, "editor. guides. bracketPairs": "active", "extensions. ignoreRecommendations": true, . . .	Оптимизация кода на разграничение прав доступа к элементам формы Maks 13.04.2026 Алгоритм из решения ниже реализован на нетиповом документе, разработанного в конфигурации КА2. Задачи, как таковой, поставлено не было, проделанное ниже исключительно моя инициатива. Было так:. . .	Контроль заполнения и очистка дат в зависимости от значения перечислений Maks 12.04.2026 Алгоритм из решения ниже реализован на примере нетипового документа "ПланированиеПерсонала", разработанного в конфигурации КА2. Задача: реализовать контроль корректности заполнения дат назначения. . .