From 8269990c3e0bdc6029492596a574991582ca27b7 Mon Sep 17 00:00:00 2001 From: ParkSuMin Date: Tue, 28 Oct 2025 14:09:37 +0300 Subject: [PATCH] =?UTF-8?q?=D0=A0=D0=B5=D0=B0=D0=BB=D0=B8=D0=B7=D0=B0?= =?UTF-8?q?=D1=86=D0=B8=D1=8F=20=D0=BC=D0=B5=D1=82=D0=BE=D0=B4=D0=B0=20?= =?UTF-8?q?=D0=9D=D1=8C=D1=8E=D1=82=D0=BE=D0=BD=D0=B0=20=D0=B4=D0=BB=D1=8F?= =?UTF-8?q?=20=D0=B2=D1=8B=D1=87=D0=B8=D1=81=D0=BB=D0=B5=D0=BD=D0=B8=D1=8F?= =?UTF-8?q?=20=D0=B2=D0=B5=D0=BA=D1=82=D0=BE=D1=80=D0=B0=20=D1=81=D0=BC?= =?UTF-8?q?=D0=B5=D1=89=D0=B5=D0=BD=D0=B8=D1=8F?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- main.cpp | 61 ++++++++++++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 57 insertions(+), 4 deletions(-) diff --git a/main.cpp b/main.cpp index 9c897fc..39448b2 100644 --- a/main.cpp +++ b/main.cpp @@ -46,12 +46,14 @@ double E1(double t) { // Якобиан 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; @@ -95,10 +97,61 @@ void nevyzka(double dt, double cur_time) { N(11) = E1(cur_time) - X(6); N(12) = E2 - X(9); + + // Перед вектором невязок стоит минус + N = -1 * N; } int main() { - Jacobian(0.1); + 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; + std::cout << Y << std::endl << N;*/ } \ No newline at end of file