Численные методы. Волновое уравнение и уравнение теплопроводности. Метод сеток

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
 
#define _SQR_(X) X*X
 
const double  XMIN = 0.0;
const double  XMAX= 1.0;
const double  YMIN = 0.0;
const double  YMAX = 1.0;
 
const int NX = 100;
const int NY = 100;
 
const double TAU = 0.0001;  // шаг по времени
const double TMAX = 10.0;    // время, до которого будут идти вычисления
 
const int SAVE_STEP = 1000;    // через сколько шагов записываем файлы с результатами
 
double hx, hy;
 
/**
 *  одномерные массивы размерности NX*NY,
 *  [i,j] -> [j*NX+i] (в памяти элементы располагаем по строкам)
 */
double *u;          // значения на текущем временном слое
double *u1;         // значения на следующем временном слое
 
/**
 *  граничные условия (элементам вдоль границы присваиваем 0 или 1)
 */
void bound_cond()
{
    for (int i = 0; i < NX; i++)
    {
        u[i]            = 0.0;
        u[(NY-1)*NX+i]  = 0.0;
    }
    for (int j = 0; j < NY; j++)
    {
        u[j*NX]         = 1.0;
        u[j*NX+NX-1]    = 0.0;
    }
}
 
/**
 *  запись результатов в текстовый файл
 */
void save(int step)
{
    char fname[50];
    sprintf(fname, "res_%010d.txt", step);
    FILE * f = fopen(fname, "w");
    for (int i = 0; i < NX; i++)
    {
        double x = XMIN+i*hx;
        for (int j = 0; j < NY; j++)
        {
            double y = YMIN+j*hy;
            fprintf(f, "%25.26f\t%25.26f\t%25.26f\n", x, y, u[j*NX+i]);
        }
    }
    fclose(f);
    printf("File '%s' saved...\n", fname);
}
 
int main(int argc, char** argv)
{
    u  = (double*)malloc(sizeof(double)*NX*NY);
    u1 = (double*)malloc(sizeof(double)*NX*NY);
 
 
    hx = (XMAX-XMIN)/NX;
    hy = (YMAX-YMIN)/NY;
 
    // начальные условия (заполняем массив нулями, т.к. температура = 0)
    memset(u, 0, sizeof(double)*NX*NY);
    bound_cond();
 
    double t = 0.0;
    int step = 0;
    while (t < TMAX)
    {
        t += TAU; step++;
        // вычисляем значения во внутренних ячейках
        for (int i = 1; i < NX-1; i++)
        {
            for (int j = 1; j < NY-1; j++)
            {
                u1[j*NX+i] = u[j*NX+i]+TAU*( (u[j*NX+i+1]-2.0*u[j*NX+i]+u[j*NX+i-1])/_SQR_(hx)+
                                              (u[(j+1)*NX+i]-2.0*u[j*NX+i]+u[(j-1)*NX+i])/_SQR_(hy) );
            }
        }
        memcpy(u, u1, sizeof(double)*NX*NY); // копируем массив, соответствующий следующему времени, в текущий массив
        bound_cond();
        if (step % SAVE_STEP == 0) save(step);
    }
 
    free(u);
    free(u1);
    return 0;
}