Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #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;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement