Lab14/Lab_14/Source.c - master - aytsan_ns/Num_meth

Lab14/Lab_14/Source.c
403 строки10 KB
aytsan-ns
14 работа
30 сен 2024, 19:49
30 сен 2024, 19:49310e36c
#define _CRT_SECURE_NO_WARNINGS
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#define PI 3.14159265358979323846

double p(double x) {
    return 1;
}

double q(double x) {
    return -sin(x);
}

double r(double x) {
    return cos(x);
}

double f(double x) {
    return 1 - cos(x);
}

double F(double x) {
    return cos(x);
}

double* TridiagonalSolve(int n, double** A, double* d) {
    
    double* a = (double*)malloc(n * sizeof(double));
    double* b = (double*)malloc(n * sizeof(double));
    double* c = (double*)malloc(n * sizeof(double));
    double* x = (double*)malloc(n * sizeof(double));

    if (a == NULL || b == NULL || c == NULL || x == NULL) {
        printf("Error");
        return NULL;
    }

    for (int i = 0; i < n; i++) {
        b[i] = A[i][i];
        if (i > 0) {
            a[i] = A[i][i - 1];
        }
        if (i < n - 1) {
            c[i] = A[i][i + 1]; 
        }
    }

    double* P = (double*)malloc(n * sizeof(double));
    double* Q = (double*)malloc(n * sizeof(double));

    if (P == NULL || Q == NULL) {
        printf("Error");
        return NULL;
    }

    P[0] = -c[0] / b[0];
    Q[0] = d[0] / b[0];

    for (int i = 1; i < n; i++) {
        double denominator = b[i] + a[i] * P[i - 1];
        P[i] = (i < n - 1) ? -c[i] / denominator : 0;
        Q[i] = (d[i] - a[i] * Q[i - 1]) / denominator;
    }

    x[n - 1] = Q[n - 1];
    for (int i = n - 2; i >= 0; i--) {
        x[i] = P[i] * x[i + 1] + Q[i];
    }

    free(a);
    free(b);
    free(c);
    free(Q);
    free(P);

    return x;
}

/*
double* GaussSolve(int n, double** A, double* b) {
    double factor;
    double* x = (double*)malloc(n * sizeof(double));
    if (x == NULL) {
        printf("Error");
        return NULL;
    }

    for (int i = 0; i < n; i++)
        x[i] = 0;
    
    for (int k = 0; k < n - 1; k++) {
        for (int i = k + 1; i < n; i++) {
            factor = A[i][k] / A[k][k];
            for (int j = k; j < n; j++) {
                A[i][j] -= factor * A[k][j];
            }
            b[i] -= factor * b[k];
        }
    }

    x[n - 1] = b[n - 1] / A[n - 1][n - 1];
    for (int i = n - 2; i >= 0; i--) {
        x[i] = b[i];
        for (int j = i + 1; j < n; j++) {
            x[i] = x[i] - A[i][j] * x[j];
        }
        x[i] /= A[i][i];
    }

    return x;
}*/

