API tools faq
paste
BorrowTheProgrammer

trud_1

BorrowTheProgrammer
Dec 4th, 2022
462
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 15.05 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <cmath>
  3. #include <vector>
  4. #include <fstream>
  5. using namespace std;
  6.  
  7. #define DT_START 1e-7
  8. #define DT_MIN 1e-8
  9. #define INTEGRATION_TIME 1e-3
  10. #define epsilon 1e-3 // for delta X
  11. #define A 10.0
  12. #define Freq 1e-4
  13. #define RB 20.0
  14. #define MFt 0.026
  15. #define Cb 2e-12
  16. #define It 1e-12
  17. #define Ru 1e6
  18. #define C3 1e-6
  19. #define L 1e-3
  20. #define DELTA1 1e-8
  21. #define DELTA2 1e-3
  22. #define INITIAL_APPROXIMATION 0.0
  23. #define MAX_NEWTON_ITER 7
  24. #define _USE_MATH_DEFINES
  25.  
  26.  
  27. /**
  28.  * @brief Функция для решения матричного уравнения вида Ax = b
  29.  *
  30.  * @param Y vector< vector<double>& Y
  31.  * @param I vector<double>& I
  32.  * @param n int
  33.  * @return vector<double>
  34.  */
  35. vector<double> SolveMatrixEquation (vector< vector<double>>& Y, vector<double>& I, int n)
  36. {
  37.     double MaxElement;
  38.     vector<double> x(n, 0.0);
  39.     int k, Index;
  40.     const double eps = 1e-12;
  41.     k = 0;
  42.  
  43.     while (k < n)
  44.     {
  45.         MaxElement = abs(Y[k][k]);
  46.         Index = k;
  47.  
  48.         for (int i = k + 1; i < n; i++)
  49.         {
  50.             if (abs(Y[i][k]) > MaxElement)
  51.             {
  52.                 MaxElement = abs(Y[i][k]);
  53.                 Index = i;
  54.             }
  55.         }
  56.  
  57.         if (MaxElement < eps)
  58.         {
  59.             std::cout << "Максимальный элемент нулевой" << endl;
  60.             exit(-1);
  61.         }
  62.  
  63.         for (int j = 0; j < n; j++)
  64.         {
  65.             double Temp = Y[k][j];
  66.             Y[k][j] = Y[Index][j];
  67.             Y[Index][j] = Temp;
  68.         }
  69.  
  70.         double Temp = I[k];
  71.         I[k] = I[Index];
  72.         I[Index] = Temp;
  73.  
  74.         for (int i = k; i < n; i++)
  75.         {
  76.             double Temp = Y[i][k];
  77.  
  78.             if (abs(Temp) < eps) {
  79.                 continue;
  80.             };
  81.  
  82.             for (int j = 0; j < n; j++)
  83.             {
  84.                 Y[i][j] = Y[i][j] / Temp;
  85.             }
  86.                
  87.             I[i] = I[i] / Temp;
  88.  
  89.             if (i == k) {
  90.                 continue;
  91.             }
  92.  
  93.             for (int j = 0; j < n; j++)
  94.             {
  95.                 Y[i][j] = Y[i][j] - Y[k][j];
  96.             }
  97.  
  98.             I[i] = I[i] - I[k];
  99.         }
  100.  
  101.         k++;
  102.     }
  103.  
  104.     for (k = n - 1; k >= 0; k--)
  105.     {
  106.         x[k] = I[k];
  107.  
  108.         for (int i = 0; i < k; i++)
  109.             I[i] = I[i] - Y[i][k] * x[k];
  110.     }
  111.    
  112.     return x;
  113. }
  114.  
  115. /**
  116.  * @brief Функция для формирования файла с данными для построения графиков
  117.  *
  118.  * @param result
  119.  * @param step
  120.  * @param n
  121.  */
  122. void MakePlotData (double* result, int step, int n)
  123. {
  124.     fstream file;
  125.     file.open("dz.txt", fstream::out);
  126.  
  127.     for(int s = 0; s < step; s++)
  128.     {
  129.         file<<"time = "<<result[s*(n+1)]<<endl;
  130.         file << "d_phi1 = " << result[s * (n + 1) + 1] << endl;
  131.         file << "d_phi2 = " << result[s * (n + 1) + 2] << endl;
  132.         file << "d_phi3 = " << result[s * (n + 1) + 3] << endl;
  133.         file << "d_phi4 = " << result[s * (n + 1) + 4] << endl;
  134.         file << "d_phi5 = " << result[s * (n + 1) + 5] << endl;
  135.         file << "int_phi1 = " << result[s * (n + 1) + 6] << endl;
  136.         file << "int_phi2 = " << result[s * (n + 1) + 7] << endl;
  137.         file << "int_phi3 = " << result[s * (n + 1) + 8] << endl;
  138.         file << "int_phi4 = " << result[s * (n + 1) + 9] << endl;
  139.         file << "int_phi5 = " << result[s * (n + 1) + 10] << endl;
  140.         file << "phi1 = " << result[s * (n + 1) + 11] << endl;
  141.         file << "phi2 = " << result[s * (n + 1) + 12] << endl;
  142.         file << "phi3 = " << result[s * (n + 1) + 13] << endl;
  143.         file << "phi4 = " << result[s * (n + 1) + 14] << endl;
  144.         file << "phi5 = " << result[s * (n + 1) + 15] << endl;
  145.         file << "IE1 = " << result[s * (n + 1) + 16] << endl;
  146.         file << "IE3 = " << result[s * (n + 1) + 17] << endl;
  147.     }
  148.  
  149.     file.close();
  150. }
  151.  
  152. int FindDeltaNorm (vector<double> DeltaX, int n)
  153. {
  154.     int count_of_x = 0;
  155.  
  156.     for (int ni = 0; ni < n; ni++)
  157.     {            
  158.         if (DeltaX[ni] < epsilon)
  159.         {
  160.             count_of_x++;
  161.         }
  162.     }
  163.  
  164.     return count_of_x;
  165. }
  166.  
  167. double CalculateEPS (double* res, int step, double delta_t, double delta_t_pred, int n)
  168. {
  169.     double eps = 0.0;
  170.     double dxdt = 0.0;
  171.     for(int i=1; i<n+1; i++)
  172.     {
  173.         dxdt = fabs(res[step*(n+1)+i]-res[(step-1)*(n+1)+i]*delta_t/delta_t_pred);
  174.         if (eps < dxdt)
  175.             eps = dxdt;
  176.     }
  177.     return eps*delta_t*delta_t/2;
  178. }
  179.  
  180. int main()
  181. {
  182.     int n = 17;
  183.     vector< vector<double> > Y(n, vector<double>(n));
  184.     vector<double> I(n, 0.0);
  185.     vector<double> DeltaX(n, 0.0);
  186.     unsigned int totalSteps = 0;
  187.  
  188.     double *Result;
  189.     Result = new double[2100000];
  190.    
  191.     int NewtonIter = 0;
  192.     double t = 0.0;
  193.  
  194.     double  phi1_past = INITIAL_APPROXIMATION;
  195.     double  phi2_past = INITIAL_APPROXIMATION;
  196.     double  phi3_past = INITIAL_APPROXIMATION;
  197.     double  phi4_past = INITIAL_APPROXIMATION;
  198.     double  phi5_past = INITIAL_APPROXIMATION;
  199.     double  int_phi1_past = INITIAL_APPROXIMATION;
  200.     double  int_phi2_past = INITIAL_APPROXIMATION;
  201.     double  int_phi3_past = INITIAL_APPROXIMATION;
  202.     double  int_phi4_past = INITIAL_APPROXIMATION;
  203.     double  int_phi5_past = INITIAL_APPROXIMATION;
  204.  
  205.     double d_phi1 = INITIAL_APPROXIMATION;
  206.     double d_phi2 = INITIAL_APPROXIMATION;
  207.     double d_phi3 = INITIAL_APPROXIMATION;
  208.     double d_phi4 = INITIAL_APPROXIMATION;
  209.     double d_phi5 = INITIAL_APPROXIMATION;
  210.     double int_phi1 = INITIAL_APPROXIMATION;
  211.     double int_phi2 = INITIAL_APPROXIMATION;
  212.     double int_phi3 = INITIAL_APPROXIMATION;
  213.     double int_phi4 = INITIAL_APPROXIMATION;
  214.     double int_phi5 = INITIAL_APPROXIMATION;
  215.     double phi1 = INITIAL_APPROXIMATION;
  216.     double phi2 = INITIAL_APPROXIMATION;
  217.     double phi3 = INITIAL_APPROXIMATION;
  218.     double phi4 = INITIAL_APPROXIMATION;
  219.     double phi5 = INITIAL_APPROXIMATION;
  220.     double IE1 = INITIAL_APPROXIMATION;
  221.     double IE3 = INITIAL_APPROXIMATION;
  222.  
  223.     Result[0] = 0.0;
  224.  
  225.     for (int i = 1; i < n + 1; i++)
  226.     {
  227.         Result[i] = INITIAL_APPROXIMATION;
  228.     }
  229.        
  230.     double dt = DT_START;
  231.     double dt_prev = 0.0;
  232.     double eps = 0.0;
  233.     int Flag = 1;
  234.     int Counter = 0;
  235.     int Step = 1;
  236.  
  237.     while (t < INTEGRATION_TIME)
  238.     {
  239.         // totalSteps++;
  240.  
  241.         while (Flag)
  242.         {    
  243.             for (auto& Row : Y)
  244.             {
  245.                 fill(Row.begin(), Row.end(), 0.0);
  246.             }
  247.  
  248.             for (auto& it : I)
  249.             {
  250.                 it = 0.0;
  251.             }
  252.  
  253.             NewtonIter++;
  254.  
  255.             Y[0][0] = 1.0;
  256.             Y[0][10] = -1.0 / dt;
  257.  
  258.             Y[1][1] = 1.0;
  259.             Y[1][11] = -1.0 / dt;
  260.  
  261.             Y[2][2] = 1.0;
  262.             Y[2][12] = -1.0 / dt;
  263.  
  264.             Y[3][3] = 1.0;
  265.             Y[3][13] = -1.0 / dt;
  266.  
  267.             Y[4][4] = 1.0;
  268.             Y[4][14] = -1.0 / dt;
  269.  
  270.             Y[5][5] = 1.0;
  271.             Y[5][10] = -dt;
  272.  
  273.             Y[6][6] = 1.0;
  274.             Y[6][11] = -dt;
  275.  
  276.             Y[7][7] = 1.0;
  277.             Y[7][12] = -dt;
  278.  
  279.             Y[8][8] = 1.0;
  280.             Y[8][13] = -dt;
  281.  
  282.             Y[9][9] = 1.0;
  283.             Y[9][14] = -dt;
  284.  
  285.             Y[10][10] = 1.0 / 1000;
  286.             Y[10][12] = -1.0 / 1000;
  287.             Y[10][15] = 1.0;
  288.             Y[10][16] = -1.0;
  289.  
  290.             Y[11][6] = 1.0 / 0.001;
  291.             Y[11][7] = -1.0 / 0.001;
  292.             Y[11][16] = 1.0;
  293.  
  294.             Y[12][6] = -1.0 / 0.001;
  295.             Y[12][7] = 1.0 / 0.001;
  296.             Y[12][10] = -1.0 / 1000;
  297.             Y[12][12] = 1.0 / 1000 + 1.0 / RB;
  298.             Y[12][13] = -1.0 / RB;
  299.  
  300.             Y[13][3] = Cb;
  301.             Y[13][4] = -Cb;
  302.             Y[13][12] = -1.0 / RB;
  303.             Y[13][13] = 1.0 / RB + 1.0 / Ru + It / MFt * exp( (phi4 - phi5) / MFt);
  304.             Y[13][14] = -1/Ru - It / MFt * exp((phi4- phi5)/MFt);
  305.  
  306.  
  307.             Y[14][3] = -Cb;
  308.             Y[14][4] = Cb + C3;
  309.             Y[14][13] = -1.0 / Ru - (It / MFt) * ( exp(phi4 - phi5) / MFt );
  310.             Y[14][14] = 1.0 / Ru + (It / MFt) * exp( (phi4 - phi5) / MFt) + 1.0 / 1000;
  311.  
  312.             Y[15][10] = 1.0;
  313.             Y[16][10] = -1.0;
  314.             Y[16][11] = 1.0;
  315.            
  316.             I[0] = -(d_phi1 - (phi1 - phi1_past) / dt);
  317.             I[1] = -(d_phi2 - (phi2 - phi2_past) / dt);
  318.             I[2] = -(d_phi3 - (phi3 - phi3_past) / dt);
  319.             I[3] = -(d_phi4 - (phi4 - phi4_past) / dt);
  320.             I[4] = -(d_phi5 - (phi5 - phi5_past) / dt);
  321.             I[5] = -(int_phi1 - (int_phi1_past + phi1 * dt));
  322.             I[6] = -(int_phi2 - (int_phi2_past + phi2 * dt));
  323.             I[7] = -(int_phi3 - (int_phi3_past + phi3 * dt));
  324.             I[8] = -(int_phi4 - (int_phi4_past + phi4 * dt));
  325.             I[9] = -(int_phi5 - (int_phi5_past + phi5 * dt));
  326.             I[10] = -(IE1 - IE3 + (phi1 - phi3) * 1.0 / 1000);
  327.             I[11] = -(IE3 + 1.0 / 0.001 * (int_phi2 - int_phi3));
  328.             I[12] = -(-1.0 / 0.001 * (int_phi2 - int_phi3) - 1.0 / 1000 * (phi1 - phi3) + 1.0 / RB * (phi3 - phi4));
  329.             I[13] = -(-1.0 / RB * (phi3 - phi4) + Cb * (d_phi4 - d_phi5) + 1.0 / Ru * (phi4 - phi5) + It * (exp((phi4 - phi5) / MFt) - 1));
  330.             I[14] = -(-Cb * (d_phi4 - d_phi5) - 1.0 / Ru * (phi4 - phi5) - It * (exp((phi4 - phi5) / MFt) - 1) + C3 * d_phi5 + 1.0 / 1000 * phi5);
  331.             I[15] = -(phi1 - 1);
  332.             I[16] = -(phi2 - phi1 - A * sin((2 * M_PI / Freq) * t));
  333.  
  334.             DeltaX = SolveMatrixEquation(Y, I, n);
  335.  
  336.             d_phi1 += DeltaX[0];
  337.             d_phi2 += DeltaX[1];
  338.             d_phi3 += DeltaX[2];
  339.             d_phi4 += DeltaX[3];
  340.             d_phi5 += DeltaX[4];
  341.             int_phi1 += DeltaX[5];
  342.             int_phi2 += DeltaX[6];
  343.             int_phi3 += DeltaX[7];
  344.             int_phi4 += DeltaX[8];
  345.             int_phi5 += DeltaX[9];
  346.             phi1 += DeltaX[10];
  347.             phi2 += DeltaX[11];
  348.             phi3 += DeltaX[12];
  349.             phi4 += DeltaX[13];
  350.             phi5 += DeltaX[14];
  351.             IE1 += DeltaX[15];
  352.             IE3 += DeltaX[16];
  353.  
  354.             if (FindDeltaNorm(DeltaX, n) == n)
  355.             {
  356.                 NewtonIter = 0;
  357.                 Flag = 0;
  358.             }
  359.             else
  360.             {
  361.                 if (NewtonIter > MAX_NEWTON_ITER)
  362.                 {
  363.                     dt = dt / 2;
  364.                     NewtonIter = 0;
  365.                     d_phi1 = Result[(Step - 1) * (n + 1) + 1];
  366.                     d_phi2 = Result[(Step - 1) * (n + 1) + 2];
  367.                     d_phi3 = Result[(Step - 1) * (n + 1) + 3];
  368.                     d_phi4 = Result[(Step - 1) * (n + 1) + 4];
  369.                     d_phi5 = Result[(Step - 1) * (n + 1) + 5];
  370.                     int_phi1 = Result[(Step - 1) * (n + 1) + 6];
  371.                     int_phi2 = Result[(Step - 1) * (n + 1) + 7];
  372.                     int_phi3 = Result[(Step - 1) * (n + 1) + 8];
  373.                     int_phi4 = Result[(Step - 1) * (n + 1) + 9];
  374.                     int_phi5 = Result[(Step - 1) * (n + 1) + 10];
  375.                     phi1 = Result[(Step - 1) * (n + 1) + 11];
  376.                     phi2 = Result[(Step - 1) * (n + 1) + 12];
  377.                     phi3 = Result[(Step - 1) * (n + 1) + 13];
  378.                     phi4 = Result[(Step - 1) * (n + 1) + 14];
  379.                     phi5 = Result[(Step - 1) * (n + 1) + 15];
  380.                     IE1 = Result[(Step - 1) * (n + 1) + 16];
  381.                     IE3 = Result[(Step - 1) * (n + 1) + 17];
  382.  
  383.                     if (dt < DT_MIN)
  384.                     {
  385.                         cout << "Решение не сходится" << endl;
  386.                         exit(0);
  387.                     }
  388.                 }
  389.             }
  390.         }
  391.  
  392.         int_phi1_past = Result[(Step - 1) * (n + 1) + 6];
  393.         int_phi2_past = Result[(Step - 1) * (n + 1) + 7];
  394.         int_phi3_past = Result[(Step - 1) * (n + 1) + 8];
  395.         int_phi4_past = Result[(Step - 1) * (n + 1) + 9];
  396.         int_phi5_past = Result[(Step - 1) * (n + 1) + 10];
  397.         phi1_past = Result[(Step - 1) * (n + 1) + 11];
  398.         phi2_past = Result[(Step - 1) * (n + 1) + 12];
  399.         phi3_past = Result[(Step - 1) * (n + 1) + 13];
  400.         phi4_past = Result[(Step - 1) * (n + 1) + 14];
  401.         phi5_past = Result[(Step - 1) * (n + 1) + 15];
  402.  
  403.         if (Step == 1 && t == 0.0){
  404.             Step = 0;
  405.         }
  406.  
  407.         Result[Step * (n + 1)] = t;
  408.         Result[Step * (n + 1) + 1] = d_phi1;
  409.         Result[Step * (n + 1) + 2] = d_phi2;
  410.         Result[Step * (n + 1) + 3] = d_phi3;
  411.         Result[Step * (n + 1) + 4] = d_phi4;
  412.         Result[Step * (n + 1) + 5] = d_phi5;
  413.         Result[Step * (n + 1) + 6] = int_phi1;
  414.         Result[Step * (n + 1) + 7] = int_phi2;
  415.         Result[Step * (n + 1) + 8] = int_phi3;
  416.         Result[Step * (n + 1) + 9] = int_phi4;
  417.         Result[Step * (n + 1) + 10] = int_phi5;
  418.         Result[Step * (n + 1) + 11] = phi1;
  419.         Result[Step * (n + 1) + 12] = phi2;
  420.         Result[Step * (n + 1) + 13] = phi3;
  421.         Result[Step * (n + 1) + 14] = phi4;
  422.         Result[Step * (n + 1) + 15] = phi5;
  423.         Result[(Step - 1) * (n + 1) + 16] = IE1;
  424.         Result[(Step - 1) * (n + 1) + 17] = IE3;
  425.         dt_prev = dt;
  426.  
  427.         dt_prev = dt;
  428.         // t += dt;
  429.  
  430.         // cout<< t <<endl;
  431.  
  432.         // dt = DT_START;
  433.  
  434.         // расчет delta_t
  435.         if(Step > 0)
  436.         {
  437.             eps = CalculateEPS(Result, Step, dt, dt_prev, n);
  438.             cout<<eps<<endl;
  439.  
  440.             if(eps < DELTA1)
  441.             {
  442.                 t += dt;
  443.                 dt_prev = dt;
  444.                 dt = 2 * dt;
  445.  
  446.                 for(int i = 1; i < n + 1; i++)
  447.                 {
  448.                     if(Step > 1)
  449.                     {
  450.                         Result[(Step - 2) * (n + 1) + i] = Result[(Step - 1) * (n + 1) + i];
  451.                         Result[(Step - 1) * (n + 1) + i] = Result[Step * (n + 1) + i];
  452.                     }
  453.                
  454.                 }
  455.             }
  456.  
  457.             if(eps > DELTA1 && eps < DELTA2)
  458.             {
  459.                 t += dt;
  460.                 dt_prev = dt;
  461.  
  462.                 for(int i = 1; i < n + 1; i++)
  463.                 {
  464.                     if(Step > 2)
  465.                     {
  466.                         Result[(Step - 2) * (n + 1) + i] = Result[(Step - 1) * (n + 1) + i];
  467.                         Result[(Step - 1) * (n + 1) + i] = Result[Step * (n + 1) + i];
  468.                     }
  469.                 }
  470.             }
  471.  
  472.             if(eps > DELTA2)
  473.             {
  474.                 t = t;
  475.                 dt = dt / 2;
  476.  
  477.                 for(int i = 1; i < n + 1; i++)
  478.                 {
  479.                     if (Step >= 1)
  480.                     {
  481.                         Result[Step * (n + 1) + i] = Result[(Step - 1) * (n + 1) + i];
  482.                     }
  483.                 }
  484.             }
  485.         }
  486.         else
  487.         {
  488.             dt =  dt / 2;
  489.             t += dt;
  490.            
  491.         }
  492.  
  493.         Step++;
  494.         Flag = 1;
  495.     }
  496.  
  497.     MakePlotData(Result, Step, n);
  498.     // cout << totalSteps << endl;
  499.  
  500.     return 0;
  501. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!