12/code/Chisl2sem3homework.cpp - master - levg/numm

12/code/Chisl2sem3homework.cpp
117 строк3 KB
Лев Черняховский
12-4
07 ноя 2024, 16:43
07 ноя 2024, 16:4398e9b73
#include <iostream>
#include <fstream>
#include <cmath>
#include <vector>
#include <algorithm>
#include <string>

using namespace std;

// Функция для уравнения: x^2 * y' + x * y + 1 = 0
double f(double x, double y) {
    return -(x * y + 1) / (x * x);
}

// Точное решение: y = (1 - ln(|x|)) / x
double exact_solution(double x) {
    return (1 - log(abs(x))) / x;
}

void fixed_step_method(double x0, double y0, double h, const char* filename) {
    ofstream file(filename);
    file << "x,y" << endl;

    while (x0 <= 3) {
        file << x0 << "," << y0 << endl;

        double k1 = h * f(x0, y0);
        double k2 = h * f(x0 + h / 2, y0 + k1 / 2);
        y0 += k2;
        x0 += h;
    }

    file.close();
}

double adaptive_euler(double x0, double y0, double h, double tol, const char* filename) {
    ofstream file(filename);
    file << "x,y,h,error" << endl;

    double max_error = 0.0;
    double f_x0_y0 = f(x0, y0);  // Начальное значение функции f(x0, y0)

    while (x0 <= 3) {
        // Вычисление большого шага
        double k1_big = h * f_x0_y0;
        double y_half_big = y0 + k1_big / 2;

        double f_x_half_big = f(x0 + h / 2, y_half_big);
        double k2_big = h * f_x_half_big;
        double y_big_step = y0 + k2_big;

        // Половинный шаг - первый подшаг
        double h_half = h / 2;
        double k1_small1 = h_half * f_x0_y0;
        double y_half_small1 = y0 + k1_small1 / 2;

        double f_x_quarter = f(x0 + h_half / 2, y_half_small1);
        double k2_small1 = h_half * f_x_quarter;
        double y_mid = y0 + k2_small1;

        // Половинный шаг - второй подшаг
        double f_x_half_small = f(x0 + h_half, y_mid);
        double k1_small2 = h_half * f_x_half_small;

        double f_x_final = f(x0 + h, y_mid + k1_small2 / 2);
        double k2_small2 = h_half * f_x_final;
        double y_small_step = y_mid + k2_small2;

        double error = abs(y_small_step - y_big_step) / 3.0;
        max_error = max(max_error, error);

        if (error > tol) {
            h *= 0.5;
        }
        else {
            file << x0 << "," << y0 << "," << h << "," << error << endl;

            x0 += h;
            y0 = y_small_step;
            if (error < tol / 4.0) {
                h *= 2.0;
            }

            // Пересчитываем значение f(x0, y0) для новых x0 и y0
            f_x0_y0 = f(x0, y0);
        }
    }

    file.close();
    return max_error;
}

int main() {
    double x0 = 1.0;
    double y0 = 1.0;

    //fixed_step_method(x0, y0, 0.1, "fixed_step_h1.csv");
    //fixed_step_method(x0, y0, 0.05, "fixed_step_h2.csv");

    double initial_h = 0.1;

    ofstream result_file("error_vs_tolerance.csv");
    result_file << "tolerance,max_error" << endl;
    vector<double> tolerances = { 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2 };

    for (double tol : tolerances) {
        string filename = "adaptive_euler_tol_" + to_string(tol)+ ".csv";
        double max_error = adaptive_euler(x0, y0, initial_h, tol, filename.c_str());

        result_file << tol << "," << max_error << endl;
        cout << "Completed for tolerance = " << tol << ", max error = " << max_error << endl;
    }

    result_file.close();

    return 0;
}