API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jun 18th, 2018
90
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 13.84 KB | None | 0 0
raw download clone embed print report
  1. #include <iostream>
  2. #include <fstream>
  3. #include <cstdlib>
  4. #include <cmath>
  5. #include <iomanip>
  6. #include <string>
  7.  
  8. double X_MIN = 0.0;
  9. double X_MAX = 1.0;
  10.  
  11. double T_MIN = 0.0;
  12. double T_MAX = 0.5;
  13.  
  14. double LAMBDA_BEZPOSREDNIE = 0.4;
  15. double LAMBDA_POSREDNIE = 1.0;
  16.  
  17. double D = 1.0;
  18.  
  19. int X_COUNT = 50;
  20. double H = X_MAX / (X_COUNT - 1);
  21.  
  22. double DT_BEZPOSREDNIE = H * H * LAMBDA_BEZPOSREDNIE / D;
  23. int T_COUNT_BEZPOSREDNIE = (int)(T_MAX/DT_BEZPOSREDNIE + 1);
  24.  
  25. double DT_POSREDNIE = H * H * LAMBDA_POSREDNIE / D;
  26. int T_COUNT_POSREDNIE = (int)(T_MAX/DT_POSREDNIE + 1);
  27.  
  28. void Algorytm_Thomasa(double **A, double *b, double *wyn, int X_COUNT){
  29.  
  30.     for(int i = 1; i < X_COUNT; i ++){
  31.  
  32.         A[i][i] = A[i][i] - A[i][i-1] * (A[i-1][i]/A[i-1][i-1]);
  33.  
  34.     }
  35.  
  36.     for(int i = 1; i < X_COUNT; i++){
  37.  
  38.         b[i] = b[i] - (b[i-1]*A[i][i-1])/A[i-1][i-1];
  39.  
  40.     }
  41.  
  42.     wyn[X_COUNT-1] = b[X_COUNT-1]/A[X_COUNT-1][X_COUNT-1];
  43.  
  44.     for(int i = X_COUNT-2; i >= 0; i--){
  45.  
  46.         wyn[i] = (b[i] - A[i][i+1]* wyn[i+1])/ A[i][i];
  47.  
  48.     }
  49.  
  50. }
  51.  
  52. double Rozwiazanie_Analityczne(double x, double t){
  53.  
  54.     return 1 + exp(-(M_PI*M_PI)*D*t) * cos(M_PI*x);
  55.  
  56. }
  57.  
  58. double **Get_Rozw_Analit(int T_COUNT, double DT, int X_COUNT, double H){
  59.  
  60.     double x = 0, t = 0;
  61.  
  62.     auto **Rozwiazanie_Analit = new double*[T_COUNT];
  63.     for (int i = 0; i < T_COUNT; i++){
  64.  
  65.         Rozwiazanie_Analit[i] = new double[X_COUNT];
  66.  
  67.     }
  68.  
  69.     for(int i = 0; i < T_COUNT; i++){
  70.  
  71.         for(int j = 0; j < X_COUNT; j++){
  72.  
  73.             Rozwiazanie_Analit[i][j] = Rozwiazanie_Analityczne(x, t);
  74.             x += H;
  75.  
  76.         }
  77.  
  78.         x = 0;
  79.         t += DT;
  80.  
  81.     }
  82.  
  83.     return Rozwiazanie_Analit;
  84.  
  85. }
  86.  
  87. double Warunek_Poczatkowy(double x){
  88.  
  89.     return 1 + cos(M_PI * x);
  90.  
  91. }
  92.  
  93. void WYKRESY_FUNKCJI(double **Analityczne, double **Metoda, int T_COUNT, double DT, const std::string &nazwa){
  94.  
  95.     auto point_interval = static_cast<int>(T_COUNT / 4);
  96.  
  97.     for(int i = 1; i < 5; i++){
  98.  
  99.         int x = 0;
  100.  
  101.         std::string name = nazwa + "_WYKRES_" + std::to_string(i) + "_" + std::to_string(i*point_interval*DT) + ".txt";
  102.         std::fstream file;
  103.         file.open(name, std::ios::out);
  104.         file.precision(5);
  105.  
  106.         double* t_a = Analityczne[point_interval*i];
  107.         double* t_m = Metoda[point_interval*i];
  108.  
  109.         for(int j = 0; j < X_COUNT; j++){
  110.  
  111.             file << std::setw(12) << j*H << std::setw(12) << t_a[j] << std::setw(12) << t_m[j] << std::endl;
  112.  
  113.         }
  114.  
  115.         file.close();
  116.  
  117.     }
  118. }
  119.  
  120. void WYKRESY_BLAD_MAX(double *blad_max, int T_COUNT, double DT, const std::string &nazwa){
  121.  
  122.     std::string name = nazwa + "_WYKRES_BLAD_MAX.txt";
  123.     std::fstream file;
  124.     file.open(name, std::ios::out);
  125.     file.precision(5);
  126.  
  127.     for(int i = 0; i < T_COUNT; i++){
  128.         file << std::setw(12) << i * DT << std::setw(12) << blad_max[i] << std::endl;
  129.     }
  130.  
  131.     file.close();
  132.  
  133. }
  134.  
  135. void Funkcja_Blad(double **Macierz, int T_COUNT, double DT, int X_COUNT, double H, const std::string &nazwa){
  136.  
  137.     double x = 0, t = 0;
  138.  
  139.     auto **Rozwiazanie_Analit = Get_Rozw_Analit(T_COUNT, DT, X_COUNT, H);
  140.  
  141.     auto *blad_max = new double[T_COUNT];
  142.  
  143.     auto **blad = new double*[T_COUNT];
  144.     for (int i = 0; i < T_COUNT; i++){
  145.  
  146.         blad[i] = new double[X_COUNT];
  147.  
  148.     }
  149.  
  150.     for (int j = 0; j < T_COUNT; j++) {
  151.  
  152.         for (int i = 1; i < X_COUNT-1 ; i++) {
  153.  
  154.             blad[j][i] = fabs(Rozwiazanie_Analit[j][i] - Macierz[j][i]); //Blad bezwzgledny
  155.  
  156.             if (blad[j][i]>blad_max[j]) blad_max[j] = blad[j][i]; //blad maksymalny
  157.             x += H;
  158.         }
  159.  
  160.         x = 0;
  161.         t += T_COUNT;
  162.  
  163.     }
  164.  
  165.     WYKRESY_FUNKCJI(Rozwiazanie_Analit, Macierz, T_COUNT, DT, nazwa);
  166.     WYKRESY_BLAD_MAX(blad_max, T_COUNT, DT, nazwa);
  167.  
  168.     delete blad, blad_max;
  169.  
  170. }
  171.  
  172. void KMB_Laasonen(){
  173.  
  174.     std::fstream file_3;
  175.     file_3.open("Rozw_Numeryczne_Thomas.txt",std:: ios::out);
  176.     file_3.flags(std::ios::left);
  177.     file_3.precision(5);
  178.  
  179.     double x = 0;
  180.     auto *b = new double[X_COUNT];
  181.     auto *wyniki = new double[X_COUNT];
  182.     auto **diag = new double *[X_COUNT];
  183.     auto **laasonen = new double *[T_COUNT_POSREDNIE];
  184.  
  185.     // Siatka Czasowo Przestrzenna
  186.     for(int i = 0; i < T_COUNT_POSREDNIE; i++){
  187.         laasonen[i] = new double[X_COUNT];
  188.     }
  189.  
  190.     // Macierz A
  191.     for(int i = 0; i < X_COUNT; i++){
  192.         diag[i] = new double[X_COUNT];
  193.     }
  194.  
  195.     // Zerowanie
  196.     for(int i = 0; i < X_COUNT; i++){
  197.         for(int j = 0; j < X_COUNT; j++){
  198.             diag[i][j] = 0.0;
  199.         }
  200.     }
  201.  
  202.     // Warunek początkowy U(x,0) = 1+cos(πx)
  203.  
  204.     for(int i = 0; i < X_COUNT; i++){
  205.         laasonen[0][i] = Warunek_Poczatkowy(i*H);
  206.     }
  207.  
  208.     // Warunki brzegowe U(0,t) = 0, U(1,t) = 0
  209.  
  210.     for(int i = 1; i < T_COUNT_POSREDNIE; i++){
  211.  
  212.         laasonen[i][0] = 1+exp(-(M_PI*M_PI)*i*DT_POSREDNIE);
  213.         laasonen[i][X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*DT_POSREDNIE);
  214.  
  215.     }
  216.  
  217.     //for(int i = 1; i < T_COUNT_POSREDNIE; i++){
  218.     //    laasonen[i][0] = 0;
  219.     //    laasonen[i][X_COUNT-1] = 0;
  220.     //}
  221.  
  222.     // Dla wszystkich poziomów czasowych
  223.     for(int i = 1; i < T_COUNT_POSREDNIE; i++){
  224.  
  225.         // Tworzenie macierzy i wektorów równania Ax=b
  226.  
  227.         diag[0][0] = -1 / H;
  228.         diag[0][1] = 1 / H;
  229.  
  230.         b[0] = 0; // Gamma = 0
  231.  
  232.         // Dla wszystkich x'ów
  233.         for(int j = 1; j < X_COUNT - 1; j++) {
  234.  
  235.             diag[j][j] = -(1 + (2 * LAMBDA_POSREDNIE));
  236.             diag[j][j + 1] = LAMBDA_POSREDNIE;
  237.             diag[j][j - 1] = LAMBDA_POSREDNIE;
  238.             b[j] = -laasonen[i - 1][j];
  239.  
  240.         }
  241.  
  242.         b[X_COUNT - 1] = 0; // Theta = 0
  243.  
  244.         diag[X_COUNT-1][X_COUNT-1] = 1 / H;
  245.         diag[X_COUNT-1][X_COUNT-2] = -1 / H;
  246.  
  247.  
  248.         // Zastosowanie algorytmu Thomasa
  249.  
  250.         Algorytm_Thomasa(diag, b, wyniki, X_COUNT);
  251.  
  252.         // Przepisanie wyników
  253.  
  254.         for(int j = 1; j < X_COUNT - 1; j++){
  255.  
  256.             laasonen[i][j] = wyniki[j];
  257.  
  258.         }
  259.     }
  260.  
  261.     for(int i = 0; i < T_COUNT_POSREDNIE; i++){
  262.         for(int j = 0; j < X_COUNT; j++){
  263.             file_3 << std::setw(12) << laasonen[i][j];
  264.         }
  265.         file_3 << std::endl;
  266.     }
  267.  
  268.     Funkcja_Blad(laasonen, T_COUNT_POSREDNIE, DT_POSREDNIE, X_COUNT, H, "LAASONEN");
  269.  
  270.     file_3.close();
  271.  
  272.     delete b, wyniki, laasonen, diag;
  273.  
  274. }
  275.  
  276. void Funkcja_Blad_Max(double **Macierz, int T_COUNT, double DT, int X_COUNT, double H, const std::string &nazwa){
  277.  
  278.     double x = 0, t = 0;
  279.  
  280.     auto **Rozwiazanie_Analit = Get_Rozw_Analit(T_COUNT, DT, X_COUNT, H);
  281.  
  282.     auto *blad_max = new double[T_COUNT];
  283.  
  284.     auto **blad = new double*[T_COUNT];
  285.     for (int i = 0; i < T_COUNT; i++){
  286.  
  287.         blad[i] = new double[X_COUNT];
  288.  
  289.     }
  290.  
  291.     for (int j = 0; j < T_COUNT; j++) {
  292.  
  293.         for (int i = 1; i < X_COUNT-1 ; i++) {
  294.  
  295.             blad[j][i] = fabs(Rozwiazanie_Analit[j][i] - Macierz[j][i]); //Blad bezwzgledny
  296.  
  297.             if (blad[j][i]>blad_max[j]) blad_max[j] = blad[j][i]; //blad maksymalny
  298.             x += H;
  299.         }
  300.  
  301.         x = 0;
  302.         t += T_COUNT;
  303.  
  304.     }
  305.  
  306.     std::string name = nazwa + "_MAX.txt";
  307.     std::fstream file;
  308.     file.open(name, std::ios::out|std::ios::app);
  309.     file.precision(6);
  310.  
  311.     file << std::setw(12) << log10(H) << std::setw(12) << log10(blad_max[T_COUNT-1]) << std::endl;
  312.  
  313.     file.close();
  314.  
  315.     delete blad_max;
  316.  
  317.     for(int i = 0; i < T_COUNT; i++){
  318.         delete blad[i];
  319.     }
  320.     delete blad;
  321.  
  322.     for (int i = 0; i < T_COUNT; i++){
  323.         delete Rozwiazanie_Analit[i];
  324.     }
  325.     delete Rozwiazanie_Analit;
  326.  
  327. }
  328.  
  329. void KMB_Laasonen_blad_tmax(){
  330.  
  331.     std::string name = "LAASONEN_MAX.txt";
  332.     std::fstream file;
  333.     file.open(name, std::ios::out|std::ios::trunc);
  334.     file.close();
  335.  
  336.     //H - 1/10 - 1/200
  337.  
  338.     for(int loop = 10; loop < 301; loop += 10){
  339.  
  340.         auto LOCAL_H = (1 / (double)loop);
  341.         auto LOCAL_X_COUNT = static_cast<int>(1 / LOCAL_H + 1);
  342.         double LOCAL_DT = LOCAL_H * LOCAL_H * LAMBDA_POSREDNIE / D;
  343.         auto LOCAL_T_COUNT = static_cast<int>(T_MAX/LOCAL_DT + 1);
  344.  
  345.         //-------------------------------------------------------
  346.         double x = 0;
  347.         double t = 0;
  348.         auto *b = new double[LOCAL_X_COUNT];
  349.         auto *wyniki = new double[LOCAL_X_COUNT];
  350.         auto **diag = new double *[LOCAL_X_COUNT];
  351.         auto **laasonen = new double *[LOCAL_T_COUNT];
  352.  
  353.  
  354.         // Siatka Czasowo Przestrzenna
  355.         for(int i = 0; i < LOCAL_T_COUNT; i++){
  356.             laasonen[i] = new double[LOCAL_X_COUNT];
  357.         }
  358.  
  359.         // Macierz A
  360.         for(int i = 0; i < LOCAL_X_COUNT; i++){
  361.             diag[i] = new double[LOCAL_X_COUNT];
  362.         }
  363.  
  364.         // Zerowanie
  365.         for(int i = 0; i < LOCAL_X_COUNT; i++){
  366.             for(int j = 0; j < LOCAL_X_COUNT; j++){
  367.                 diag[i][j] = 0.0;
  368.             }
  369.         }
  370.  
  371.         // Warunek początkowy U(x,0) = 1+cos(πx)
  372.  
  373.         for(int i = 0; i < LOCAL_X_COUNT; i++){
  374.             laasonen[0][i] = Warunek_Poczatkowy(i*LOCAL_H);
  375.         }
  376.  
  377.         // Warunki brzegowe U(0,t) = 0, U(1,t) = 0
  378.  
  379.         for(int i = 1; i < LOCAL_T_COUNT; i++){
  380.  
  381.             laasonen[i][0] = 1+exp(-(M_PI*M_PI)*i*LOCAL_DT);
  382.             laasonen[i][LOCAL_X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*LOCAL_DT);
  383.  
  384.         }
  385.  
  386.         // Dla wszystkich poziomów czasowych
  387.         for(int i = 1; i < LOCAL_T_COUNT; i++){
  388.  
  389.             // Tworzenie macierzy i wektorów równania Ax=b
  390.  
  391.             diag[0][0] = -1 / LOCAL_H;
  392.             diag[0][1] = 1 / LOCAL_H;
  393.  
  394.             b[0] = 0; // Gamma = 0
  395.  
  396.             // Dla wszystkich x'ów
  397.             for(int j = 1; j < LOCAL_X_COUNT - 1; j++) {
  398.  
  399.                 diag[j][j] = -(1 + (2 * LAMBDA_POSREDNIE));
  400.                 diag[j][j + 1] = LAMBDA_POSREDNIE;
  401.                 diag[j][j - 1] = LAMBDA_POSREDNIE;
  402.                 b[j] = -laasonen[i - 1][j];
  403.  
  404.             }
  405.  
  406.             b[LOCAL_X_COUNT - 1] = 0; // Theta = 0
  407.  
  408.             diag[LOCAL_X_COUNT-1][LOCAL_X_COUNT-1] = 1 / LOCAL_H;
  409.             diag[LOCAL_X_COUNT-1][LOCAL_X_COUNT-2] = -1 / LOCAL_H;
  410.  
  411.  
  412.             // Zastosowanie algorytmu Thomasa
  413.  
  414.             Algorytm_Thomasa(diag, b, wyniki, LOCAL_X_COUNT);
  415.  
  416.             // Przepisanie wyników
  417.  
  418.             for(int j = 1; j < LOCAL_X_COUNT - 1; j++){
  419.  
  420.                 laasonen[i][j] = wyniki[j];
  421.  
  422.             }
  423.  
  424.         }
  425.  
  426.         // W Tym miejscu mamy wyniki dla każdej iteracji
  427.  
  428.         Funkcja_Blad_Max(laasonen,LOCAL_T_COUNT, LOCAL_DT, LOCAL_X_COUNT, LOCAL_H, "LAASONEN");
  429.  
  430.         delete b, wyniki;
  431.  
  432.         for(int i = 0; i < LOCAL_X_COUNT; i++){
  433.             delete diag[i];
  434.         }
  435.         delete diag;
  436.  
  437.         for(int i = 0; i < LOCAL_T_COUNT; i++){
  438.             delete laasonen[i];
  439.         }
  440.         delete laasonen;
  441.  
  442.     }
  443.  
  444. }
  445.  
  446. void Bezposrednia(){
  447.  
  448.     std::fstream file_3;
  449.     file_3.open("Rozw_Numeryczne_Bezpo.txt",std:: ios::out);
  450.     file_3.flags(std::ios::left);
  451.     file_3.precision(5);
  452.  
  453.     double t = 0, x = 0;
  454.  
  455.     auto **bezpo = new double*[T_COUNT_BEZPOSREDNIE];
  456.     for(int i = 0; i < T_COUNT_BEZPOSREDNIE; i++){
  457.         bezpo[i] = new double[X_COUNT];
  458.     }
  459.  
  460.     // Warunek poczatkowy
  461.     for(int i = 0; i < X_COUNT; i++){
  462.         bezpo[0][i] = Warunek_Poczatkowy(i*H);
  463.     }
  464.  
  465.     for(int i = 1; i < T_COUNT_BEZPOSREDNIE; i++){
  466.  
  467.         bezpo[i][0] = 1+exp(-(M_PI*M_PI)*i*DT_BEZPOSREDNIE);
  468.         bezpo[i][X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*DT_BEZPOSREDNIE);
  469.  
  470.     }
  471.  
  472.     for(int i = 0; i < T_COUNT_BEZPOSREDNIE-1; i++) {
  473.  
  474.         for(int j = 1; j < X_COUNT-1; j++) {
  475.  
  476.             bezpo[i+1][j] = LAMBDA_BEZPOSREDNIE * bezpo[i][j-1] + (1-2*LAMBDA_BEZPOSREDNIE) * bezpo[i][j] + LAMBDA_BEZPOSREDNIE * bezpo[i][j+1];
  477.  
  478.         }
  479.  
  480.         //bezpo[i+1][0] = bezpo[i+1][1];
  481.         //bezpo[i+1][X_COUNT-1] = bezpo[i+1][X_COUNT-2];
  482.  
  483.     }
  484.  
  485.     for(int i = 0; i < T_COUNT_BEZPOSREDNIE; i++){
  486.         for(int j = 0; j < X_COUNT; j++){
  487.             file_3 << std::setw(12) << bezpo[i][j];
  488.         }
  489.         file_3 << std::endl;
  490.     }
  491.  
  492.     Funkcja_Blad(bezpo, T_COUNT_BEZPOSREDNIE, DT_BEZPOSREDNIE, X_COUNT, H, "BEZPOSREDNIA");
  493.  
  494.     file_3.close();
  495.  
  496. }
  497.  
  498. void Bezposrednia_blad_tmax(){
  499.  
  500.     std::string name = "BEZPOSREDNIA_MAX.txt";
  501.     std::fstream file;
  502.     file.open(name, std::ios::out|std::ios::trunc);
  503.     file.close();
  504.  
  505.     for(int loop = 10; loop < 301; loop += 10){
  506.  
  507.         auto LOCAL_H = (1 / (double)loop);
  508.         auto LOCAL_X_COUNT = static_cast<int>(1 / LOCAL_H + 1);
  509.         double LOCAL_DT = LOCAL_H * LOCAL_H * LAMBDA_BEZPOSREDNIE / D;
  510.         auto LOCAL_T_COUNT = static_cast<int>(T_MAX/LOCAL_DT + 1);
  511.  
  512.         //-----------------------------------------------------
  513.         double t = 0, x = 0;
  514.  
  515.         auto **bezpo = new double*[LOCAL_T_COUNT];
  516.         for(int i = 0; i < LOCAL_T_COUNT; i++){
  517.             bezpo[i] = new double[LOCAL_X_COUNT];
  518.         }
  519.  
  520.         // Warunek poczatkowy
  521.         for(int i = 0; i < LOCAL_X_COUNT; i++){
  522.             bezpo[0][i] = Warunek_Poczatkowy(i*LOCAL_H);
  523.         }
  524.  
  525.         for(int i = 1; i < LOCAL_T_COUNT; i++){
  526.  
  527.             bezpo[i][0] = 1+exp(-(M_PI*M_PI)*i*LOCAL_DT);
  528.             bezpo[i][LOCAL_X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*LOCAL_DT);
  529.  
  530.         }
  531.  
  532.         for(int i = 0; i < LOCAL_T_COUNT-1; i++) {
  533.  
  534.             for(int j = 1; j < LOCAL_X_COUNT-1; j++) {
  535.  
  536.                 bezpo[i+1][j] = LAMBDA_BEZPOSREDNIE * bezpo[i][j-1] + (1-2*LAMBDA_BEZPOSREDNIE) * bezpo[i][j] + LAMBDA_BEZPOSREDNIE * bezpo[i][j+1];
  537.  
  538.             }
  539.  
  540.         }
  541.         //-----------------------------------------------------
  542.  
  543.         Funkcja_Blad_Max(bezpo, LOCAL_T_COUNT, LOCAL_DT, LOCAL_X_COUNT, LOCAL_H, "BEZPOSREDNIA");
  544.  
  545.         for(int i = 0; i < LOCAL_T_COUNT; i++){
  546.             delete bezpo[i];
  547.         }
  548.         delete bezpo;
  549.  
  550.     }
  551.  
  552. }
  553.  
  554. int main() {
  555.     KMB_Laasonen();
  556.     KMB_Laasonen_blad_tmax();
  557.  
  558.     Bezposrednia();
  559.     Bezposrednia_blad_tmax();
  560.     return 0;
  561. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
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!