Ready code for defence

.gitignore

@@ -361,3 +361,5 @@ MigrationBackup/
 
 # Fody - auto-generated XML schema
 FodyWeavers.xsd
+
+test*.cpp

MemMapr_DZ_2.vcxproj

@@ -102,6 +102,7 @@
       <SDLCheck>true</SDLCheck>
       <PreprocessorDefinitions>_DEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <ConformanceMode>true</ConformanceMode>
+      <AdditionalIncludeDirectories>include/</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>
@@ -116,20 +117,22 @@
       <SDLCheck>true</SDLCheck>
       <PreprocessorDefinitions>NDEBUG;_CONSOLE;%(PreprocessorDefinitions)</PreprocessorDefinitions>
       <ConformanceMode>true</ConformanceMode>
+      <AdditionalIncludeDirectories>include/</AdditionalIncludeDirectories>
     </ClCompile>
     <Link>
       <SubSystem>Console</SubSystem>
       <GenerateDebugInformation>true</GenerateDebugInformation>
     </Link>
   </ItemDefinitionGroup>
-  <ItemGroup>
-    <ClCompile Include="main.cpp" />
-  </ItemGroup>
   <ItemGroup>
     <None Include="packages.config" />
   </ItemGroup>
   <ItemGroup>
-    <ClInclude Include="Solver.hpp" />
+    <ClInclude Include="include\Header.hpp" />
+    <ClInclude Include="include\Params.hpp" />
+  </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="src\main.cpp" />
   </ItemGroup>
   <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
   <ImportGroup Label="ExtensionTargets">

MemMapr_DZ_2.vcxproj.filters

@@ -14,17 +14,20 @@
       <Extensions>rc;ico;cur;bmp;dlg;rc2;rct;bin;rgs;gif;jpg;jpeg;jpe;resx;tiff;tif;png;wav;mfcribbon-ms</Extensions>
     </Filter>
   </ItemGroup>
-  <ItemGroup>
-    <ClCompile Include="main.cpp">
-      <Filter>Исходные файлы</Filter>
-    </ClCompile>
-  </ItemGroup>
   <ItemGroup>
     <None Include="packages.config" />
   </ItemGroup>
   <ItemGroup>
-    <ClInclude Include="Solver.hpp">
+    <ClInclude Include="include\Header.hpp">
+      <Filter>Файлы заголовков</Filter>
+    </ClInclude>
+    <ClInclude Include="include\Params.hpp">
       <Filter>Файлы заголовков</Filter>
     </ClInclude>
   </ItemGroup>
+  <ItemGroup>
+    <ClCompile Include="src\main.cpp">
+      <Filter>Исходные файлы</Filter>
+    </ClCompile>
+  </ItemGroup>
 </Project>

Solver.hpp

@@ -1,2 +0,0 @@
-#pragma once
-

include/Header.hpp

@@ -0,0 +1,25 @@
+#include <iostream>
+#include <vector>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+#define M_PI 3.14159265358979323846
+
+enum Variables {
+    dU_C1 = 0,
+    dU_Cb,
+    dI_L,
+    U_C1,
+    U_Cb,
+    I_L,
+    phi_1,
+    phi_2,
+    phi_3,
+    phi_4,
+    phi_5,
+    I_E1,
+    I_E2
+};
+
+constexpr char RESULT[] = "result.dat";

include/Params.hpp

@@ -0,0 +1,22 @@
+#define DIMENTION 13 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+// <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD>
+#define MFt 0.026
+
+#define C1 0.1e-6
+#define Cb 2e-12
+
+#define E2 5
+#define L1 0.1
+#define It 1e-12
+
+#define Rb 20.0
+#define Ru 1000000
+#define R1 1000.0
+#define R2 10000.0
+
+#define T 1e-3 // <20><><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+#define N_iter 7 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+#define ACR 1e-2 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>
+#define EPS_MAX 0.2 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD>
+#define EPS_MIN 1e-8 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><> <20><><EFBFBD><EFBFBD>
+#define SMN 1e-10 // <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD> <20><><EFBFBD> <20><> <20><><EFBFBD><EFBFBD><EFBFBD><EFBFBD><EFBFBD>

main.cpp

