14/code/Chisl2sem4homework.cpp - master - levg/numm

14/code/Chisl2sem4homework.cpp
183 строки6 KB
Лев Черняховский
14-1
10 ноя 2024, 03:55
10 ноя 2024, 03:5591c0a7b
#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <string>

const double a = 0.2;
const double b = 1.0;

double alpha0 = 1.0, alpha1 = 1.0;
double A = 6.0;
double beta0 = 1.0;
double B = 2.0;

double u_a = 0.0;
double du_a = A / alpha1;
double v_a = A / alpha0;
double dv_a = 0.0;

struct Solution {
    double x;
    double y;
    double step_size;
};

double p(double x) { return x * x * (x + 1); }
double q(double x) { return -1; }
double r(double x) { return -2; }
double f(double x) { return 1 / (x * x); }

double exactSolution(double x) {
    return 1 + 1 / x;
}

double ode(double x, double y, double dy, bool homogeneous) {
    if (homogeneous) {
        return (-q(x) * dy - r(x) * y) / p(x);
    }
    else {
        return (f(x) - q(x) * dy - r(x) * y) / p(x);
    }
}

std::vector<Solution> rungeKuttaFixed(double x0, double y0, double dy0, double h, bool homogeneous) {
    std::vector<Solution> result;
    double x = x0, y = y0, dy = dy0;

    while (x <= b) {
        double k1 = h * dy;
        double l1 = h * ode(x, y, dy, homogeneous);
        double k2 = h * (dy + 0.5 * l1);
        double l2 = h * ode(x + 0.5 * h, y + 0.5 * k1, dy + 0.5 * l1, homogeneous);
        double k3 = h * (dy + 0.5 * l2);
        double l3 = h * ode(x + 0.5 * h, y + 0.5 * k2, dy + 0.5 * l2, homogeneous);
        double k4 = h * (dy + l3);
        double l4 = h * ode(x + h, y + k3, dy + l3, homogeneous);

        double y_new = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
        double dy_new = dy + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;

        result.push_back({ x, y, h });
        x += h;
        y = y_new;
        dy = dy_new;
    }
    return result;
}

void writeToFile(const std::string& filename, const std::vector<Solution>& data, bool include_error) {
    std::ofstream file(filename);
    file << "x\ty\ty_exact\terror\tstep_size\n";
    for (const auto& entry : data) {
        double y_exact = exactSolution(entry.x);
        double error = include_error ? std::abs(entry.y - y_exact) : 0.0;
        file << std::fixed << std::setprecision(6)
            << entry.x << "\t" << entry.y << "\t" << y_exact << "\t"
            << error << "\t" << entry.step_size << "\n";
    }
    file.close();
}

std::vector<Solution> rungeKuttaAdaptive(double x0, double y0, double dy0, double tol, bool homogeneous) {
    std::vector<Solution> result;
    double x = x0, y = y0, dy = dy0;
    double h = 0.1;

    double ode1= ode(x, y, dy, homogeneous);
    while (x < b) {
        double k1 = h * dy;
        double l1 = h * ode1;
        double k2 = h * (dy + 0.5 * l1);
        double l2 = h * ode(x + 0.5 * h, y + 0.5 * k1, dy + 0.5 * l1, homogeneous);
        double k3 = h * (dy + 0.5 * l2);
        double l3 = h * ode(x + 0.5 * h, y + 0.5 * k2, dy + 0.5 * l2, homogeneous);
        double k4 = h * (dy + l3);
        double l4 = h * ode(x + h, y + k3, dy + l3, homogeneous);

        double y_full = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
        double dy_full = dy + (l1 + 2 * l2 + 2 * l3 + l4) / 6.0;

        double half_h = h / 2.0;
        double y_half = y, dy_half = dy;

        for (int i = 0; i < 2; ++i) {
            double k1_half = half_h * dy_half;
            double l1_half = half_h * ode(x + i * half_h, y_half, dy_half, homogeneous);
            double k2_half = half_h * (dy_half + 0.5 * l1_half);
            double l2_half = half_h * ode(x + (i + 0.5) * half_h, y_half + 0.5 * k1_half, dy_half + 0.5 * l1_half, homogeneous);
            double k3_half = half_h * (dy_half + 0.5 * l2_half);
            double l3_half = half_h * ode(x + (i + 0.5) * half_h, y_half + 0.5 * k2_half, dy_half + 0.5 * l2_half, homogeneous);
            double k4_half = half_h * (dy_half + l3_half);
            double l4_half = half_h * ode(x + (i + 1) * half_h, y_half + k3_half, dy_half + l3_half, homogeneous);

            y_half += (k1_half + 2 * k2_half + 2 * k3_half + k4_half) / 6.0;
            dy_half += (l1_half + 2 * l2_half + 2 * l3_half + l4_half) / 6.0;
        }

        double error = std::abs(y_half - y_full);

        if (error < tol) {
            result.push_back({ x, y, h });
            x += h;
            y = y_full;
            dy = dy_full;
            h = std::min(h * 1.5, b - x);
            ode1 = ode(x, y, dy, homogeneous);
        }
        else {
            h *= 0.5;
        }
    }
    return result;
}

void runAdaptiveStepStudy(double tol) {
    std::string filename = "solution_adaptive_tol_" + std::to_string(tol) + ".txt";

    auto u_data = rungeKuttaAdaptive(a, u_a, du_a, tol, true);
    auto v_data = rungeKuttaAdaptive(a, v_a, dv_a, tol, false);

    std::vector<Solution> y_data;
    double C = (B - beta0 * v_data.back().y) / (beta0 * u_data.back().y);

    for (size_t i = 0; i < u_data.size(); ++i) {
        double x = u_data[i].x;
        double y_numeric = C * u_data[i].y + v_data[i].y;
        y_data.push_back({ x, y_numeric, u_data[i].step_size });
    }

    writeToFile(filename, y_data, true);
}

void runFixedStepStudy(double h) {
    std::string filename = "solution_fixed_h_" + std::to_string(h) + ".txt";

    auto u_data = rungeKuttaFixed(a, u_a, du_a, h, true);
    auto v_data = rungeKuttaFixed(a, v_a, dv_a, h, false);

    std::vector<Solution> y_data;
    double C = (B - beta0 * v_data.back().y) / (beta0 * u_data.back().y);

    for (size_t i = 0; i < u_data.size(); ++i) {
        double x = u_data[i].x;
        double y_numeric = C * u_data[i].y + v_data[i].y;
        y_data.push_back({ x, y_numeric, h });
    }

    writeToFile(filename, y_data, true);
}

int main() {
    //runFixedStepStudy(0.1);
    //runFixedStepStudy(0.05);

    double tolerances[] = { 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1.0 };
    for (double tol : tolerances) {
        runAdaptiveStepStudy(tol);
    }
    std::cout << "All data written to files for fixed and adaptive steps." << std::endl;

    return 0;
}