скx

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
/*double Q(double x) // построение источника
{
 if (abs(x - 0.5) < 0.1)
 return 1;
 else
 return 0;
}*/
void appr(int n, double *a, double *b, double *c, double *f)
//построение аппроксимации
{
 for (int i = 0; i < n - 1; i++)
 {
 a[i] = 1.0*n*n;
 b[i] = 1.0*n*n;
 c[i] = - 2.0*n*n;
 f[i] = 0;
 }
 c[0] = -1;
 c[n - 1] = -3.0*n*n; //???
 f[0] = - 1.0 / (n) *n*n;
}
void invert(int n, double *d, double *inv, double *a, double *b, double
*c, double *q, double *p) // умножение обратной к ТД-матрице на вектор
{
 q[0] = -b[0]/c[0];
 p[0] = d[0]/c[0];
 for (int i = 1; i < n-1; i++) // прямой ход прогонки
 {
 q[i] = -b[i] / (c[i] + a[i-1] * q[i-1]);
 p[i] = (d[i] - a[i-1] * p[i-1]) / (c[i] + a[i-1] * q[i-1]);
 }
 inv[n-1] = (d[n-1] - p[n-2] * a[n-2])/(c[n-1] + a[n-2] * q[n-2]);
// обратный ход прогонки
 for (int i = n-2; i > -1; i--)
 inv[i] = q[i] * inv[i+1] + p[i];
}
int main()
{
 int n = 150; //число шагов сетки по х
 int nt = 225000; //число шагов сетки по времени
 double t1 = 2;
 double t2 = 5; // заданные моменты времени
 double dt = t2 / nt;
 double *u_0, *u, *a, *b, *c, *f, *ut, *a_new, *b_new, *c_new, *d,
*inv, *q, *p;
 // трёхдиагональная матрица аппроксимации:
 a = new double[n - 1];
 b = new double[n - 1];
 c = new double[n];
// модифицированная матрица
 a_new = new double[n - 1];
 b_new = new double[n - 1];
 c_new = new double[n];
 q = new double[n]; // вспомогательные
 p = new double[n]; // коэффициенты для прогонки
 f = new double[n]; // функция источника + гр. условий на сетке
 d = new double[n]; // вспомогательные
 inv = new double[n]; // векторы
 u_0 = new double[n]; //начальное распределение
 ut = new double[nt];
 for (int i = 0; i < n; i++)
 u_0[i] = 0;
 appr(n, a, b, c, f);
 u = u_0; //текущее распределение по х
 for (int k = 0; k < nt; k++)
 {
 ut[k] = (u[n / 2] + u[n / 2 - 1]) / 2; // температура в середине
стержня в данный момент
 for (int i = 0; i < n - 1; i++) // делаем новые матрицу и
вектор, чтобы проще записывать обращение
 {
 d[i] = u[i] + dt*f[i]; // новый вектор
 a_new[i] = -dt*a[i]/2; // новая матрица A' = E - dt/2*A,
 b_new[i] = -dt*b[i]/2; // где A(a,c,b) - ТД-матрица
аппроксимации
 c_new[i] = 1 - dt*c[i]/2;
 }
 c_new[n-1] = 1 - dt*c[n-1]/2;
 d[n-1] = u[n-1] + dt*f[n-1];
 invert(n, d, inv, a_new, b_new, c_new, q, p); // домножаем d на
обратную к A'
for (int i = 0; i < n; i++)
 u[i] = 2*inv[i] - d[i];
 if (abs(k*dt - t1) < dt / 2)
 {
 ofstream u_t1_out("CrNi_t1.txt");
 for (int i = 0; i < n; i++)
 u_t1_out << (i + 0.5) / n << ' ' << u[i] << '\n';
 u_t1_out.close();
 }
 }
 ofstream u_t2_out("CrNi_t2.txt"); // распределение по стержню в момент
времени t2
 ofstream u_x_out("CrNi_x.txt"); // распределение по времени в середине
стержня
 for (int k = 0; k < nt; k++)
 u_x_out << k * dt << ' ' << ut[k] << '\n';
 for (int i = 0; i < n; i++)
 u_t2_out << (i + 0.5) / n << ' ' << u[i] << '\n';
 u_x_out.close();
 u_t2_out.close();
 return 0;
}