double* FiniteDifference(double a, double b, int N) {
    double h = (b - a) / (N - 1);
    double y0 = 1;
    double x1 = a + h;

    double** A = (double**)malloc((N - 1) * sizeof(double*));
    for (int i = 0; i < N - 1; i++) {
        A[i] = (double*)malloc((N - 1) * sizeof(double));
    }

    double* B = (double*)malloc((N - 1) * sizeof(double));
    if (B == NULL) {
        printf("Error3");
        return NULL;
    }

    for (int k = 1; k < N - 1; k++) {
        double x = a + k * h;
        if (x == x1) {
            double A1 = p(x) * (-2 / (h * h)) + r(x);
            double A2 = p(x) * (1 / (h * h)) + q(x) / (2 * h);
            B[k - 1] = f(x) - y0 * (p(x) * (1 / (h * h)) - q(x) / (2 * h));

            for (int i = 0; i < N - 1; i++)
                A[k - 1][i] = 0;
            A[k - 1][k - 1] = A1;
            A[k - 1][k] = A2;
        }
        else {
            if (x == b) {
                double A0 = -1;
                double A2 = 1;
                B[k - 1] = -2*h;

                for (int i = 0; i < N - 1; i++)
                    A[k - 1][i] = 0;
                A[k - 1][k - 2] = A0;
                A[k - 1][k] = A2;
            }
            else {
                double A0 = p(x) * (1 / (h * h)) - q(x) / (2 * h);
                double A1 = p(x) * (-2 / (h * h)) + r(x);
                double A2 = p(x) * (1 / (h * h)) + q(x) / (2 * h);
                B[k - 1] = f(x);
                for (int i = 0; i < N - 1; i++)
                    A[k - 1][i] = 0;
                A[k - 1][k - 1] = A1;
                A[k - 1][k - 2] = A0;
                A[k - 1][k] = A2;
            }
        }
    }
    double* yn = TridiagonalSolve(N - 1, A, B);
    
    double* y = (double*)malloc(N * sizeof(double));
    if (y == NULL) {
        printf("Error");
        return NULL;
    }
    y[0] = y0;
    for (int i = 1; i < N; i++) {
        y[i] = yn[i - 1];
    }

    free(B);
    for (int i = 0; i < N - 1; i++) {
        free(A[i]);
    }
    free(A);
    free(yn);

    return y;
}

double* CreateMesh(double a, double b, int N) {
    double* points = (double*)malloc(sizeof(double) * N);
    if (points == NULL) {
        printf("Error");
        return NULL;
    }
    for (int i = 0; i < N; i++) {
        points[i] = (a + b) / 2.0 + (b - a) / 2.0 * sin(PI * (2.0 * i - N + 1) / (2.0 * (N - 1)));
    }
    double* mesh = (double*)malloc(sizeof(double) * N);
    if (mesh == NULL) {
        printf("Error");
        return NULL;
    }
    mesh[0] = points[0];
    mesh[N - 1] = points[N - 1];
    int k = (N - 1) / 2;
    for (int i = 1; i < (N / 2); i++) {
        mesh[i] = a + (b - a) / 2 * (sin(i * PI / N));
        mesh[N - 1 - i] = b - (b - a) / 2 * (sin(i * PI / N));
    }
    free(points);

    return mesh;
}

