ParkSuMin/MemMAPR-DZ
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Изменение имён переменных

Browse Source 
Для более точного соответствия отчёту
This commit is contained in:
ParkSuMin
2025-11-19 22:50:27 +03:00
parent ed2b7b2d60
commit 8567ba171d
1 changed files with 40 additions and 39 deletions
Download Patch File Download Diff File Expand all files Collapse all files

79
src/main.cpp
View File

@@ -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++;