@@ -1,157 +0,0 @@
-#include <iostream>
-#include <Eigen/Dense>
-#define M_PI       3.14159265358979323846   // pi
-
-using namespace Eigen;
-
-#define STATE_VARS 3*2					// Сами переменные и их проивзодные
-#define PHI 5							// Количество узлов в схеме
-#define EDS 2							// Количество источиков напряжения
-#define MATRIX_DIM STATE_VARS+PHI+EDS	// Размерность якобиана
-
-const double TK = 1e-3; // время расчета
-const double T_Start = 1e-10; // начальный шаг
-const double SMN = 1e-15; // минимальный шаг
-const double ACR = 1e-2; // точность
-
-// Параметры модели
-double TAU = 5e-9;
-double MFt = 0.026;
-
-double C1 = 0.1;
-double Cb = 2e-12;
-
-double L = 0.001;
-
-double R1 = 1000;
-double R2 = 10000;
-double Ru = 1e+6;
-double Rb = 20;
-
-double It = 1e-12;
-
-double IL_prev = 0;
-double UC1_prev = 0;
-double UCb_prev = 0;
-
-double E2 = 5.;
-// Синусоидальный источник тока
-double E1(double t) {
-	const double amplitude = 10.0;
-	const double period = 1e-4;
-	const double frequency = 2 * M_PI / period;
-	const double phase = 0.0;
-	return amplitude * sin(frequency * t + phase);
-}
-
-// Якобиан
-MatrixXd Y = MatrixXd::Zero(MATRIX_DIM, MATRIX_DIM);
-// Вектор невязок
-VectorXd N = VectorXd::Zero(MATRIX_DIM);
-// Вектор поправок
-VectorXd dX = VectorXd::Zero(MATRIX_DIM);
-// Вектор фазовых переменных
-VectorXd X = VectorXd::Zero(MATRIX_DIM);
-// Копия вектора фазовых переменных (для обработки случая неудачного шага)
-VectorXd X_backup = VectorXd::Zero(MATRIX_DIM);
-
-void Jacobian(double dt) {
-	Y(0, 0) = 1, Y(0, 3) = -1. / dt;
-	Y(1, 1) = 1, Y(1, 4) = -1. / dt;
-	Y(2, 2) = 1, Y(2, 5) = -1. / dt;
-
-	Y(3, 3) = 1., Y(3, 0) = -1.;
-	Y(4, 4) = 1., Y(4, 8) = -1., Y(4, 9) = 1.;
-	Y(5, 2) = L, Y(5, 7) = -1., Y(5, 10) = 1.;
-
-	Y(6, 6) = 1. / R1, Y(6, 7) - 1. / R1, Y(6, 11) = -1.;
-	Y(7, 5) = 1., Y(7, 6) = -1. / R1, Y(7, 7) = 1. / R1 + 1. / Rb, Y(7, 8) = -1. / Rb;
-	Y(8, 1) = Cb, Y(8, 7) = -1. / Rb, Y(8, 8) = 1. / Rb + 1. / Ru, Y(8, 9) = -1. / Ru;
-	Y(9, 1) = -Cb, Y(9, 8) = -1. / Ru, Y(9, 9) = 1. / Ru, Y(9, 11) = 1.;
-	Y(10, 0) = C1, Y(10, 5) = -1., Y(10, 10) = 1. / R2;
-
-	Y(11, 6) = -1.;
-	Y(12, 9) = -1.;
-}
-
-// 0 - 2 производные
-// 3 - 5 переменные состояния
-// 6 - 10 узлы
-// 11 - 12 токи на ЭДС
-
-// Расчёт вектора невязки
-void nevyzka(double dt, double cur_time) {
-	N(0) = X(0) - (X(3) - UC1_prev) / dt;
-	N(1) = X(1) - (X(4) - UCb_prev) / dt;
-	N(2) = X(2) - (X(5) - IL_prev) / dt;
-
-	N(3) = X(3) - X(10);
-	N(4) = X(4) - (X(8) - X(9));
-	N(5) = X(5) * L - (X(7) - X(10));
-	
-	N(6) = -X(11) + (X(6) - X(7)) / R1;
-	N(7) = -(X(6) - X(7)) / R1 + X(5) + (X(7) - X(8)) / Rb;
-	N(8) = -(X(7) - X(8)) / Rb + It + (X(8) - X(9)) / Ru + Cb * X(1);
-	N(9) = -It - (X(8) - X(9)) / Ru - Cb * X(1) + X(11);
-	N(10) = -X(5) + C1 * X(0) + X(10) / R2;
-
-	N(11) = E1(cur_time) - X(6);
-	N(12) = E2 - X(9);
-
-	// Перед вектором невязок стоит минус
-	N = -1 * N;
-}
-
-int main() {
-	double dt = T_Start;			// Шаг по времени
-	double dt_prev = dt;			// Предыдущий шаг по времени
-	double eps_max = ACR;			// Максимальное отклонение
-	double eps_min = 1e-10;			// Минимальное отклонение
-	int counter = 0;				// Счетчик шагов по времени
-	double cur_time = 0;			// Текущее время
-	double eps_cur;
-	double eps;
-	double step_c = 0;
-
-	while (dt <= TK) {
-		// Счётчик итераций
-		int iteration = 0;
-		bool flag = true;
-		while (flag) {
-			Jacobian(dt);
-			nevyzka(dt, cur_time);
-
-			#ifdef _DEBUG
-			std::cout << Y << std::endl << N;
-			#endif
-			
-			dX = Y.fullPivLu().solve(N);
-			#ifdef _DEBUG
-			std::cout << '\n' << dX;
-			#endif
-
-			X += dX;
-			iteration++;
-			double maximum = X.maxCoeff();
-			if (maximum <= ACR) {
-				flag = false;
-			}
-			else if (iteration == 7 && flag) {
-				iteration = 0;
-				dt /= 2.;
-				if (dt < SMN) {
-					std::cerr << "Sorry, but i can not solve this" << std::endl;
-					exit(-1);
-				}
-				X = X_backup;
-			}
-		}
-		// TODO: реализовать вычисление уточнённое решения, сравнить его с вычисленным, 
-		// возможно изменить шаг моделирования
-		
-	}
-	// Для примера
-	/*Jacobian(0.1);
-	nevyzka(0.1, 0);
-	std::cout << Y << std::endl << N;*/
-}