double* FiniteDifferenceForMesh(double a, double b, int N, double* mesh) {
    double* h = (double*)malloc(sizeof(double) * (N - 1));
    if (h == NULL) {
        printf("Error1");
        return NULL;
    }
    for (int i = 0; i < (N - 1); i++) {
        h[i] = mesh[i + 1] - mesh[i];
    }

    double y0 = 1;
    double x1 = mesh[1];

    double** A = (double**)malloc((N - 1) * sizeof(double*));
    for (int i = 0; i < N - 1; i++) {
        A[i] = (double*)malloc((N - 1) * sizeof(double));
    }

    double* B = (double*)malloc((N - 1) * sizeof(double));
    if (B == NULL) {
        printf("Error2");
        return NULL;
    }

    for (int k = 1; k < N - 1; k++) {
        double x = mesh[k];
        if (x == x1) {
            double A1 = p(x) * (-2 / (h[k - 1] * h[k])) + r(x);
            double A2 = p(x) * (1 / (h[k - 1] * h[k])) + q(x) / (h[k] + h[k - 1]);
            B[k - 1] = f(x) - y0 * (p(x) * (1 / (h[k] * h[k - 1])) - q(x) / (h[k - 1] + h[k]));

            for (int i = 0; i < N - 1; i++)
                A[k - 1][i] = 0;
            A[k - 1][k - 1] = A1;
            A[k - 1][k] = A2;
        }
        else {
            if (x == b) {
                double A0 = 1;
                double A1 = -1;
                B[k - 1] = -h[k-1];

                for (int i = 0; i < N - 1; i++)
                    A[k - 1][i] = 0;
                A[k - 1][k - 2] = A0;
                A[k - 1][k - 1] = A1;
            }
            else {
                double A0 = p(x) * (1 / (h[k] * h[k - 1])) - q(x) / (h[k] + h[k - 1]);
                double A1 = p(x) * (-2 / (h[k] * h[k - 1])) + r(x);
                double A2 = p(x) * (1 / (h[k] * h[k - 1])) + q(x) / (h[k] + h[k - 1]);
                B[k - 1] = f(x);
                for (int i = 0; i < N - 1; i++)
                    A[k - 1][i] = 0;
                A[k - 1][k - 1] = A1;
                A[k - 1][k - 2] = A0;
                A[k - 1][k] = A2;
            }
        }
    }
    double* yn = TridiagonalSolve(N - 1, A, B);

    double* y = (double*)malloc(N * sizeof(double));
    if (y == NULL) {
        printf("Error");
        return NULL;
    }
    y[0] = y0;
    for (int i = 1; i < N; i++) {
        y[i] = yn[i - 1];
    }

    free(h);
    free(B);
    for (int i = 0; i < N - 1; i++) {
        free(A[i]);
    }
    free(A);
    free(yn);

    return y;
}
/*
double ComputeErrorEps(double a, double b, double eps) {
    int N = 2;
    double* y1;
    double* y2;
    y1 = FiniteDifference(a, b, N);
    do {
        N *= 2;
        y2 = FiniteDifference(a, b, N);

        double max_error = 0;
        for (int i = 0; i < N / 2; i++) {
            double error = fabs(y1[i] - y2[2 * i]);
            if (error > max_error) max_error = error;
        }

        if (max_error < eps) return max_error;

        y1 = y2;
    } while (1);
}
*/
void Task1() {
    double a = 0, b = PI / 2;
    int N1 = 4, N2 = 7;
    double h1 = (b - a) / (N1 - 1);
    double h2 = (b - a) / (N2 - 1);

    double* y1 = FiniteDifference(a, b, N1);
    double* y2 = FiniteDifference(a, b, N2);

    FILE* result1 = fopen("result1.txt", "w");
    FILE* result2 = fopen("result2.txt", "w");
    FILE* error1 = fopen("error1.txt", "w");
    FILE* error2 = fopen("error2.txt", "w");

    for (int i = 0; i < N1; i++) {
        double x = a + i * h1;
        fprintf(result1, "%.10lf %.10lf\n", x, y1[i]);
        fprintf(error1, "%.10lf %.10lf\n", x, fabs(y1[i] - F(x)));
    }
    for (int i = 0; i < N2; i++) {
        double x = a + i * h2;
        fprintf(result2, "%.10lf %.10lf\n", x, y2[i]);
        fprintf(error2, "%.10lf %.10lf\n", x, fabs(y2[i] - F(x)));
    }

    _fcloseall();
}

void Task2() {
    double a = 0, b = PI / 2;
    int N = 20;
    double* mesh = CreateMesh(a, b, N);

    double* y = FiniteDifferenceForMesh(a, b, N, mesh);

    FILE* result = fopen("result.txt", "w");
    FILE* error = fopen("error.txt", "w");

    for (int i = 0; i < N; i++) {
        double x = mesh[i];
        fprintf(result, "%.10lf %.10lf\n", x, y[i]);
        fprintf(error, "%.10lf %.10lf\n", x, fabs(y[i] - F(x)));
    }
    _fcloseall();
}

void Task3() {
    double a = 0, b = PI / 2;
    
    FILE* result_error_eps = fopen("result_error_eps.txt", "w");

    int N = 2;
    double* y1;
    double* y2;
    y1 = FiniteDifference(a, b, N);

    for (double eps = 0.1; eps>= 1e-4; eps/=10) {
        do {
            N *= 2;
            y2 = FiniteDifference(a, b, N);

            double max_error = 0;
            for (int i = 0; i < N / 2; i++) {
                double error = fabs(y1[i] - y2[2 * i]) / 3;
                if (error > max_error) max_error = error;
            }

            if (max_error < eps) {
                fprintf(result_error_eps, "%.10lf %.10lf\n", eps, max_error);
                y1 = y2;
                printf("%i ", N);
                break;
            }

            y1 = y2;
        } while (1);
    }
    _fcloseall();
}

void main() {
    //Task1();
    //Task2();
    Task3();
}