#include<iostream> #include<math.h> #include <fstream> #include "calerf.h" void rozw_analityczne (DOUBLE **analit,INT ile_x, INT ile_t, DOUBLE h, DOUBLE dt); void lassonen(DOUBLE **lasson,INT ile_x, INT ile_t, DOUBLE h, DOUBLE dt); void blad (DOUBLE **blad, DOUBLE **analit, DOUBLE **lasson,INT ile_x, INT ile_t, DOUBLE h, DOUBLE dt); void thomas( DOUBLE *a, DOUBLE *d, DOUBLE *c, DOUBLE *b, DOUBLE *u, INT n ); USING namespace std; DOUBLE D = 1; //wspolczynnik ciepla DOUBLE t_max = 2; // MAX t (na osi y) DOUBLE x_max = 6*SQRT(D*t_max); // MAX na osi x DOUBLE lambda = 1; //lambda dla metody bezposredniej class metoda { public: DOUBLE **laasonen, **analit, **blad, *blad_max; metoda (); INT ile_x; DOUBLE h; DOUBLE dt; INT ile_t; void rozw_analityczne (); void rozw_laasonen (); void thomas ( DOUBLE *up, DOUBLE *diag, DOUBLE *low, DOUBLE *b, DOUBLE *wyn); void f_blad (); void zapis_do_pliku_m (CHAR* nazwa, DOUBLE** macierz); void zapis_do_pliku_w (CHAR* nazwa, DOUBLE* wektor); }; metoda::metoda () //konstruktor wylicza przedzialy i alokuje pamiec na macierze rozwiazan { ile_x = 40; h = x_max / ile_x; dt = h*h*lambda / D; ile_t = static_cast<int> ( (t_max/dt) + 1 ); laasonen = NEW DOUBLE*[ile_t]; FOR (INT i=0;i<ile_t;i++) laasonen[i] = NEW DOUBLE[ile_x]; analit = NEW DOUBLE*[ile_t]; FOR (INT i=0;i<ile_t;i++) analit[i] = NEW DOUBLE[ile_x]; blad = NEW DOUBLE*[ile_t]; FOR (INT i=0;i<ile_t;i++) blad[i] = NEW DOUBLE[ile_x]; blad_max = NEW DOUBLE[ile_t]; } void metoda::rozw_analityczne () { DOUBLE t = 0,x = 0; FOR (INT j = 0; j < ile_t; j++) { FOR (INT i = 0; i < ile_x ; i++) { analit[j][i] = erfc (x / (2*SQRT(D*t))); // rozwiazanie analityczne x = x+h; } x = 0; t = t+dt; } zapis_do_pliku_m ("rozw_analityczne_l.txt", analit); } void metoda::rozw_laasonen () { DOUBLE x=0,t=0; DOUBLE *up = NEW DOUBLE[ile_x]; DOUBLE *b = NEW DOUBLE[ile_x]; DOUBLE *low = NEW DOUBLE[ile_x]; DOUBLE *diag = NEW DOUBLE[ile_x]; DOUBLE *wyn = NEW DOUBLE[ile_x]; FOR( INT i = 0; i < ile_x; i++ ) //warunek poczatkowy laasonen[0][i] = 0; FOR( INT i = 0; i < ile_t; i++ ) // warunki brzegowe laasonen[i][0] = 1; FOR( INT i = 0; i < ile_t; i++ ) laasonen[i][ile_x-1] = 0; /*Wypelnianie macierzy */ FOR( INT k = 1; k < ile_t; k++ ) { up[0] = 0; diag[0] = 1; b[0] = laasonen[k-1][0]; FOR( INT i = 1; i < ile_x-1; i++ ) { up[i] = lambda; diag[i] = -( 1 + (2*lambda) ); low[i] = lambda; b[i] = -laasonen[k-1][i]; } up[ile_x-1] = 0; diag[ile_x-1] = 1; b[ile_x-1] = laasonen[k-1][ile_x-1]; thomas( up, diag, low, b, wyn); FOR( INT i = 1; i < ile_x-1; i++ ) laasonen[k][i] = wyn[i]; // wyliczenie wektora rozwiazan Thomasem i zapisanie do kolejnego wiersza } zapis_do_pliku_m ("rozw_metoda_laasonen.txt", laasonen); } void metoda::thomas( DOUBLE *up, DOUBLE *diag, DOUBLE *low, DOUBLE *b, DOUBLE *wyn) { INT i; low[0] = low[0] / diag[0]; b[0] = b[0] / diag[0]; FOR( i = 1; i < ile_x; i++ ) { DOUBLE id = (diag[i] - low[i-1] * up[i]); low[i] = low[i] / id; b[i] = ( b[i] - b[i-1] * up[i] ) / id; } // wyliczanie wyniku wyn[ile_x - 1] = b[ile_x - 1]; FOR( i = ile_x - 2; i >= 0; i-- ) wyn[i] = b[i] - low[i] * wyn[i+1]; } void metoda::f_blad () { DOUBLE x = 0,t = 0; FOR (INT j = 0; j < ile_t; j++) { FOR (INT i = 0; i < ile_x ; i++) { blad[j][i]=fabs(analit[j][i] - laasonen[j][i]); IF (blad[j][i]>blad_max[j]) blad_max[j] = blad[j][i]; //blad bezwzgledny x=x+h; } t=t+dt; x=0; } zapis_do_pliku_m ("blad_bezwzgledny_dla_laasonen.txt", blad); zapis_do_pliku_w ("blad_max_dla_laasonen.txt", blad_max); } void metoda::zapis_do_pliku_m (CHAR* nazwa, DOUBLE** macierz) { ofstream file( nazwa ); DOUBLE x = 0, t = 0; file<<nazwa<<", h = "<<h<<", dt = "<<dt<<endl<<endl; FOR (INT j = 0; j < ile_t; j++) { file<<"Dla poziomu czasowego = "<<t<<endl<<endl; FOR (INT i = 0; i < ile_x ; i++) { file<<"wartosc przestrzenna: "<<x<<",\t U(x,t) = "<<macierz[j][i]<<endl; x=x+h; } file<<"\n"<<endl; t = t+dt; x = 0; } file.close(); } void metoda::zapis_do_pliku_w (CHAR* nazwa, DOUBLE* wektor) { ofstream file( nazwa ); DOUBLE t = 0; file<<nazwa<<", dt = "<<dt<<endl<<endl; FOR (INT j = 0; j < ile_t; j++) { file<<"Dla poziomu czasowego = "<<t<<"\t blad_max = "<<wektor[j]<<endl; t = t+dt; } file.close(); } INT main () { metoda nowa; nowa.rozw_analityczne (); nowa.rozw_laasonen (); nowa.f_blad (); cout<<"Program rozwiazuje rownanie rozniczkowe czastkowe opisujace trasport ciepla"<<endl; cout<<"METODA POSREDNIA LASSONEN ORAZ ALGORYTMEM THOMASA"<<endl<<endl; cout<<"Wzor rownania: dU(x,t) / dt = D d^2U(x,t) / dx^2"<<endl<<endl; cout<<"Wspolrzedna przestrzenna okreslona na zbiorze: [0; +nieskonczonosci)"<<endl; cout<<"Wspolrzedna czasowa okreslona na zbiorze: [0; t_max]"<<endl; cout<<"Warunek poczatkowy: U(x,0) = 0"<<endl; cout<<"Warunki brzegowe U(0,t) = 1, U(+nieskonczonosci,t) = 0"<<endl<<endl;; cout<<"t_max = 2, D = 1, lambda = 1 \nprzedzial nieskonczonosc zastapiony przedzialem >= od 6*sqrt(D*t_max)"<<endl; cout<<"Przedzial ten wynosi: "<<x_max<<endl<<endl; cout<<"Ilosc podzialow na osi x opisujacej przestrzej wybrano: "<<nowa.ile_x<<endl; cout<<"Wyliczony krok na osi przestrzennej h = "<<nowa.h<<endl; cout<<"Krok na osi czasowej wyliczono ze wzoru: lambda = D*dt/h^2 i wyniosl dt = "<<nowa.dt<<endl; cout<<"Ilosc podzialow na osi t opisujacej czas: "<<nowa.ile_t<<endl; cout<<"\n\n\nPoszczegolne wyniki dla kolejnych poziomow czasowych zapisano w plikach oraz blad bezwzgledy"; cout<<" w porownaniu do rozwiazania analitycznego"<<endl; system ("pause"); }