src/main.cpp

@@ -0,0 +1,176 @@
+#include <Eigen/Dense>
+#include "Header.hpp"
+#include "Params.hpp"
+
+double dt = 1e-5; // Начальный шаг по времени
+double prev_dt = dt;
+double t_cur = 0;
+
+using namespace std;
+using namespace Eigen;
+
+double E1(double t) {
+    const double amplitude = 10.0;
+    const double period = 1e-4;
+    const double frequency = 2 * M_PI / period;
+    const double phase = 0.0;
+    return amplitude * sin(frequency * t + phase);
+}
+
+// Заполнение матрицы
+void fill_Jacobian(MatrixXd& Y, const VectorXd& cur_X) {
+    double d_Id = It / MFt * exp(cur_X[U_Cb] / MFt);
+
+    for (int i = 0; i < 3; i++) {
+        Y(i, i) = 1.0;
+        Y(i, i + 3) = -1.0 / dt;
+    }
+
+    Y(U_C1, U_C1) = 1.;
+    Y(U_C1, phi_5) = -1.;
+
+    Y(U_Cb, U_Cb) = 1.;
+    Y(U_Cb, phi_3) = -1.;
+    Y(U_Cb, phi_4) = 1.;
+
+    Y(I_L, dI_L) = L1;
+    Y(I_L, phi_2) = -1.;
+    Y(I_L, phi_5) = 1.;
+
+    Y(phi_1, phi_1) = 1. / R1;
+    Y(phi_1, phi_2) = -1. / R1;
+    Y(phi_1, I_E1) = 1.;
+
+    Y(phi_2, I_L) = 1.;
+    Y(phi_2, phi_1) = -1. / R1;
+    Y(phi_2, phi_2) = 1. / R1 + 1. / Rb;
+    Y(phi_2, phi_3) = -1. / Rb;
+
+    Y(phi_3, dU_Cb) = Cb;
+    Y(phi_3, phi_2) = -1. / Rb;
+    Y(phi_3, phi_3) = 1. / Rb + 1. / Ru + d_Id;
+    Y(phi_3, phi_4) = -1. / Ru - d_Id;
+
+    Y(phi_4, dU_Cb) = -Cb;
+    Y(phi_4, phi_3) = -1. / Ru - d_Id;
+    Y(phi_4, phi_4) = 1. / Ru + d_Id;
+    Y(phi_4, I_E2) = 1.;
+
+    Y(phi_5, dU_C1) = C1;
+    Y(phi_5, I_L) = -1.;
+    Y(phi_5, phi_5) = 1. / R2;
+
+    Y(I_E1, phi_1) = -1.;
+
+    Y(I_E2, phi_4) = -1.;
+}
+// Заполнение вектора невязок
+void fill_Nevyzka(VectorXd& b, const VectorXd& cur_X, const VectorXd& prev_X) {
+    double Id = It * exp(cur_X[U_Cb] / MFt);
+
+    // Первые три уравнения (производные)
+    b[dU_C1] = cur_X[dU_C1] - (cur_X[U_C1] - prev_X[U_C1]) / dt;
+    b[dU_Cb] = cur_X[dU_Cb] - (cur_X[U_Cb] - prev_X[U_Cb]) / dt;
+    b[dI_L] = cur_X[dI_L] - (cur_X[I_L] - prev_X[I_L]) / dt;
+
+    // Уравнения для переменных состояния
+    b[U_C1] = cur_X[U_C1] - cur_X[phi_5];
+    b[U_Cb] = cur_X[U_Cb] - (cur_X[phi_3] - cur_X[phi_4]);
+    b[I_L] = L1 * cur_X[dI_L] - (cur_X[phi_2] - cur_X[phi_5]);
+
+    // Уравнения Кирхгофа для узлов
+    b[phi_1] = cur_X[I_E1] + (cur_X[phi_1] - cur_X[phi_2]) / R1;
+    b[phi_2] = -(cur_X[phi_1] - cur_X[phi_2]) / R1 + cur_X[I_L] + (cur_X[phi_2] - cur_X[phi_3]) / Rb;
+    b[phi_3] = -(cur_X[phi_2] - cur_X[phi_3]) / Rb + Id + (cur_X[phi_3] - cur_X[phi_4]) / Ru + Cb * cur_X[dU_Cb];
+    b[phi_4] = -Id - (cur_X[phi_3] - cur_X[phi_4]) / Ru - Cb * cur_X[dU_Cb] + cur_X[I_E2];
+    b[phi_5] = -cur_X[I_L] + C1 * cur_X[dU_C1] + cur_X[phi_5] / R2;
+
+    // Уравнения для источников
+    b[I_E1] = E1(t_cur) - cur_X[phi_1];
+    b[I_E2] = E2 - cur_X[phi_4];
+
+    b = -1 * b;
+}
+
+int main() {
+    ofstream result_file(RESULT);
+    VectorXd delta_X(DIMENTION);
+    VectorXd cur_X = VectorXd::Zero(DIMENTION);
+    VectorXd prev_X = VectorXd::Zero(DIMENTION);
+    VectorXd prev_prev_X = VectorXd::Zero(DIMENTION);
+    MatrixXd Y = MatrixXd::Zero(DIMENTION, DIMENTION);
+    VectorXd b(DIMENTION);
+    double step_count = 0;              // Счетчик числа шагов
+    double it_count = 0;                // Счетчик числа итераций метода Ньютона
+    int counter = 0;                    // Счетчик числа шагов
+    while (t_cur < T) {
+        bool is_end = false;            // Индикатор завершения итераций метода Ньютона
+        size_t n = 0;                   // Счетчик числа итераций метода Ньютона
+        while (!is_end) {
+            // Заполнение матрицы узловых проводимостей и вектора невязок значениями
+            fill_Jacobian(Y, cur_X);
+            fill_Nevyzka(b, cur_X, prev_X);
+
+            // Решение СЛАУ
+            delta_X = Y.partialPivLu().solve(b);
+            cur_X += delta_X;
+
+            // Проверка точности
+            is_end = delta_X.cwiseAbs().maxCoeff() <= ACR;
+            n++;
+            // Если требуемая точность не достигнута
+            if (!is_end) {
+                // Если превышено максимальное число итераций
+                if (n > N_iter) {
+                    n = 0;
+                    dt *= 0.5;
+                    cur_X = prev_X;
+                    if (dt < SMN) {
+                        throw domain_error("Решение не сходится, dt < SMN");
+                    }
+                }
+            }
+        }
+        double cur_delta = 0.0;
+        // Расчет текущей ошибки метода Ньютона
+        for (int i = dU_C1; i < DIMENTION; i++) {
+            //if (i == dU_Cb) continue;
+            double tmp = 0.5 * dt * dt * abs(((cur_X(i) - prev_X(i)) / dt - (prev_X(i) - prev_prev_X(i)) / prev_dt));
+            cur_delta = (tmp > cur_delta) ? tmp : cur_delta;
+        }
+        // Если интегрирование неудовлетворительно по точности
+        if (cur_delta > EPS_MAX && dt > SMN) {
+            // Уменьшаем шаг по времени и отбрасываем результаты текущего шага
+            dt *= 0.5;
+            cur_X = prev_X;
+        }
+
+        else {
+            // Сохранение значений с предыдущего шага
+            prev_prev_X = prev_X;
+            prev_X = cur_X;
+            prev_dt = dt;
+            // Вывод значения переменной на текущем временном шаге
+            result_file
+                << t_cur << ' '
+                << cur_X(phi_1) << ' '
+                << cur_X(phi_2) << ' '
+                << cur_X(phi_3) << ' '
+                << cur_X(phi_4) << ' '
+                << cur_X(phi_5) << std::endl;
+            // Шаг по времени
+            t_cur += dt;
+            counter++;
+            step_count += dt;
+            it_count += n;
+            cout << counter << " n_it=" << n << " t=" << t_cur << endl;
+            if (cur_delta < EPS_MIN) {
+                dt *= 2;
+            }
+        }
+    }
+    cout << "Число шагов по времени: " << counter << endl;
+    cout << "Средний шаг по времени: " << step_count / counter << endl;
+    cout << "Среднее число итераций: " << it_count / counter << endl;
+    return 0;
+}