diff --git a/main.cpp b/main.cpp index 0e10110..9c897fc 100644 --- a/main.cpp +++ b/main.cpp @@ -28,7 +28,6 @@ double R2 = 10000; double Ru = 1e+6; double Rb = 20; -double I = 1; double It = 1e-12; double IL_prev = 0; @@ -36,7 +35,7 @@ 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; @@ -49,7 +48,10 @@ double E1(double t) { MatrixXd Y = MatrixXd::Zero(MATRIX_DIM, MATRIX_DIM); // Вектор невязок -VectorXd x = VectorXd::Zero(MATRIX_DIM); +VectorXd N = VectorXd::Zero(MATRIX_DIM); + +// Вектор поправок +VectorXd X = VectorXd::Zero(MATRIX_DIM); void Jacobian(double dt) { Y(0, 0) = 1, Y(0, 3) = -1. / dt; @@ -70,7 +72,33 @@ void Jacobian(double dt) { 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); +} + int main() { Jacobian(0.1); - std::cout << Y; + nevyzka(0.1, 0); + std::cout << Y << std::endl << N; } \ No newline at end of file