Вывод информации о токе катушки

Для более точного соответствия отчёту

.gitignore

@@ -360,4 +360,7 @@ MigrationBackup/
 .ionide/
 
 # Fody - auto-generated XML schema
-FodyWeavers.xsd
+FodyWeavers.xsd
+
+test*.cpp
+*.dat

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 <chrono>
+
+#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,184 @@
+#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;
+
+VectorXd Phi = VectorXd::Zero(DIMENTION);
+VectorXd delta_Phi(DIMENTION);
+VectorXd prev_Phi = VectorXd::Zero(DIMENTION);
+VectorXd prev_prev_Phi = VectorXd::Zero(DIMENTION);
+MatrixXd Y = MatrixXd::Zero(DIMENTION, DIMENTION);
+VectorXd I(DIMENTION);
+
+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() {
+    double d_Id = It / MFt * exp(Phi[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() {
+    double Id = It * exp(Phi[U_Cb] / MFt);
+
+    // Первые три уравнения (производные)
+    I[dU_C1] = Phi[dU_C1] - (Phi[U_C1] - prev_Phi[U_C1]) / dt;
+    I[dU_Cb] = Phi[dU_Cb] - (Phi[U_Cb] - prev_Phi[U_Cb]) / dt;
+    I[dI_L] = Phi[dI_L] - (Phi[I_L] - prev_Phi[I_L]) / dt;
+
+    // Уравнения для переменных состояния
+    I[U_C1] = Phi[U_C1] - Phi[phi_5];
+    I[U_Cb] = Phi[U_Cb] - (Phi[phi_3] - Phi[phi_4]);
+    I[I_L] = L1 * Phi[dI_L] - (Phi[phi_2] - Phi[phi_5]);
+
+    // Уравнения Кирхгофа для узлов
+    I[phi_1] = Phi[I_E1] + (Phi[phi_1] - Phi[phi_2]) / R1;
+    I[phi_2] = -(Phi[phi_1] - Phi[phi_2]) / R1 + Phi[I_L] + (Phi[phi_2] - Phi[phi_3]) / Rb;
+    I[phi_3] = -(Phi[phi_2] - Phi[phi_3]) / Rb + Id + (Phi[phi_3] - Phi[phi_4]) / Ru + Cb * Phi[dU_Cb];
+    I[phi_4] = -Id - (Phi[phi_3] - Phi[phi_4]) / Ru - Cb * Phi[dU_Cb] + Phi[I_E2];
+    I[phi_5] = -Phi[I_L] + C1 * Phi[dU_C1] + Phi[phi_5] / R2;
+
+    // Уравнения для источников
+    I[I_E1] = E1(t_cur) - Phi[phi_1];
+    I[I_E2] = E2 - Phi[phi_4];
+
+    I = -1 * I;
+}
+
+int main() {
+    ofstream result_file(RESULT);
+    double step_count = 0;              // Счетчик числа шагов
+    double it_count = 0;                // Счетчик числа итераций метода Ньютона
+    int counter = 0;                    // Счетчик числа шагов
+    
+    auto start = std::chrono::steady_clock::now();
+    while (t_cur < T) {
+        bool is_end = false;            // Индикатор завершения итераций метода Ньютона
+        size_t n = 0;                   // Счетчик числа итераций метода Ньютона
+        while (!is_end) {
+            // Заполнение матрицы узловых проводимостей и вектора невязок значениями
+            fill_Jacobian();
+            fill_Nevyzka();
+
+            // Решение СЛАУ
+            delta_Phi = Y.partialPivLu().solve(I);
+            Phi += delta_Phi;
+
+            // Проверка точности
+            is_end = delta_Phi.cwiseAbs().maxCoeff() <= ACR;
+            n++;
+            // Если требуемая точность не достигнута
+            if (!is_end) {
+                // Если превышено максимальное число итераций
+                if (n > N_iter) {
+                    n = 0;
+                    dt *= 0.5;
+                    Phi = prev_Phi;
+                    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(((Phi(i) - prev_Phi(i)) / dt - (prev_Phi(i) - prev_prev_Phi(i)) / prev_dt));
+            cur_delta = (tmp > cur_delta) ? tmp : cur_delta;
+        }
+        // Если интегрирование неудовлетворительно по точности
+        if (cur_delta > EPS_MAX && dt > SMN) {
+            // Уменьшаем шаг по времени и отбрасываем результаты текущего шага
+            dt *= 0.5;
+            Phi = prev_Phi;
+        }
+
+        else {
+            // Сохранение значений с предыдущего шага
+            prev_prev_Phi = prev_Phi;
+            prev_Phi = Phi;
+            prev_dt = dt;
+            // Вывод значения переменной на текущем временном шаге
+            result_file
+                << t_cur << ' '
+                << Phi(phi_1) << ' '
+                << Phi(phi_2) << ' '
+                << Phi(phi_3) << ' '
+                << Phi(phi_4) << ' '
+                << Phi(phi_5) << ' '
+                << Phi(I_L) << 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;
+            }
+        }
+    }
+    auto end = std::chrono::steady_clock::now();
+    auto sec = std::chrono::duration<double>(end - start);
+
+    std::cout << "Time: " << sec.count() << " sec\n";
+    cout << "Time steps: " << counter << endl;
+    cout << "Average step by time: " << step_count / counter << endl;
+    cout << "Average iteration count: " << it_count / counter << endl;
+    return 0;
+}

test.ipynb

@@ -0,0 +1,95 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a4654478",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pandas as pd\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "\n",
+    "data = pd.read_csv('result.dat', sep=r'\\s+', header=None)\n",
+    "\n",
+    "n_cols = data.shape[1]\n",
+    "print(f\"Загружено {n_cols} столбцов.\")\n",
+    "\n",
+    "if n_cols < 3:\n",
+    "    raise ValueError(\"Недостаточно столбцов (ожидалось минимум: t, φ, I_L)\")\n",
+    "\n",
+    "time = data[0]\n",
+    "\n",
+    "num_phi = n_cols - 2\n",
+    "print(f\"Обнаружено {num_phi} узлов φ.\")\n",
+    "\n",
+    "phi_cols = list(range(1, 1 + num_phi))\n",
+    "\n",
+    "colors = plt.cm.tab10(np.linspace(0, 1, 10))\n",
+    "\n",
+    "# ---------------------------------------------------------\n",
+    "# Графики φ₁…φ_N\n",
+    "# ---------------------------------------------------------\n",
+    "for idx, col in enumerate(phi_cols, start=1):\n",
+    "    label = f'$\\\\phi_{idx}, В$'\n",
+    "\n",
+    "    fig, ax = plt.subplots(figsize=(8, 4))\n",
+    "\n",
+    "    ax.plot(time, data[col],\n",
+    "            color=colors[idx % len(colors)],\n",
+    "            linewidth=1.8,\n",
+    "            label=label)\n",
+    "\n",
+    "    ax.set_xlabel('$t$, с', fontsize=12)\n",
+    "    ax.set_ylabel(label, fontsize=12)\n",
+    "    ax.grid(True, linestyle='--', linewidth=0.5, alpha=0.7)\n",
+    "\n",
+    "    fig.tight_layout()\n",
+    "    #fig.savefig(f\"program-{idx}.png\", dpi=300)\n",
+    "    plt.show()\n",
+    "\n",
+    "# ---------------------------------------------------------\n",
+    "# График тока I_L\n",
+    "# ---------------------------------------------------------\n",
+    "iL_col = n_cols - 1\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(8, 4))\n",
+    "\n",
+    "ax.plot(time, data[iL_col],\n",
+    "        color='black',\n",
+    "        linewidth=1.8,\n",
+    "        label='$I_L$')\n",
+    "\n",
+    "ax.set_xlabel('$t$, с', fontsize=12)\n",
+    "ax.set_ylabel('$I_L$, A', fontsize=12)\n",
+    "ax.grid(True, linestyle='--', linewidth=0.5, alpha=0.7)\n",
+    "\n",
+    "fig.tight_layout()\n",
+    "#fig.savefig(\"program-IL.png\", dpi=300)\n",
+    "plt.show()\n"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.12.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}