Решение волнового уравнения методом конечных разностей

#include "basic_libs.h"
Matrix wave_equation_solver( size_t n, double L, double t0, double t_end, double tau,
                            auto tension, auto density, auto force,auto boundary_condition_left,
                            auto boundary_condition_right,Matrix& u0, const std::string& filename)
{
    double dx = L / (n - 1);
 
    double max_speed = 1.0;
    for (double x = 0.0; x <= L; x += dx) {
        double T = tension(x);
        double rho = density(x);
        if (T <= 0 || rho <= 0) {
            throw std::runtime_error("Tension and density must be positive.");
        }
        max_speed = std::max(max_speed, std::sqrt(T / rho));
    }
    if (tau > dx / max_speed) {
        throw std::runtime_error("Time step is too large (CFL condition violated).");
    }
 
    auto wave_equation_rhs = [&](Matrix u, double t) {
        Matrix dudt(2 * n, Vector(1, 0.0));
 
        for (size_t i = 0; i < n; ++i) {
            dudt[i][0] = u[n + i][0];
        }
 
        for (size_t i = 1; i < n - 1; ++i) {
            double x = i * dx;
            double T_left  = tension(x - dx / 2);
            double T_right = tension(x + dx / 2);
            double rho = density(x);
 
            double d2u_dx2 = (T_right * (u[i + 1][0] - u[i][0]) -
                              T_left  * (u[i][0] - u[i - 1][0])) / (dx * dx);
 
            dudt[n + i][0] = d2u_dx2 / rho + force(x, t) / rho;
        }
 
        dudt[n + 0][0] = boundary_condition_left(t);
        dudt[n + n - 1][0] = (boundary_condition_right(t)-u[n-1][0])/(dx*dx);
 
        return dudt;
    };
 
    return Ralston(wave_equation_rhs, u0, t0, t_end, tau, filename);
}
 
void test_wave()
{
    size_t n = 101;
    double L = 1.0;
    double t0 = 0.0;
    double t_end = 1.0;
    double tau = 0.0001;
 
    auto tension = [](double x) {
        return 1.0;
    };
 
    auto density = [](double x) {
        return 1.0;
    };
 
    auto force = [](double x, double t) {
        return 0.0;
    };
 
    auto boundary_condition_left = [](double t) {
        return 0.0;
    };
 
    auto boundary_condition_right = [](double t) {
        return 0.0;
    };
    auto initial_displacement = [L](double x) {
        return cos(1*M_PI * x/(2*L))+cos(3*M_PI * x/(2*L))+cos(5*M_PI * x/(2*L))+cos(7*M_PI * x/(2*L))+cos(9*M_PI * x/(2*L))+cos(11*M_PI * x/(2*L))+cos(13*M_PI * x/(2*L))+cos(15*M_PI * x/(2*L));
    };
    auto initial_velocity = [L](double x) {
        return 0.0;
    };
 
    Matrix u0(2*n, Vector(1, 0.0));
    for (size_t i = 0; i < n; ++i) {
        double x = i * L / (n - 1);
        u0[i][0] = initial_displacement(x);
        u0[n+i][0] = initial_velocity(x);
    }
 
    std::string filename = "wave_solution_test.txt";
    Matrix result = wave_equation_solver(n, L, t0, t_end, tau, tension, density, force,
                                         boundary_condition_left, boundary_condition_right, u0, filename);
 
    std::cout << "Solution saved to " << filename << std::endl;
}Matrix Ralston(const Function f, const Matrix& u0, double t0, double t_end, double tau, const std::string& filename) {
    Matrix u = u0;
    double t = t0;
    Matrix result;
    double t_between_write = 0.01;
 
    size_t n = u0.size();
 
    std::vector<std::pair<double, Matrix>> data_points;
    data_points.push_back({t, u});
 
    do {
        Matrix k1 = f(u, t);
        Matrix u1 = u;
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < u[i].size(); ++j) {
                u1[i][j] = u[i][j] + tau * 0.4 * k1[i][j];
            }
        }
 
        Matrix k2 = f(u1, t + 0.4 * tau);
        Matrix u2 = u;
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < u[i].size(); ++j) {
                u2[i][j] = u[i][j] + tau * (0.29697761 * k1[i][j] + 0.15875964 * k2[i][j]);
            }
        }
 
        Matrix k3 = f(u2, t + 0.45573725 * tau);
        Matrix u3 = u;
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < u[i].size(); ++j) {
                u3[i][j] = u[i][j] + tau * (0.21810040 * k1[i][j] - 3.05096516 * k2[i][j] + 3.83286476 * k3[i][j]);
            }
        }
 
        Matrix k4 = f(u3, t + tau);
        Matrix u_new = u;
        for (size_t i = 0; i < n; ++i) {
            for (size_t j = 0; j < u[i].size(); ++j) {
                u_new[i][j] = u[i][j] + tau * (0.17476028 * k1[i][j] - 0.55148066 * k2[i][j] + 1.20553560 * k3[i][j] + 0.17118478 * k4[i][j]);
            }
        }
 
        u = u_new;
        if (t < t_end && t + tau > t_end) {
            tau = t_end - t;
        }
        t += tau;
        data_points.push_back({t, u});
 
    } while (std::abs(t - t_end) > 1e-10);
 
    std::ofstream outputFile(filename);
    if (!outputFile.is_open()) {
        throw std::runtime_error("Error opening file: " + filename);
    }
 
    outputFile << std::fixed << std::setprecision(10);
    outputFile << "# Time and values\n";
 
    double write_time = t0;
    while (write_time <= t_end) {
        size_t i = 0;
        while (i + 1 < data_points.size() && data_points[i + 1].first < write_time) {
            ++i;
        }
 
        double t1 = data_points[i].first;
        double t2 = data_points[i + 1].first;
        Matrix u1 = data_points[i].second;
        Matrix u2 = data_points[i + 1].second;
 
        double alpha = (t2 - t1) > 1e-10 ? (write_time - t1) / (t2 - t1) : 0.0;
        Matrix interpolated_u(u1.size(), Vector(u1[0].size(), 0.0));
        for (size_t j = 0; j < u1.size(); ++j) {
            for (size_t k = 0; k < u1[j].size(); ++k) {
                interpolated_u[j][k] = u1[j][k] + alpha * (u2[j][k] - u1[j][k]);
            }
        }
 
        outputFile << write_time;
        for (const auto& row : interpolated_u) {
            for (double val : row) {
                outputFile << " " << val;
            }
        }
        outputFile << std::endl;
 
        result.push_back({write_time, interpolated_u[0][0], interpolated_u[1][0]});
        write_time += t_between_write;
    }
 
    if (std::abs(result.back()[0] - t_end) > 1e-10) {
        result.push_back({t_end, u[0][0], u[1][0]});
    }
 
    if (std::abs(data_points.back().first - t_end) <= 1e-10) {
        outputFile << t_end;
        for (const auto& row : data_points.back().second) {
            for (double val : row) {
                outputFile << " " << val;
            }
        }
        outputFile << std::endl;
    }
 
    return result;
}