task4.cpp

#pragma GCC optimize("Ofast")
#pragma GCC target("avx")
#include <bits/stdc++.h>
#define watch(x) std::cout << std::setw(12) << #x << " = " << std::setw(12) << x
#define watchln(x) watch(x) << std::endl
typedef double real;
typedef std::vector<real> Vec;
typedef std::vector<Vec> Mat;
const real EPS = real(1e-14);

Mat transpose(const Mat& m);
real dot(const Vec& a, const Vec& b);
Vec operator*(const Vec& a, const Vec& b);
Vec operator-(const Vec& a, const Vec& b);
Vec gauss(Mat A, Vec v);
Vec operator*(const Vec& a, real coeff);
Vec operator*(real coeff, const Vec& a) { return a * coeff; }

class Timer {
    std::chrono::time_point<std::chrono::steady_clock> timePoint;
    size_t value;
public:
    void start() { timePoint = std::chrono::steady_clock::now(); }
    void finish() {
        auto curr = std::chrono::steady_clock::now();
        auto elapsed = std::chrono::duration_cast<std::chrono::milliseconds>(curr - timePoint);
        value = elapsed.count();
    }
    size_t get() const { return value; }
};

Timer timer;

struct Solver {
    std::string title;
    real dx, ds, delta;
    Vec x, s, u, ud, v, z;
    Mat K, a;
    Solver(std::string title_, real coeff_, real xa, real xb, int nx, real sa, real sb, int ns)
        : title(title_), dx((xb - xa) / (nx - 1)), ds((sb - sa) / (ns - 1))
    {
        x = make_grid(xa, xb, nx);
        s = make_grid(sa, sb, ns);
        K.assign(ns, Vec(nx));
        for (int is = 0; is < ns; is++) {
            for (int ix = 0; ix < nx; ix++) {
                K[is][ix] = func_K(x[ix], s[is]);
            }
        }
        std::mt19937 gen(123);
        std::uniform_real_distribution<real> dist(-coeff_, coeff_);
        u.resize(nx);
        ud.resize(nx);
        for (int ix = 0; ix < nx; ix++) {
            u[ix] = func_u(x[ix]);
            ud[ix] = u[ix] + dist(gen);
        }
        v.resize(ns);
        auto integrate = make_integral_coeff(nx);
        delta = std::sqrt(dot((u-ud)*(u-ud), integrate) * dx);
        a.assign(ns, Vec(ns, 0));
        for (int is = 0; is < ns; is++) {
            v[is] = dot(K[is] * u, integrate) * dx;
            for (int it = 0; it < ns; it++) {
                a[is][it] = dot(K[is] * K[it], integrate) * dx * ds;
                //std::cout << a[is][it] << " ";
            }
            //std::cout << std::endl;
        }

        auto denum = std::sqrt(dot(ud*ud, integrate) * dx);
        std::cout << "Solver created, delta is " << delta << ", delta / || ud || = " << delta / denum << std::endl;
        auto norm_u = std::sqrt(dot(u*u, integrate) * dx);
        std::cout << "||ud||=" << denum << std::endl;
        std::cout << "||u||=" << norm_u << std::endl;
        std::cout << "delta/||u||=" << delta / norm_u << std::endl;
    }

    void solve(real alpha) {
        for (int i = 0; i < (int)a.size(); i++) { a[i][i] += alpha; }
        const auto ids2 = 1 / (ds * ds);
        for (int i = 1; i+1 < (int)a.size(); i++) { a[i][i] += 2 * alpha * ids2; };
        a[0][0] += alpha * ids2;
        a.back().back() += alpha * ids2;
        for (int i = 1; i < (int)a.size(); i++) {
            a[i-1][i] -= alpha * ids2;
            a[i][i-1] -= alpha * ids2;
        }
        z = gauss(a, v);
        for (int i = 1; i < (int)a.size(); i++) {
            a[i-1][i] += alpha * ids2;
            a[i][i-1] += alpha * ids2;
        }
        a[0][0] -= alpha * ids2;
        a.back().back() -= alpha * ids2;
        for (int i = 1; i+1 < (int)a.size(); i++) { a[i][i] -= 2 * alpha * ids2; };
        for (int i = 0; i < (int)a.size(); i++) { a[i][i] -= alpha; }
    }

    real diff() {
        const int nx = (int)x.size(), ns = (int)s.size();
        auto int_x = make_integral_coeff(nx);
        auto int_s = make_integral_coeff(ns);
        Vec tmp(nx);
        K = transpose(K);
        for (int ix = 0; ix < nx; ix++) {
            tmp[ix] = dot(K[ix] * z, int_s) * ds - ud[ix];
        }
        K = transpose(K);
        return std::sqrt(dot(tmp * tmp, int_x) * dx) - delta;
    }

    void dump_u() {
        auto file = fopen((title + "_u.txt").c_str(), "wt");
        for (int i = 0; i < (int)x.size(); i++) {
            fprintf(file, "%f %f\n", (double)x[i], (double)u[i]);
        }
        fclose(file);
        std::cout << "File " << (title + "_u.txt") << " is created" << std::endl;
    }

