МКЭ

#include <stdio.h>
#include <iostream>
#include <vector>

using namespace std;

class finite_elements {
private:
    vector<double> A; // Диагональ слева
    vector<double> B; // Диагональ справа
    vector<double> C; // Главная диагональ
    vector<double> F; // Вектор правой части

    int n;  // Количество элементов
    int n1; // Количество узлов

public:
    finite_elements(int m) {
        n = m;
        n1 = m + 1;
        A = B = C = F = vector<double>(n1, 0);

        filling();
        progonka();
    }

    void progonka() {
        vector<double> X(n1);
        vector<double> alpha(n1 + 1, 0);
        vector<double> beta(n1 + 1, 0);
        alpha[1] = -B[0] / C[0];
        beta[1] = F[0] / C[0];
        for (int i = 1; i < (A.size() - 1); i++) {
            alpha[i + 1] = -B[i] / (A[i] * alpha[i] + C[i]);
            beta[i + 1] = (F[i] - A[i] * beta[i]) / (A[i] * alpha[i] + C[i]);
        }
        int n = A.size() - 1;
        X[n] = (F[n] - A[n] * beta[n]) / (C[n] + A[n] * alpha[n]);
        for (int i = (A.size() - 1); i > 0; i--)
            X[i - 1] = alpha[i] * X[i] + beta[i];

        for (int i = 0; i <= n; i++)
            cout << X[i] << endl;
        system("pause");
    }

    double get_k(double x) {
        if ((x > 0.3) && (x < 0.6))
            return 0.001;
        else return 1;
        //return 1;
    }

    /**
     * Значение функции в точке X
     */
    double func(double x) {
        return 0;
    }

    double inc_integral(double left, double right) {
        double sum = 0;
        double h = (right - left) / 100.f;
        double x = left - 0.5*h;
        for (int i = 0; i < 100; i++) {
            x += h; // значение х
            sum += (func(x) * n * x) * h;
        }
        return sum;
    }

    double dec_integral(double left, double right) {
        double sum = 0;
        double h = (right - left) / 100.f;
        double x = left - 0.5*h;
        for (int i = 0; i < 100; i++) {
            x += h; // значение х
            sum += (func(x) * n * (right - x)) * h;
        }
        return sum;
    }

    void filling() {
        for (int i = 0; i < n; i++) {
            F[i]     += dec_integral((double)i / n, (double)(i + 1) / n);
            F[i + 1] += inc_integral((double)i / n, (double)(i + 1) / n);
        }
        for (int i = 0; i < n; i++) {
            C[i] += n * get_k((i + 0.5) / n);
            B[i] = -1 * n * get_k((i + 0.5) / n);
            A[i + 1] = -1 * n * get_k((i + 0.5) / n);
            C[i + 1] += n * get_k((i + 0.5) / n);
        }
        C[0] = 1;
        B[0] = 0;
        F[0] = 20;
        A[n] = 0;
        C[n] = 1;
        F[n] = 50;
    }
};

int main() {
    finite_elements A(10);
    return 0;
}