From 8567ba171d4e5ba6c02a2aff92aa7f6a3cc98088 Mon Sep 17 00:00:00 2001 From: ParkSuMin Date: Wed, 19 Nov 2025 22:50:27 +0300 Subject: [PATCH] =?UTF-8?q?=D0=98=D0=B7=D0=BC=D0=B5=D0=BD=D0=B5=D0=BD?= =?UTF-8?q?=D0=B8=D0=B5=20=D0=B8=D0=BC=D1=91=D0=BD=20=D0=BF=D0=B5=D1=80?= =?UTF-8?q?=D0=B5=D0=BC=D0=B5=D0=BD=D0=BD=D1=8B=D1=85?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Для более точного соответствия отчёту --- src/main.cpp | 79 ++++++++++++++++++++++++++-------------------------- 1 file changed, 40 insertions(+), 39 deletions(-) diff --git a/src/main.cpp b/src/main.cpp index d3cf078..3353f79 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -9,6 +9,13 @@ 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; @@ -18,8 +25,8 @@ double E1(double t) { } // Заполнение матрицы -void fill_Jacobian(MatrixXd& Y, const VectorXd& cur_X) { - double d_Id = It / MFt * exp(cur_X[U_Cb] / MFt); +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; @@ -65,41 +72,35 @@ void fill_Jacobian(MatrixXd& Y, const VectorXd& cur_X) { 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); +void fill_Nevyzka() { + double Id = It * exp(Phi[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; + 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; // Уравнения для переменных состояния - 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]); + 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]); // Уравнения Кирхгофа для узлов - 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; + 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; // Уравнения для источников - b[I_E1] = E1(t_cur) - cur_X[phi_1]; - b[I_E2] = E2 - cur_X[phi_4]; + I[I_E1] = E1(t_cur) - Phi[phi_1]; + I[I_E2] = E2 - Phi[phi_4]; - b = -1 * b; + I = -1 * I; } 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; // Счетчик числа шагов @@ -110,15 +111,15 @@ int main() { size_t n = 0; // Счетчик числа итераций метода Ньютона while (!is_end) { // Заполнение матрицы узловых проводимостей и вектора невязок значениями - fill_Jacobian(Y, cur_X); - fill_Nevyzka(b, cur_X, prev_X); + fill_Jacobian(); + fill_Nevyzka(); // Решение СЛАУ - delta_X = Y.partialPivLu().solve(b); - cur_X += delta_X; + delta_Phi = Y.partialPivLu().solve(I); + Phi += delta_Phi; // Проверка точности - is_end = delta_X.cwiseAbs().maxCoeff() <= ACR; + is_end = delta_Phi.cwiseAbs().maxCoeff() <= ACR; n++; // Если требуемая точность не достигнута if (!is_end) { @@ -126,7 +127,7 @@ int main() { if (n > N_iter) { n = 0; dt *= 0.5; - cur_X = prev_X; + Phi = prev_Phi; if (dt < SMN) { throw domain_error("Решение не сходится, dt < SMN"); } @@ -137,29 +138,29 @@ int main() { // Расчет текущей ошибки метода Ньютона 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)); + 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; - cur_X = prev_X; + Phi = prev_Phi; } else { // Сохранение значений с предыдущего шага - prev_prev_X = prev_X; - prev_X = cur_X; + prev_prev_Phi = prev_Phi; + prev_Phi = Phi; 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; + << Phi(phi_1) << ' ' + << Phi(phi_2) << ' ' + << Phi(phi_3) << ' ' + << Phi(phi_4) << ' ' + << Phi(phi_5) << std::endl; // Шаг по времени t_cur += dt; counter++;