McKormack new

#define _CRT_SECURE_NO_WARNINGS
#include <iostream>
#include<stdio.h>
#include <math.h>
#include <fstream>
using namespace std;

int N;
double h, tau;

//Точное решение

double U(double X, double T) {
    if (X + T <= 1) {
        if (T >= X) return 1;
        return (1 - 2 * X);
    }
    else if (X + T > 1)
        return -1;
    else return 0;
}

double max(double *u2[], int N, int M) {
    double max = u2[0][0];
    for (int i = 1; i < N + 1; i++) {
        for (int j = 0; j < M + 1; j++)
            if (max < u2[i][j]) {
                max = u2[i][j];
            }
    }
    return max;
}

double f2(double *u[], int N, int M, double h, double tau) {
    u[0][0] = 1;//точка разрыва (для гладкого условия просто значение 1)
    u[N][M] = 1;//точка разрыва, перенесённая по характеристике на правую границу
    for (int i = 1; i < N + 1; i++) {//Н.У
        u[i][0] = cos(3.14159*(i*h));
        //u[i][0] = 0;
    }
    for (int i = 1; i < M + 1; i++) {//К.У слева
        u[0][i] = exp(-i * tau);
        //u[0][i] = 1;
    }
    for (int i = 1; i < M; i++) {// К.У справа
        u[N][i] = cos(3.14159*(1 - i * tau));
        //u[N][i] = 0;
    }
    for (int i = 1; i < N; i++) {// двухслойная схема для второго (первого после нулевого) слоя
        u[i][1] = u[i - 1][0] + (tau - h) / (h + tau)*u[i][0] - (tau - h) / (h + tau)*u[i - 1][1];
    }
    double **u2 = new double*[N + 1];//погрешность в узлах
    for (int k = 0; k < N + 1; k++) {
        u2[k] = new double[M + 1];
    }
    for (int i = 0; i < N + 1; i++) {
        for (int j = 0; j < M + 1; j++) {
            u2[i][j] = fabs(U(i, j) - u[i][j]);
        }
    }
    double maxpogr;
    maxpogr = max(u2, N, M);
    for (int d = 0; d < N + 1; d++) {
        delete[] u2[d];
    }
    return maxpogr;
}

//Функция дл создания массива
double **cr_arr(size_t a, size_t b) {
    double **m = new double*[a];
    m[0] = new double[a*b];
    for (size_t i = 1; i != a; ++i) {
        m[i] = m[i - 1] + b;

    }
    return m;
}

//Заполняем массив чиселками
double ** fill_ar() {
    double **u = cr_arr(N + 1, N + 1);
    //у1-с волной, у2 - с чертой
    double *u1 = new double[N + 1], *u2 = new double[N + 1];
    for (int i = 0; i < N + 1; i++) {//Н.У
        u[i][0] = -2 * i * h + 1;
    }
    for (int i = 1; i < N + 1; i++) {//К.У слева
        u[0][i] = 1;
    }
    for (int i = 1; i < N + 1; i++) {// К.У справа
        u[N][i] = -1;
    }
    for (int m = 0; m < N; m++) {
        u1[0] = tau * (u[0][m] * u[0][m] - u[1][m] * u[1][m]) / (2.*h) + u[0][m];
        for (int j = 1; j < N; j++) {
            u1[j] = tau * (u[j][m] * u[j][m] - u[j + 1][m] * u[j + 1][m]) / (2.*h) + u[j][m];
            u2[j] = u[j][m] - tau * u1[j] * u1[j] / (2.*h) + tau * u1[j - 1] * u1[j - 1] / (2.*h);
            u[j][m + 1] = (u1[j] + u2[j]) / 2.;
        }
    }
    return(u);
}

int main() {
    FILE *ux8, *ut8, *ux, *ut, *ut_err, *ux_err;
    ux8 = fopen("ux8.txt", "w+");
    ut8 = fopen("ut8.txt", "w+");
    ux = fopen("ux.txt", "w+");
    ut = fopen("ut.txt", "w+");
    ut_err = fopen("ut_err.txt", "w+");
    ux_err = fopen("ux_err.txt", "w+");
    double **u, ux_er = 0, ut_er = 0, u_err = 0;
    N = 10;
    for (int k = 0; k < 3; k++) {
        h = 1 / (double)(N), tau = 1 / (double)N;
        u = fill_ar();
        for (int i = 0; i <N + 1; i++) {
            for (int j = 0; j < N + 1; j++) {
                if (h*i == 0.5)
                    fprintf(ut, "%0.5lf %0.5lf \n", j*tau, u[i][j]);
                if (tau*j == 0.5)
                    fprintf(ux, "%0.5lf %0.5lf \n", i*h, u[i][j]);
                if (tau*j == 0.8)
                    fprintf(ux8, "%0.5lf %0.5lf \n", i*h, u[i][j]);
                if (h*i == 0.8)
                    fprintf(ut8, "%0.5lf %0.5lf \n", j*tau, u[i][j]);
                if (fabs(u[i][j] - U(i*h, j*tau)) >= u_err)
                    u_err = fabs(u[i][j] - U(i*h, j*tau));
            }
            ux_er = u_err;
        }
        for (int j = 0; j < N + 1; j++) {
            for (int i = 0; i < N + 1; i++) {
                if (fabs(u[i][j] - U(i*h, j*tau)) >= u_err)
                    u_err = fabs(u[i][j] - U(i*h, j*tau));
            }
            ut_er = u_err;
        }
        fprintf(ux_err, "N = %d  %0.5lf\n", N, ux_er);
        fprintf(ut_err, "N = %d  %0.5lf\n", N, ut_er);
        N = N * 10;

    }
    fclose(ux);
    fclose(ut);
    fclose(ux8);
    fclose(ut8);
    return 0;
}