SHARE
TWEET

Untitled

a guest Apr 26th, 2019 64 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include <vector>
  3. #include <algorithm>
  4. #include <math.h>
  5.  
  6. auto even(const std::vector<double>& u_0, const std::vector<double>& u_1, double h, double tau){
  7.     auto N = u_0.size();
  8.     std::vector<double> u_2(N, 0);
  9.     u_2[N - 1] = 1.;
  10.     for (int i = N - 2; i > 0; --i) {
  11.         u_2[i] = 2 * h * h / tau * ((u_1[i] + h / 2.) - u_1[i]) +
  12.                 u_1[i] + (u_1[i] + h / 2.) - (u_0[i] + h / 2.);
  13.     }
  14.     u_2[0] = u_2[1] - h; // Depending on border condition approximation
  15.     return u_2;
  16. }
  17.  
  18. auto odd(const std::vector<double>& u_0, const std::vector<double>& u_1, double h, double tau){
  19.     auto N = u_0.size();
  20.     std::vector<double> u_2(N, 0);
  21.     u_2[N - 1] = 1.;
  22.     for (int i = 1; i < N - 1; ++i) {
  23.         u_2[i] = 2 * h * h / tau * (u_1[i] - (u_1[i - 1] + h / 2.)) +
  24.                 u_1[i] + (u_1[i - 1] + h / 2.) - (u_0[i - 1] + h / 2);
  25.     }
  26.     u_2[0] = u_2[1] - h; // Depending on border condition approximation
  27.     return u_2;
  28. }
  29.  
  30.  
  31.  
  32. double func(double x) {
  33.     return x - 1. + 4. * cos(3. * M_PI * x);
  34.     //return 1.;
  35. }
  36.  
  37. int main() {
  38.     int x_0 = 0;
  39.     int x_1 = 1;
  40.     int t_0 = 0;
  41.     int t_1 = 10;    // End time
  42.     double N = 10; // Space grid size
  43.     double T = 100; // Time grid size
  44.     double h = (x_1 - x_0) / N; // Space step
  45.     double tau = (t_1 - t_0) / T; // Time step
  46.  
  47.     std::vector<double> u_0(N, 0);
  48.     std::vector<double> u_1(N, 0);
  49.  
  50.     for (int i = 0; i < N; ++i) {
  51.         u_0[i] = func(i * h); // Initial at t = 0
  52.         u_1[i] = func(i * h); // Second layer at t = 0 + dt
  53.     }
  54.  
  55.  
  56.     std::cout << "u(" << 0 << ") = ";
  57.     for (int i = 0; i < N; ++i) {
  58.         std::cout << u_0[i] << " ";
  59.     }
  60.     std::cout << std::endl;
  61.  
  62.     std::cout << "u(" << tau << ") = ";
  63.     for (int i = 0; i < N; ++i) {
  64.         std::cout << u_1[i] << " ";
  65.     }
  66.     std::cout << std::endl;
  67.  
  68.     std::cout << t_0 + 2 * tau << std::endl;
  69.     std::vector<double> u_2(N, 0);
  70.     int k = 3; // First layer to count in cycle
  71.     for (double t = t_0 + 2 * tau; t <= t_1; t += tau) {
  72.         if (k % 2 == 0) {
  73.             u_2 = even(u_0, u_1, h, tau);
  74.         } else {
  75.             u_2 = odd(u_0, u_1, h, tau);
  76.         }
  77.         k++;
  78.         u_0 = u_1;
  79.         u_1 = u_2;
  80.         std::cout << "u(" << t << ") = ";
  81.         for (int i = 0; i < N; ++i) {
  82.             std::cout << u_2[i] << " ";
  83.         }
  84.         std::cout << std::endl;
  85.     }
  86.     return 0;
  87. }
RAW Paste Data
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
 
Top