• API
• FAQ
• Tools
• Archive
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.

Top