кранк-никлсон

gistfile1.cpp

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

using namespace std;

typedef double Scalar;
typedef vector<Scalar> Vector;

struct ConstArr{
    Scalar val;
    Scalar operator[](int i){return val;}
    ConstArr(int val): val(val){}
};

template<typename T>
struct FuncArr{
    vector<bool> used;
    vector<Scalar> arr;
    int l, r, n;
    Scalar step;
    T f;
    Scalar operator[](int i){
        int j = i - l;
        if(j < 0) throw 1;
        return used[j] ? arr[j] : arr[j] = f(i * step);
    }
    FuncArr(T f, int l, int r, Scalar step):
            f(f), l(l), r(r), n(ceil((r - l + 1) / step)), step(step){
        used.resize(n);
        arr.resize(n);
    }
};

template<typename F>
Vector tdma(const int n, const Vector &a, Vector c, const Vector &b, F f){
    Vector x(n);
    Scalar m = 0;
    Vector f1(n);
    f1[0] = f[0];
    for(int i = 1; i < n; i++){
        m = a[i] / c[i - 1];
        c[i] -= m * b[i - 1];
        f1[i] = f[i] - m * f1[i - 1];
    }

   x[n-1] = f1[n-1] / c[n-1];
   for(int i = n - 2; i >= 0; i--)
       x[i] = (f1[i] - b[i] * x[i+1]) / c[i];
   return x;
}

Vector get_tdma_const(int n, int i, Scalar val){
    Vector v(n, val);
    if(i == 0) v[0] = 0;
    else if(i == 1) v[n-1] = 0;
    return v;
}

Vector get_thomas_const(int n, int i, Scalar val){
    Vector v(n, val);
    int a[3] = {0, 1, 0};
    v[0] = v[n-1] = a[i];
    return v;
}

typedef vector<Vector> Solution;

template<class TF, class CU0, class CU1, class CU2>
struct OneDimDiffusionEq{
//u_t = a^2 * u_xx + f(x,t)
//u(x,0) = u0(x)
//u(0, t) = u1(t)
//u(l, t) = u2(t)
//    => u0(0) = u1(0) and u0(l) = u2(0) !!!
    Scalar a;
    Scalar l;
    TF f;
    CU0 u0;
    CU1 u1;
    CU2 u2;
};

template<class TF, class CU0, class CU1, class CU2>
struct OneDimDiffusionEqSolver{
//u_t = a^2 * u_xx + f(x,t)
//u(x,0) = u0(x)
//u(0, t) = u1(t)
//u(l, t) = u2(t)
//    => u0(0) = u1(0) and u0(l) = u2(0) !!!
    Scalar a;
    Scalar l;
    TF f;
    CU0 u0;
    CU1 u1;
    CU2 u2;

    Scalar T, dx, dt;

    struct TdmaRight{
        TF &f;
        Solution &v;
        int n, sz;
        Scalar r, dt;
        CU1 &u1;
        CU2 &u2;
        Scalar operator[](int i){
            if(!i)
                return u1[n];
            else if(i == sz - 1)
                return u2[n];
            return r * v[n][i-1] + 2 * (1 - r) * v[n][i] + r * v[n][i+1] + dt * (f(i, n) + f(i, n + 1));
        }
        TdmaRight(TF &f, Solution &v, Scalar r, Scalar dt, int n, int sz, CU1 &u1, CU2 &u2):
                f(f), v(v), r(r), dt(dt), n(n), sz(sz), u1(u1), u2(u2){}
    };

    Solution sol_crank_nickleson();
    OneDimDiffusionEqSolver(Scalar a, Scalar l, TF f, CU0 u0, CU1 u1, CU2 u2, Scalar T, Scalar dx, Scalar dt)
            :a(a), l(l), f(f), u0(u0), u1(u1), u2(u2), T(T), dx(dx), dt(dt){}
};

template<class F, class A, class B, class C, class D>
OneDimDiffusionEqSolver<F, A, B, C>
get_OneDimDiffusionEqSolver(Scalar a, Scalar l, F f, A u0, B u1, C u2, D T, D dx, D dt){
    return OneDimDiffusionEqSolver<F, A, B, C>(a, l, f, u0, u1, u2, T, dx, dt);
}

template<class TF, class CU0, class CU1, class CU2>
Solution OneDimDiffusionEqSolver<TF, CU0, CU1, CU2>::sol_crank_nickleson(){
    const Scalar r = dt * a * a / (dx * dx);

    const int xn = dx * l, tn = dt * T;
    Vector c[3];
    Scalar vals[3] = {-r, 2 * (1 + r), -r};
    for(int i = 0; i < 3; i++) c[i] = get_thomas_const(xn, i, vals[i]);
    
    Solution v(tn + 1, Vector(xn));

    for(int i = 0; i < xn; i++)
        v[0][i] = u0[i];
    for(int i = 1; i <= tn; i++)
        v[i] = tdma(xn, c[0], c[1], c[2], TdmaRight(f, v, r, dt, i-1, xn, u1, u2));
    return v;
}

Scalar f0(int x, int t){
    return (x == 5) ? 25 : 1;
}

Scalar u0(Scalar x){
    return -x * (x - (1.0/9 + 9));
}

Scalar u1(Scalar t){
    return t;
}

template<Scalar f(Scalar)>
struct Shatai{
    Scalar d;
    Scalar operator()(Scalar x){
        return f(x - d);
    }
    Shatai(Scalar d): d(d){}
};

int main(){
    Scalar tn = 10, xn = 10;
    ConstArr ZERO = 0, ONE = 1, TWO = 25;
    FuncArr<Scalar(*)(Scalar)> U0(u0, 0, xn, 1), U1(u1, 0, tn, 1);
    Solution v = get_OneDimDiffusionEqSolver(1, xn, f0, U0, U1, ONE, tn, 1.0, 1.0).sol_crank_nickleson();
    for(int i = 0; i < v.size(); i++){
        for(int j = 0; j < v[i].size(); j++)
            printf("%.1f\t", v[i][j]);
        cout << endl;
    }
    printf("\nшатать начальные условия\n");
    const int ndxs = 7;
    Scalar dxs[ndxs] = {0.0001, 0.001, 0.1, 0.5, 1.0, 10.0, 100.0};
    printf("ошибка  \tневязка \tневязка/ошибка\n");
    for(int i = 0; i < ndxs; i++){
        FuncArr<Shatai<u0> > SU0(Shatai<u0>(dxs[i]), 0, xn, 1);
        FuncArr<Shatai<u1> > SU1(Shatai<u1>(dxs[i]), 0, tn, 1);

        Solution sv = get_OneDimDiffusionEqSolver(1, xn, f0, SU0, SU1, ONE, tn, 1.0, 1.0).sol_crank_nickleson();

        Scalar m = 0;
        for(int j = 0; j < v.size(); j++)
            for(int k = 0; k < v[i].size(); k++)
                m = max(m, abs(v[j][k] - sv[j][k]));
        printf("%f\t%f\t%f\n", dxs[i], m, m/dxs[i]);
        
        
        
    }
}