Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <fstream>
- #include <cstdlib>
- #include <cmath>
- #include <iomanip>
- #include <string>
- double X_MIN = 0.0;
- double X_MAX = 1.0;
- double T_MIN = 0.0;
- double T_MAX = 0.5;
- double LAMBDA_BEZPOSREDNIE = 0.4;
- double LAMBDA_POSREDNIE = 1.0;
- double D = 1.0;
- int X_COUNT = 50;
- double H = X_MAX / (X_COUNT - 1);
- double DT_BEZPOSREDNIE = H * H * LAMBDA_BEZPOSREDNIE / D;
- int T_COUNT_BEZPOSREDNIE = (int)(T_MAX/DT_BEZPOSREDNIE + 1);
- double DT_POSREDNIE = H * H * LAMBDA_POSREDNIE / D;
- int T_COUNT_POSREDNIE = (int)(T_MAX/DT_POSREDNIE + 1);
- void Algorytm_Thomasa(double **A, double *b, double *wyn, int X_COUNT){
- for(int i = 1; i < X_COUNT; i ++){
- A[i][i] = A[i][i] - A[i][i-1] * (A[i-1][i]/A[i-1][i-1]);
- }
- for(int i = 1; i < X_COUNT; i++){
- b[i] = b[i] - (b[i-1]*A[i][i-1])/A[i-1][i-1];
- }
- wyn[X_COUNT-1] = b[X_COUNT-1]/A[X_COUNT-1][X_COUNT-1];
- for(int i = X_COUNT-2; i >= 0; i--){
- wyn[i] = (b[i] - A[i][i+1]* wyn[i+1])/ A[i][i];
- }
- }
- double Rozwiazanie_Analityczne(double x, double t){
- return 1 + exp(-(M_PI*M_PI)*D*t) * cos(M_PI*x);
- }
- double **Get_Rozw_Analit(int T_COUNT, double DT, int X_COUNT, double H){
- double x = 0, t = 0;
- auto **Rozwiazanie_Analit = new double*[T_COUNT];
- for (int i = 0; i < T_COUNT; i++){
- Rozwiazanie_Analit[i] = new double[X_COUNT];
- }
- for(int i = 0; i < T_COUNT; i++){
- for(int j = 0; j < X_COUNT; j++){
- Rozwiazanie_Analit[i][j] = Rozwiazanie_Analityczne(x, t);
- x += H;
- }
- x = 0;
- t += DT;
- }
- return Rozwiazanie_Analit;
- }
- double Warunek_Poczatkowy(double x){
- return 1 + cos(M_PI * x);
- }
- void WYKRESY_FUNKCJI(double **Analityczne, double **Metoda, int T_COUNT, double DT, const std::string &nazwa){
- auto point_interval = static_cast<int>(T_COUNT / 4);
- for(int i = 1; i < 5; i++){
- int x = 0;
- std::string name = nazwa + "_WYKRES_" + std::to_string(i) + "_" + std::to_string(i*point_interval*DT) + ".txt";
- std::fstream file;
- file.open(name, std::ios::out);
- file.precision(5);
- double* t_a = Analityczne[point_interval*i];
- double* t_m = Metoda[point_interval*i];
- for(int j = 0; j < X_COUNT; j++){
- file << std::setw(12) << j*H << std::setw(12) << t_a[j] << std::setw(12) << t_m[j] << std::endl;
- }
- file.close();
- }
- }
- void WYKRESY_BLAD_MAX(double *blad_max, int T_COUNT, double DT, const std::string &nazwa){
- std::string name = nazwa + "_WYKRES_BLAD_MAX.txt";
- std::fstream file;
- file.open(name, std::ios::out);
- file.precision(5);
- for(int i = 0; i < T_COUNT; i++){
- file << std::setw(12) << i * DT << std::setw(12) << blad_max[i] << std::endl;
- }
- file.close();
- }
- void Funkcja_Blad(double **Macierz, int T_COUNT, double DT, int X_COUNT, double H, const std::string &nazwa){
- double x = 0, t = 0;
- auto **Rozwiazanie_Analit = Get_Rozw_Analit(T_COUNT, DT, X_COUNT, H);
- auto *blad_max = new double[T_COUNT];
- auto **blad = new double*[T_COUNT];
- for (int i = 0; i < T_COUNT; i++){
- blad[i] = new double[X_COUNT];
- }
- for (int j = 0; j < T_COUNT; j++) {
- for (int i = 1; i < X_COUNT-1 ; i++) {
- blad[j][i] = fabs(Rozwiazanie_Analit[j][i] - Macierz[j][i]); //Blad bezwzgledny
- if (blad[j][i]>blad_max[j]) blad_max[j] = blad[j][i]; //blad maksymalny
- x += H;
- }
- x = 0;
- t += T_COUNT;
- }
- WYKRESY_FUNKCJI(Rozwiazanie_Analit, Macierz, T_COUNT, DT, nazwa);
- WYKRESY_BLAD_MAX(blad_max, T_COUNT, DT, nazwa);
- delete blad, blad_max;
- }
- void KMB_Laasonen(){
- std::fstream file_3;
- file_3.open("Rozw_Numeryczne_Thomas.txt",std:: ios::out);
- file_3.flags(std::ios::left);
- file_3.precision(5);
- double x = 0;
- auto *b = new double[X_COUNT];
- auto *wyniki = new double[X_COUNT];
- auto **diag = new double *[X_COUNT];
- auto **laasonen = new double *[T_COUNT_POSREDNIE];
- // Siatka Czasowo Przestrzenna
- for(int i = 0; i < T_COUNT_POSREDNIE; i++){
- laasonen[i] = new double[X_COUNT];
- }
- // Macierz A
- for(int i = 0; i < X_COUNT; i++){
- diag[i] = new double[X_COUNT];
- }
- // Zerowanie
- for(int i = 0; i < X_COUNT; i++){
- for(int j = 0; j < X_COUNT; j++){
- diag[i][j] = 0.0;
- }
- }
- // Warunek początkowy U(x,0) = 1+cos(πx)
- for(int i = 0; i < X_COUNT; i++){
- laasonen[0][i] = Warunek_Poczatkowy(i*H);
- }
- // Warunki brzegowe U(0,t) = 0, U(1,t) = 0
- for(int i = 1; i < T_COUNT_POSREDNIE; i++){
- laasonen[i][0] = 1+exp(-(M_PI*M_PI)*i*DT_POSREDNIE);
- laasonen[i][X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*DT_POSREDNIE);
- }
- //for(int i = 1; i < T_COUNT_POSREDNIE; i++){
- // laasonen[i][0] = 0;
- // laasonen[i][X_COUNT-1] = 0;
- //}
- // Dla wszystkich poziomów czasowych
- for(int i = 1; i < T_COUNT_POSREDNIE; i++){
- // Tworzenie macierzy i wektorów równania Ax=b
- diag[0][0] = -1 / H;
- diag[0][1] = 1 / H;
- b[0] = 0; // Gamma = 0
- // Dla wszystkich x'ów
- for(int j = 1; j < X_COUNT - 1; j++) {
- diag[j][j] = -(1 + (2 * LAMBDA_POSREDNIE));
- diag[j][j + 1] = LAMBDA_POSREDNIE;
- diag[j][j - 1] = LAMBDA_POSREDNIE;
- b[j] = -laasonen[i - 1][j];
- }
- b[X_COUNT - 1] = 0; // Theta = 0
- diag[X_COUNT-1][X_COUNT-1] = 1 / H;
- diag[X_COUNT-1][X_COUNT-2] = -1 / H;
- // Zastosowanie algorytmu Thomasa
- Algorytm_Thomasa(diag, b, wyniki, X_COUNT);
- // Przepisanie wyników
- for(int j = 1; j < X_COUNT - 1; j++){
- laasonen[i][j] = wyniki[j];
- }
- }
- for(int i = 0; i < T_COUNT_POSREDNIE; i++){
- for(int j = 0; j < X_COUNT; j++){
- file_3 << std::setw(12) << laasonen[i][j];
- }
- file_3 << std::endl;
- }
- Funkcja_Blad(laasonen, T_COUNT_POSREDNIE, DT_POSREDNIE, X_COUNT, H, "LAASONEN");
- file_3.close();
- delete b, wyniki, laasonen, diag;
- }
- void Funkcja_Blad_Max(double **Macierz, int T_COUNT, double DT, int X_COUNT, double H, const std::string &nazwa){
- double x = 0, t = 0;
- auto **Rozwiazanie_Analit = Get_Rozw_Analit(T_COUNT, DT, X_COUNT, H);
- auto *blad_max = new double[T_COUNT];
- auto **blad = new double*[T_COUNT];
- for (int i = 0; i < T_COUNT; i++){
- blad[i] = new double[X_COUNT];
- }
- for (int j = 0; j < T_COUNT; j++) {
- for (int i = 1; i < X_COUNT-1 ; i++) {
- blad[j][i] = fabs(Rozwiazanie_Analit[j][i] - Macierz[j][i]); //Blad bezwzgledny
- if (blad[j][i]>blad_max[j]) blad_max[j] = blad[j][i]; //blad maksymalny
- x += H;
- }
- x = 0;
- t += T_COUNT;
- }
- std::string name = nazwa + "_MAX.txt";
- std::fstream file;
- file.open(name, std::ios::out|std::ios::app);
- file.precision(6);
- file << std::setw(12) << log10(H) << std::setw(12) << log10(blad_max[T_COUNT-1]) << std::endl;
- file.close();
- delete blad_max;
- for(int i = 0; i < T_COUNT; i++){
- delete blad[i];
- }
- delete blad;
- for (int i = 0; i < T_COUNT; i++){
- delete Rozwiazanie_Analit[i];
- }
- delete Rozwiazanie_Analit;
- }
- void KMB_Laasonen_blad_tmax(){
- std::string name = "LAASONEN_MAX.txt";
- std::fstream file;
- file.open(name, std::ios::out|std::ios::trunc);
- file.close();
- //H - 1/10 - 1/200
- for(int loop = 10; loop < 301; loop += 10){
- auto LOCAL_H = (1 / (double)loop);
- auto LOCAL_X_COUNT = static_cast<int>(1 / LOCAL_H + 1);
- double LOCAL_DT = LOCAL_H * LOCAL_H * LAMBDA_POSREDNIE / D;
- auto LOCAL_T_COUNT = static_cast<int>(T_MAX/LOCAL_DT + 1);
- //-------------------------------------------------------
- double x = 0;
- double t = 0;
- auto *b = new double[LOCAL_X_COUNT];
- auto *wyniki = new double[LOCAL_X_COUNT];
- auto **diag = new double *[LOCAL_X_COUNT];
- auto **laasonen = new double *[LOCAL_T_COUNT];
- // Siatka Czasowo Przestrzenna
- for(int i = 0; i < LOCAL_T_COUNT; i++){
- laasonen[i] = new double[LOCAL_X_COUNT];
- }
- // Macierz A
- for(int i = 0; i < LOCAL_X_COUNT; i++){
- diag[i] = new double[LOCAL_X_COUNT];
- }
- // Zerowanie
- for(int i = 0; i < LOCAL_X_COUNT; i++){
- for(int j = 0; j < LOCAL_X_COUNT; j++){
- diag[i][j] = 0.0;
- }
- }
- // Warunek początkowy U(x,0) = 1+cos(πx)
- for(int i = 0; i < LOCAL_X_COUNT; i++){
- laasonen[0][i] = Warunek_Poczatkowy(i*LOCAL_H);
- }
- // Warunki brzegowe U(0,t) = 0, U(1,t) = 0
- for(int i = 1; i < LOCAL_T_COUNT; i++){
- laasonen[i][0] = 1+exp(-(M_PI*M_PI)*i*LOCAL_DT);
- laasonen[i][LOCAL_X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*LOCAL_DT);
- }
- // Dla wszystkich poziomów czasowych
- for(int i = 1; i < LOCAL_T_COUNT; i++){
- // Tworzenie macierzy i wektorów równania Ax=b
- diag[0][0] = -1 / LOCAL_H;
- diag[0][1] = 1 / LOCAL_H;
- b[0] = 0; // Gamma = 0
- // Dla wszystkich x'ów
- for(int j = 1; j < LOCAL_X_COUNT - 1; j++) {
- diag[j][j] = -(1 + (2 * LAMBDA_POSREDNIE));
- diag[j][j + 1] = LAMBDA_POSREDNIE;
- diag[j][j - 1] = LAMBDA_POSREDNIE;
- b[j] = -laasonen[i - 1][j];
- }
- b[LOCAL_X_COUNT - 1] = 0; // Theta = 0
- diag[LOCAL_X_COUNT-1][LOCAL_X_COUNT-1] = 1 / LOCAL_H;
- diag[LOCAL_X_COUNT-1][LOCAL_X_COUNT-2] = -1 / LOCAL_H;
- // Zastosowanie algorytmu Thomasa
- Algorytm_Thomasa(diag, b, wyniki, LOCAL_X_COUNT);
- // Przepisanie wyników
- for(int j = 1; j < LOCAL_X_COUNT - 1; j++){
- laasonen[i][j] = wyniki[j];
- }
- }
- // W Tym miejscu mamy wyniki dla każdej iteracji
- Funkcja_Blad_Max(laasonen,LOCAL_T_COUNT, LOCAL_DT, LOCAL_X_COUNT, LOCAL_H, "LAASONEN");
- delete b, wyniki;
- for(int i = 0; i < LOCAL_X_COUNT; i++){
- delete diag[i];
- }
- delete diag;
- for(int i = 0; i < LOCAL_T_COUNT; i++){
- delete laasonen[i];
- }
- delete laasonen;
- }
- }
- void Bezposrednia(){
- std::fstream file_3;
- file_3.open("Rozw_Numeryczne_Bezpo.txt",std:: ios::out);
- file_3.flags(std::ios::left);
- file_3.precision(5);
- double t = 0, x = 0;
- auto **bezpo = new double*[T_COUNT_BEZPOSREDNIE];
- for(int i = 0; i < T_COUNT_BEZPOSREDNIE; i++){
- bezpo[i] = new double[X_COUNT];
- }
- // Warunek poczatkowy
- for(int i = 0; i < X_COUNT; i++){
- bezpo[0][i] = Warunek_Poczatkowy(i*H);
- }
- for(int i = 1; i < T_COUNT_BEZPOSREDNIE; i++){
- bezpo[i][0] = 1+exp(-(M_PI*M_PI)*i*DT_BEZPOSREDNIE);
- bezpo[i][X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*DT_BEZPOSREDNIE);
- }
- for(int i = 0; i < T_COUNT_BEZPOSREDNIE-1; i++) {
- for(int j = 1; j < X_COUNT-1; j++) {
- bezpo[i+1][j] = LAMBDA_BEZPOSREDNIE * bezpo[i][j-1] + (1-2*LAMBDA_BEZPOSREDNIE) * bezpo[i][j] + LAMBDA_BEZPOSREDNIE * bezpo[i][j+1];
- }
- //bezpo[i+1][0] = bezpo[i+1][1];
- //bezpo[i+1][X_COUNT-1] = bezpo[i+1][X_COUNT-2];
- }
- for(int i = 0; i < T_COUNT_BEZPOSREDNIE; i++){
- for(int j = 0; j < X_COUNT; j++){
- file_3 << std::setw(12) << bezpo[i][j];
- }
- file_3 << std::endl;
- }
- Funkcja_Blad(bezpo, T_COUNT_BEZPOSREDNIE, DT_BEZPOSREDNIE, X_COUNT, H, "BEZPOSREDNIA");
- file_3.close();
- }
- void Bezposrednia_blad_tmax(){
- std::string name = "BEZPOSREDNIA_MAX.txt";
- std::fstream file;
- file.open(name, std::ios::out|std::ios::trunc);
- file.close();
- for(int loop = 10; loop < 301; loop += 10){
- auto LOCAL_H = (1 / (double)loop);
- auto LOCAL_X_COUNT = static_cast<int>(1 / LOCAL_H + 1);
- double LOCAL_DT = LOCAL_H * LOCAL_H * LAMBDA_BEZPOSREDNIE / D;
- auto LOCAL_T_COUNT = static_cast<int>(T_MAX/LOCAL_DT + 1);
- //-----------------------------------------------------
- double t = 0, x = 0;
- auto **bezpo = new double*[LOCAL_T_COUNT];
- for(int i = 0; i < LOCAL_T_COUNT; i++){
- bezpo[i] = new double[LOCAL_X_COUNT];
- }
- // Warunek poczatkowy
- for(int i = 0; i < LOCAL_X_COUNT; i++){
- bezpo[0][i] = Warunek_Poczatkowy(i*LOCAL_H);
- }
- for(int i = 1; i < LOCAL_T_COUNT; i++){
- bezpo[i][0] = 1+exp(-(M_PI*M_PI)*i*LOCAL_DT);
- bezpo[i][LOCAL_X_COUNT-1] = 1-exp(-(M_PI*M_PI)*i*LOCAL_DT);
- }
- for(int i = 0; i < LOCAL_T_COUNT-1; i++) {
- for(int j = 1; j < LOCAL_X_COUNT-1; j++) {
- bezpo[i+1][j] = LAMBDA_BEZPOSREDNIE * bezpo[i][j-1] + (1-2*LAMBDA_BEZPOSREDNIE) * bezpo[i][j] + LAMBDA_BEZPOSREDNIE * bezpo[i][j+1];
- }
- }
- //-----------------------------------------------------
- Funkcja_Blad_Max(bezpo, LOCAL_T_COUNT, LOCAL_DT, LOCAL_X_COUNT, LOCAL_H, "BEZPOSREDNIA");
- for(int i = 0; i < LOCAL_T_COUNT; i++){
- delete bezpo[i];
- }
- delete bezpo;
- }
- }
- int main() {
- KMB_Laasonen();
- KMB_Laasonen_blad_tmax();
- Bezposrednia();
- Bezposrednia_blad_tmax();
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement