weno.cpp

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <vector>

template<typename T>
struct wrapped {
    std::vector<T> v;
    wrapped(std::size_t n, const T &e = T()) : v(n, e) { }
    std::size_t size() const { return v.size(); }
    std::ptrdiff_t idx(std::ptrdiff_t i) const {
        const std::size_t s = size();
        if (i < 0)
            return i + s;
        if (i >= s)
            return i - s;
        return i;
    }
    const T &operator[](std::ptrdiff_t i) const { return v[idx(i)]; }
    T &operator[](std::ptrdiff_t i) { return v[idx(i)]; }
};

struct poly {
    double a, b, c;
    poly() { a = b = c = 0; }
    poly(double a, double b, double c) : a(a), b(b), c(c) { }
    double operator()(double x) const {
        return a * x * x + b * x + c;
    }
    poly operator*(double v) const {
        return poly(a * v, b * v, c * v);
    }
    poly operator+(const poly &o) const {
        return poly(a + o.a, b + o.b, c + o.c);
    }
    poly &operator+=(const poly &o) {
        a += o.a;
        b += o.b;
        c += o.c;
        return *this;
    }
    static poly interpolate(double ul, double uc, double ur) {
        double a = 0.5 * (ul - 2 * uc + ur);
        double b = 0.5 * (ur - ul);
        return poly(a, b, uc);
    }
    // Return q(x) = p(x - s)
    poly shift(double s) {
        return poly(a, b - 2 * a * s, (*this)(-s));
    }
    double smoothness() const {
        return 2 * a + a * a / 3 + b * b;
    }
};

struct poly5 {
    double a, b, c, d, e;
    poly5() {
        a = b = c = d = e = 0;
    }
    void add(const poly &g, const poly &p) {
        a += g.a * p.a;
        b += g.a * p.b + g.b * p.a;
        c += g.a * p.c + g.b * p.b + g.c * p.a;
        d += g.b * p.c + g.c * p.b;
        e += g.c * p.c;
    }
    double operator()(double x) const {
        return e + x * (d + x * (c + x * (b + x * a)));
    }
};

struct rat {
    poly5 num;
    poly denom;
    rat(const poly5 &num = poly5(), const poly &denom = poly()) : num(num), denom(denom) { }
    double operator()(const double x) {
        return num(x) / denom(x);
    }
};

template<typename Cont>
rat weno(const Cont &u, std::ptrdiff_t i) {
    const poly pl = poly::interpolate(u[i-2], u[i-1], u[i]).shift(-1);
    const poly pc = poly::interpolate(u[i-1], u[i], u[i+1]);
    const poly pr = poly::interpolate(u[i], u[i+1], u[i+2]).shift(1);

    const poly gl(1. / 12, -1. / 4, 1. / 6);
    const poly gc(-1. / 6, 0, 2. / 3);
    const poly gr(1. / 12, 1. / 4, 1. / 6);

    poly5 num;
    poly  denom;

    const double eps = 1e-6;
    double al = std::pow(pl.smoothness() + eps, -2);
    double ac = std::pow(pc.smoothness() + eps, -2);
    double ar = std::pow(pr.smoothness() + eps, -2);

//    al = ac = ar = 1;

    num.add(gl * al, pl);
    num.add(gc * ac, pc);
    num.add(gr * ar, pr);

    denom += gl * al;
    denom += gc * ac;
    denom += gr * ar;

    return rat(num, denom);
}

double u0(double x) {

    while (x < 0)
        x += 1;
    while (x > 1)
        x -= 1;

    return 0.001e-100 * sin(156 * x) + ((x > .2) && (x < .4) ? 1 : 0);

//    return sin(2 * 3.14159265358929 * x);
}

int main() {
    const size_t N = 100;
    wrapped<double> u(N);
    wrapped<rat> p(u.size());

    const double h = 1. / N;

    for (size_t i = 0; i < N; i++) {
        double x = (i + .5) * h;
        u[i] = u0(x);
    }

    const double C = .15;
    const double dt = C * h;
    double t = 0;
    for (int steps = 0; steps < 20 * N; steps ++) {
        for (size_t i = 0; i < N; i++)
            p[i] = weno(u, i);
        for (size_t i = 0; i < N; i++)
            u[i] += C * (p[i-1](.5) - p[i](.5));
        t += dt;
    }

    for (size_t i = 0; i < u.size(); i++) {
        const int M = 20; // output sample points per interval
        for (int j = 0; j < M; j++) {
            double x = (j + .5) / M;
            std::cout << i + x << " " << p[i](x - .5) << " " << u0((i + x) * h - t) << std::endl;
        }
        std::cout << std::endl;
    }

    return 0;
}