Format

Уж больно это сильно мешало при чтении

MemMAPR-MKE.cpp

@@ -6,7 +6,7 @@ int main() {
     Solver slv(5., 4., 0., 1., 30, 0, 10);
     std::cout << "Linear element:";
     slv.Execute_Linear(0, 5);
-    //std::cout << "\nCubic  element:" << std::endl;
+    std::cout << "\nCubic element:";
-    //slv.Execute_Cubic(0, 5);
+    slv.Execute_Cubic(0, 5);
     return 0;
 }

Solver.cpp

@@ -5,8 +5,7 @@ using namespace Eigen;
 Solver::Solver(double _A, double _B, double _C, double _D, int _N, int _l, int _u) {
     A = _A, B = _B, C = _C, D = _D, N = _N;
     upper = _u, lower = _l;
-	L = upper - lower;
+	L = (double)(upper - lower) / N;
-	dx = L / N;
 }
 
 void Solver::Execute_Linear(double val1, double val2) {
@@ -14,13 +13,13 @@ void Solver::Execute_Linear(double val1, double val2) {
     MatrixXd local = MatrixXd::Zero(2, 2);
     VectorXd local_load(2);
 
-    local(0, 0) = -A / dx - B / 2. + C * dx / 3.;
+    local(0, 0) = -A / L - B / 2. + C * L / 3.;
-    local(0, 1) = A / dx + B / 2. + C * dx / 6.;
+    local(0, 1) = A / L + B / 2. + C * L / 6.;
-    local(1, 0) = A / dx - B / 2. + C * dx / 6.;
+    local(1, 0) = A / L - B / 2. + C * L / 6.;
-    local(1, 1) = -A / dx + B / 2. + C * dx / 3.;
+    local(1, 1) = -A / L + B / 2. + C * L / 3.;
 
-    local_load(0) = -D * dx / 2.;
+    local_load(0) = -D * L / 2.;
-    local_load(1) = -D * dx / 2.;
+    local_load(1) = -D * L / 2.;
 
     MatrixXd ansamb = MatrixXd::Zero(N + 1, N + 1);
     VectorXd global_load = VectorXd::Zero(N + 1);
@@ -46,17 +45,13 @@ void Solver::Execute_Linear(double val1, double val2) {
 #endif
 
     double u_right = val2;
-
-    // Clear first and last rows
     ansamb.row(0).setZero();
     ansamb.row(N).setZero();
 
-    // u'(0) = u(0)
+    ansamb(0, 0) = L + 1;
-    ansamb(0, 0) = dx + 1;
     ansamb(0, 1) = -1;
     global_load(0) = 0;
 
-    // u(10) = u_right = 5
     ansamb(N, N) = 1;
     global_load(N) = u_right;
 
@@ -71,7 +66,7 @@ void Solver::Execute_Linear(double val1, double val2) {
     std::cout << solution << std::endl;
 
     std::ofstream file("matrix_linear.txt");
-    for (int i = 0; i < N; i++) {
+    for (int i = 0; i < solution.size(); i++) {
         file << solution(i) << ' ';
     }
     file << std::endl;
@@ -87,22 +82,22 @@ void Solver::Execute_Cubic(double val1, double val2) {
 
     // Assemble matrix
     for (int i = 0; i < mat_dim - 3; i += 3) {
-        Amat(i, i + 0) -= A * 37.0 / 10.0 / dx;
+        Amat(i, i + 0) -= A * 37.0 / 10.0 / L;
-        Amat(i, i + 1) -= A * (-189.0) / 40.0 / dx;
+        Amat(i, i + 1) -= A * (-189.0) / 40.0 / L;
-        Amat(i, i + 2) -= A * 27.0 / 20.0 / dx;
+        Amat(i, i + 2) -= A * 27.0 / 20.0 / L;
-        Amat(i, i + 3) -= A * (-13.0) / 40.0 / dx;
+        Amat(i, i + 3) -= A * (-13.0) / 40.0 / L;
-        Amat(i + 1, i + 0) -= A * (-189.0) / 40.0 / dx;
+        Amat(i + 1, i + 0) -= A * (-189.0) / 40.0 / L;
-        Amat(i + 1, i + 1) -= A * 54.0 / 5.0 / dx;
+        Amat(i + 1, i + 1) -= A * 54.0 / 5.0 / L;
-        Amat(i + 1, i + 2) -= A * (-297.0) / 40.0 / dx;
+        Amat(i + 1, i + 2) -= A * (-297.0) / 40.0 / L;
-        Amat(i + 1, i + 3) -= A * 27.0 / 20.0 / dx;
+        Amat(i + 1, i + 3) -= A * 27.0 / 20.0 / L;
-        Amat(i + 2, i + 0) -= A * 27.0 / 20.0 / dx;
+        Amat(i + 2, i + 0) -= A * 27.0 / 20.0 / L;
-        Amat(i + 2, i + 1) -= A * (-297.0) / 40.0 / dx;
+        Amat(i + 2, i + 1) -= A * (-297.0) / 40.0 / L;
-        Amat(i + 2, i + 2) -= A * 54.0 / 5.0 / dx;
+        Amat(i + 2, i + 2) -= A * 54.0 / 5.0 / L;
-        Amat(i + 2, i + 3) -= A * (-189.0) / 40.0 / dx;
+        Amat(i + 2, i + 3) -= A * (-189.0) / 40.0 / L;
-        Amat(i + 3, i + 0) -= A * (-13.0) / 40.0 / dx;
+        Amat(i + 3, i + 0) -= A * (-13.0) / 40.0 / L;
-        Amat(i + 3, i + 1) -= A * 27.0 / 20.0 / dx;
+        Amat(i + 3, i + 1) -= A * 27.0 / 20.0 / L;
-        Amat(i + 3, i + 2) -= A * (-189.0) / 40.0 / dx;
+        Amat(i + 3, i + 2) -= A * (-189.0) / 40.0 / L;
-        Amat(i + 3, i + 3) -= A * 37.0 / 10.0 / dx;
+        Amat(i + 3, i + 3) -= A * 37.0 / 10.0 / L;
 
         Amat(i + 0, i + 0) += B * (-1.0) / 2.0;
         Amat(i + 0, i + 1) += B * 57.0 / 80.0;
@@ -122,17 +117,16 @@ void Solver::Execute_Cubic(double val1, double val2) {
         Amat(i + 3, i + 3) += B * 1.0 / 2.0;
     }
 
-    // Assembdxe vector
+    // AssembLe vector
     for (int i = 0; i < mat_dim - 3; i += 3) {
-        b(i) -= D * dx / 8.0;
+        b(i) -= D * L / 8.0;
-        b(i + 1) -= D * 3.0 * dx / 8.0;
+        b(i + 1) -= D * 3.0 * L / 8.0;
-        b(i + 2) -= D * 3.0 * dx / 8.0;
+        b(i + 2) -= D * 3.0 * L / 8.0;
-        b(i + 3) -= D * dx / 8.0;
+        b(i + 3) -= D * L / 8.0;
     }
 
-
     Amat.row(0).setZero();
-    Amat(0, 0) = dx / 3.0 + 1;
+    Amat(0, 0) = L / 3.0 + 1;
     Amat(0, 1) = -1;
     b(0) = 0;
 
@@ -141,7 +135,7 @@ void Solver::Execute_Cubic(double val1, double val2) {
     b(mat_dim - 1) = val2;
 
     // Решение системы
-    VectorXd solution = Amat.colPivHouseholderQr().solve(b);
+    VectorXd solution = Amat.fullPivLu().solve(b);
     std::cout << "\nSolution:" << std::endl;
     std::cout << solution << std::endl;
     

Solver.h

@@ -2,7 +2,8 @@
 
 class Solver {
 private:
-	double A, B, C, D, L, dx;
+	double A, B, C, D;
+	double L;
 	int N, upper, lower;
 public:
 	Solver(double _A, double _B, double _C, double _D, int _N, int _l, int _u);