Merge pull request #19 from BlackAngel2108/lubaw

chebyshev test

Author: Amazingkivas

Poisson_Equation/headers/Solvers.hpp

@@ -388,8 +388,8 @@ namespace numcpp
         FP Mmin;
         FP Mmax;
     public:
-        ChebyshevIteration(std::vector<FP> initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr<IMatrix> system_matrix, std::vector<FP> b) :
-            ISolver(std::move(initial_approximation), max_iterations, required_precision, std::move(system_matrix), std::move(b)) {}
+        ChebyshevIteration(std::vector<FP> initial_approximation, size_t max_iterations, FP required_precision, std::unique_ptr<IMatrix> system_matrix, std::vector<FP> b, FP min = 0, FP max = 0) :
+            ISolver(std::move(initial_approximation), max_iterations, required_precision, std::move(system_matrix), std::move(b)), Mmin(min), Mmax(max) {}
 
         ChebyshevIteration() = default;
         ~ChebyshevIteration() = default;
@@ -402,28 +402,38 @@ namespace numcpp
         std::vector<FP> solve() const override
         {
             std::vector<FP> approximation = std::move(initial_approximation);
-            double size = (*system_matrix).size();
-            //FP h = sqrt(1.0 / (*system_matrix).at(0, 1));
-            //FP k = sqrt(1.0 / (*system_matrix).at(0, sqrt(size)));
 
             FP k_cheb = 2.0;
             FP tau0 = 1.0 / ((Mmin + Mmax) / 2.0 + (Mmax - Mmin) / 2 * cos(PI / (2.0 * k_cheb) * (1.0 + 2.0 * 0.0)));
             FP tau1 = 1.0 / ((Mmin + Mmax) / 2.0 + (Mmax - Mmin) / 2 * cos(PI / (2.0 * k_cheb) * (1.0 + 2.0 * 1.0)));
 
-
+            std::cout<<"tau1 = "<<tau0<<" tau2 = "<<tau1<<"\n";
             for (size_t i = 0; i < max_iterations; ++i)
             {
                 std::vector<FP> residual = (*system_matrix) * approximation - b;
-
-                FP residual_norm = norm(residual);
-                if (residual_norm <= required_precision) break;
+                std::vector<FP> saved_approximation = approximation;
 
                 if (i % 2 == 0)
-                    approximation = vector_FMA(residual, -tau0, approximation);
+                    approximation = vector_FMA(residual, tau0, approximation);
                 else
-                    approximation = vector_FMA(residual, -tau1, approximation);
+                    approximation = vector_FMA(residual, tau1, approximation);
+                
+                FP approximation_error = error(saved_approximation, approximation);
+                if (approximation_error <= required_precision) break;
             }
-
+            // for (size_t i = 0; i < max_iterations; ++i)
+            // {
+            //     std::vector<FP> residual = (*system_matrix) * approximation - b;
+
+            //     FP residual_norm = norm(residual);
+            //     if (residual_norm <= required_precision) break;
+
+            //     if (i % 2 == 0)
+            //         approximation = vector_FMA(residual, -tau0, approximation);
+            //     else
+            //         approximation = vector_FMA(residual, -tau1, approximation);
+            //     std::cout << residual_norm<< "  "<<  required_precision <<"\n";
+            // }
             return approximation;
         }
     };

Poisson_Equation/src/test.cpp

@@ -118,17 +118,53 @@ void test_task_custom_grid(size_t m, size_t n, std::unique_ptr<numcpp::ISolver>
     auto [duration, solution] = estimate_time(solver);
 }
 
+void test_Chebyshev(size_t n, size_t m) {
+    // size_t n = 2048;
+    // size_t m = 2048;
+    std::cout << "n = " << n << " m = " << m << '\n';
+
+    std::array<double, 4> corners = { -1.0, -1.0, 1.0, 1.0 };
+    std::vector<FP> init_app((n - 1) * (m - 1), 0.0);
+
+    FP h = (corners[2] - corners[0]) / n;
+    FP k = (corners[3] - corners[1]) / m;
+
+    FP Mmin = 4.0 / pow(h, 2) * pow(sin(numcpp::PI / 2.0 / n), 2) + 4.0 / pow(k, 2) * pow(sin(numcpp::PI / 2.0 / m), 2);
+    FP Mmax = 4.0 / pow(h, 2) * pow(sin(numcpp::PI * (n - 1) / 2.0 / n), 2) + 4.0 / pow(k, 2) * pow(sin(numcpp::PI * (m - 1) / 2.0 / m), 2);
+    std::cout << "Mmin = " << Mmin << " Mmax = " << Mmax << '\n';
+    auto u = [](double x, double y) { return exp(1 - pow(x, 2) - pow(y, 2)); };
+    auto f = [](double x, double y) { return -4 * exp(1 - pow(x, 2) - pow(y, 2)) * (pow(y, 2) + pow(x, 2) - 1); };
+
+    auto mu1 = [](double y) { return exp(-pow(y, 2)); };
+    auto mu2 = [](double y) { return exp(-pow(y, 2)); };
+    auto mu3 = [](double x) { return exp(-pow(x, 2)); };
+    auto mu4 = [](double x) { return exp(-pow(x, 2)); };
+
+    numcpp::DirichletProblemSolver<numcpp::Regular> dirichlet_task;
+    dirichlet_task.set_fraction(m, n);
+    dirichlet_task.set_corners(corners);
+    dirichlet_task.set_u(u);
+    dirichlet_task.set_f(f);
+    dirichlet_task.set_boundary_conditions({ mu1, mu2, mu3, mu4 });
+
+    dirichlet_task.set_solver(std::make_unique<numcpp::ChebyshevIteration>(init_app, 1000000000, 0.0000000000001, nullptr, std::vector<FP>(), Mmin, Mmax));
+
+    auto [duration, solution] = estimate_time(dirichlet_task);
+}
+
 int main() 
 {
-    size_t m1 = 1024;
-    size_t n1 = 1024;
+    // size_t m1 = 1024;
+    // size_t n1 = 1024;
+
+    // auto LS_solver1 = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.000000001, nullptr, std::vector<FP>());
+    // test_task_custom_grid(m1, n1, std::move(LS_solver1));
 
-    auto LS_solver1 = std::make_unique<numcpp::ConGrad>(std::vector<FP>(), 1000000, 0.000000001, nullptr, std::vector<FP>());
-    test_task_custom_grid(m1, n1, std::move(LS_solver1));
+    // size_t m2 = 200;
+    // size_t n2 = 200;
 
-    size_t m2 = 200;
-    size_t n2 = 200;
+    // auto LS_solver2 = std::make_unique<numcpp::MinRes>(std::vector<FP>(), 1000000, 0.000001, nullptr, std::vector<FP>());
+    // test_task(m2, n2, std::move(LS_solver2));
 
-    auto LS_solver2 = std::make_unique<numcpp::MinRes>(std::vector<FP>(), 1000000, 0.000001, nullptr, std::vector<FP>());
-    test_task(m2, n2, std::move(LS_solver2));
+    test_Chebyshev(1024,1024);
 }