    void dump_ud() {
        auto file = fopen((title + "_ud.txt").c_str(), "wt");
        for (int i = 0; i < (int)x.size(); i++) {
            fprintf(file, "%f %f\n", (double)x[i], (double)ud[i]);
        }
        fclose(file);
        std::cout << "File " << (title + "_ud.txt") << " is created" << std::endl;
    }

    void dump_z() {
        auto file = fopen((title + "_z.txt").c_str(), "wt");
        for (int i = 0; i < (int)s.size(); i++) {
            fprintf(file, "%f %f\n", (double)s[i], (double)z[i]);
        }
        fclose(file);
        std::cout << "File " << (title + "_z.txt") << " is created" << std::endl;
    }

    void dump_ro() {
        auto file = fopen((title + "_ro.txt").c_str(), "wt");
        for (real alpha = 1e-10; alpha < 1; alpha*=10) {
            solve(alpha);
            fprintf(file, "%f %f\n", alpha, diff());
        }
        fclose(file);
        std::cout << "File " << (title + "_ro.txt") << " is created" << std::endl;
    }

    Vec make_grid(real a, real b, int n) {
        Vec res(n);
        real step = (b - a) / (n - 1);
        res[0] = a;
        for (int i = 1; i < n; i++) res[i] = res[i-1] + step;
        assert(std::abs(b-res.back()) < EPS * std::max(real(1), std::abs(b)));
        return res;
    }

    Vec make_integral_coeff(int n) {
        Vec res(n, 1);
        res.front() /= 2;
        res.back()  /= 2;
        return res;
    }

    real func_K(real x, real s) {
        return 1/(1+100*(x-s)*(x-s));
    }

    real func_u(real x) {
        return 0.01 * (10 * (1-x) * (std::atan(10*(1-x)) + std::atan(10*x)) - 0.5 * std::log((1+100*(1-x)*(1-x))/(1+100*x*x)));
    }
};

int main() {
    timer.start();
    Solver solver("task4", 0.02, -1, 1, 501, 0, 1, 501);
    real L = 1e-6, R = 1;
    for (int iter = 0; iter < 100; iter++) {
        real m = (L + R) / 2;
        solver.solve(m);
        if (solver.diff() >= 0) R = m;
        else L = m;
        if (std::abs(solver.diff()) < 0.001 * solver.delta) break;
    }
    real alpha = (L+R) / 2;
    std::cout << "Opt alpha = " << alpha << std::endl;
    solver.solve(alpha);
    std::cout << "Diff is " << solver.diff() << std::endl;
    solver.dump_u();
    solver.dump_ud();
    solver.dump_z();
    solver.dump_ro();
    timer.finish();
    std::cerr << "Done. Runtime is " << timer.get() << " ms" << std::endl;
    return 0;
}

real dot(const Vec& a, const Vec& b) {
    assert(a.size() == b.size());
    return std::inner_product(a.begin(), a.end(), b.begin(), real(0));
}

Vec operator*(const Vec& a, real coeff) {
    Vec res = a;
    for (auto &it : res) it *= coeff;
    return res;
}

Vec operator*(const Vec& a, const Vec& b) {
    assert(a.size() == b.size());
    Vec res(a.size());
    for (int i = 0; i < (int)a.size(); i++) {
        res[i] = a[i] * b[i];
    }
    return res;
}

Vec operator-(const Vec& a, const Vec& b) {
    assert(a.size() == b.size());
    Vec res(a.size());
    for (int i = 0; i < (int)a.size(); i++) {
        res[i] = a[i] - b[i];
    }
    return res;
}

Mat transpose(const Mat& m) {
    Mat res(m[0].size(), Vec(m.size()));
    for (int i = 0; i < (int)res.size(); i++) {
        for (int j = 0; j < (int)res[i].size(); j++) {
            res[i][j] = m[j][i];
        }
    }
    return res;
}

Vec gauss(Mat A, Vec v) {
    int nRows = (int)A.size(), nCols = (int)A[0].size();
    assert(nRows == nCols);
    for (int col = 0, row = 0; col < nCols; col++,row++) {
        int best = row;
        for (int r = row; r < nRows; r++) {
            if (std::abs(A[r][col]) > std::abs(A[best][col])) {
                best = r;
            }
        }
        std::swap(A[row], A[best]);
        std::swap(v[row], v[best]);
        //assert(std::abs(A[row][col]) > 1e-9);
        for (int r = 0; r < nRows; r++) {
            if (r == row) continue;
            const real coeff = -A[r][col] / A[row][col];
            for (int c = col; c < nCols; c++) {
                A[r][c] += A[row][c] * coeff;
            }
            v[r] += v[row] * coeff;
        }
    }
    Vec answ(nCols);
    for (int i = 0; i < nCols; i++) {
        answ[i] = v[i] / A[i][i];
    }
    return answ;
}