Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <iostream>
- #include <stdlib.h>
- #include <cstdlib>
- #include <fstream>
- #include <String>
- #include <Math.h>
- #include <sys/timeb.h>
- using namespace std;
- fstream file,file2,file3,file4,file5,file6,file7,file8,file9,file10;
- //deklaracja stałych:
- const double D = 1.0;
- const double lambda_bez = 1.0; //dla metod bezposredniej
- const double lam_pos = 1.0; //dla metod posrednich
- const double h = 0.1; //krok dla x
- const double dt = (lam_pos*h*h)/D; //krok dla t
- const double PI=3.14159265359;
- const double x_max=1; //maksymalny przedzial x
- const double x_min=0; //minimalny przedzial x
- const double t_max=0.5; //maksymalny przedzial x
- const double t_min=0; //minimalny przedzial x
- //wymiary macierzy do obliczeń numerycznych:
- //const int N = (int)((2*x_max/h)+1);//liczba kolumn (os x)
- //const int N = 40; //liczba kolumn(x)
- //const int M = (int)((t_max/dt)+1);//liczba wierszy (os t)
- const int M=(t_max-t_min)/dt;
- const int N=(x_max-x_min)/h;
- //**************************************************************************************************
- //tworzy nową macierz MxN
- double **macierz(int m, int n)
- {
- double **a;
- a=new double *[m];
- for(int i=0; i<m; i++)
- a[i]=new double[n];
- return a;
- }
- //**************************************************************************************************
- //usuwa macierz
- void usun_macierz(double **a, int m)
- {
- for(int i=m-1; i>=0; i--)
- delete []a[i];
- delete []a;
- }
- //**************************************************************************************************
- //zapisywanie macierzy do pliku
- void zapisz(double **macierz, fstream &plik)
- {
- //plik.setf(ios::scientific, ios::floatfield);
- //plik.precision(4);
- for(int j = 0; j<M; j++){
- for(int i=0;i<N;i++){
- plik <<i*h<<" "<<j*dt<<" " <<macierz[j][i] << endl;
- }
- }
- plik.close();
- }
- //**************************************************************************************************
- void zapiszW(double *wek, fstream &plik)
- {
- for(int j = 0; j<M; j++){
- plik << j*h << wek[j] << "\n";
- }
- plik.close();
- }
- /*********************************************************************************************/
- /*metoda obliczająca czas wykonywania obliczen*/
- //**************************************************************************************************
- //Metoda rozwiązująca analitycznie równanie rózniczkowe
- void analitycznie(double **A, double *x, double *t)
- {
- //warunek początkowy
- for(int i=0;i<M;i++){
- A[0][i]=0;
- }
- //warunek brzegowy
- for(int i=0;i<N;i++){
- A[i][M-1]=0;
- A[i][0]=0;
- }
- //wypełnienie całej macierzy:
- /*for(int i =1; i<M; i++)
- {
- for(int j =1; j<N-1; j++)
- {
- A[j][i] = (1.0 - exp(-1.0*PI*PI*D*j*dt))*sin(PI*i*h);
- }
- }
- */
- double y=0;
- for (int n = 1;n<N-1;n++)
- {
- for (int m = 1;m < M; m++)
- {
- A[m][n] = (1.0 - exp(-1.0*PI*PI*D*m*dt))*sin(PI*n*h);
- }
- }
- //zapisanie rozwiązania do pliku:
- zapisz(A, file);
- }
- //**************************************************************************************************
- //klasyczna metoda bezposrednia
- //**************************************************************************************************
- void blad(double *BLAD, double **KMB, double **A, fstream &plik)
- {
- double temp;
- for(int i =0; i<M; i++)
- {
- double max = 0.0;
- for(int j =0; j<N-1; j++){
- temp = fabs(KMB[i][j] - A[i][j]);
- if(temp>max){
- max = temp;
- }
- }
- BLAD[i] = max;
- }
- zapiszW(BLAD, plik);
- }
- //**************************************************************************************************
- double norma(double *dane){
- double w = dane[0];
- for( register int i = 1; i < N; i++ ){
- if( w < dane[i] )
- w = dane[i] ;
- }
- return w;
- }
- //**************************************************************************************************
- //funkcja dokonująca częściowego wyboru elementów podstawowych
- void szukaj_zamien(double **A, int *indeksy, int start){
- int max = start;
- for(register int i = start+1; i < M; i++ )
- if( fabs( A[indeksy[i]][start] ) > A[indeksy[max]][start] )
- max = i;
- //sprawdzanie czy udało się znaleść odpowiedni element
- if( A[indeksy[max]][start] == 0 ){
- printf("Bladna macierz A\n");
- exit(1);
- }
- indeksy[start] = max;
- indeksy[max] = start;
- }
- //**************************************************************************************************
- //funkcja dokonująca dekompozycji LU
- void rozklad_LU(double **LU, double **A, int *indeksy){
- double wsp;
- for(int i = 0; i < N; i++ ){
- for(int j = 0; j < N; j++ ){
- LU[indeksy[i]][j] = A[indeksy[i]][j];
- }
- }
- for( int i = 0; i < N; i++ ){
- for(int j = i+1; j < N; j++ ){
- if(LU[indeksy[i]][i] == 0 ) szukaj_zamien(LU,indeksy,i);
- wsp = LU[indeksy[j]][i]/LU[indeksy[i]][i];
- for(int k = i; k < M; k++ ){
- LU[indeksy[j]][k] = LU[indeksy[j]][k] - LU[indeksy[i]][k]*wsp;
- if( i < j ){
- LU[indeksy[j]][i] = wsp;
- }
- }
- }
- }
- }
- //**************************************************************************************************
- //funkcja dokonująca obliczeń wektora y z równania Ly=b
- void obliczanie_y(double **LU, double *Y, double *B, int *indeksy, const int n){
- double S;
- for(register int i = 0; i < n; i++ ){
- S = 0;
- for(register int j = 0; j < i; j++ ){
- S += Y[j]*LU[indeksy[i]][j];
- }
- Y[i] = B[indeksy[i]] - S;
- }
- }
- //**************************************************************************************************
- //funkcja dokonująca obliczeń wektora x z równania Ux=y
- void obliczanie_x(double **LU, double *X, double *Y, int *indeksy, const int n){
- double S;
- register int j;
- for( register int i = n-1; i >= 0; i-- ){
- S = 0;
- for( j = n-1; j > i; j-- )
- S += X[j]*LU[indeksy[i]][j];
- X[i] = (Y[i] - S)/LU[indeksy[i]][i];
- }
- }
- //**************************************************************************************************
- //funkcja inicjalizująca tablicę indeksów potrzebna dla powyższych funkcji.
- void init_kolumny(int *indeksy, const int n ){
- for( register int i = 0; i < n; i++ )
- indeksy[i] = i;
- }
- //**************************************************************************************************
- void crank_nicolson_LU(double **A,double *x, double *t, double lam_pos){
- register int i,j,k;
- double *b,*z;
- double **LU = new double * [N];
- double **WSP = new double * [N];
- for( i = 0 ; i < N ; i++ ){
- LU[i] = new double [N];
- WSP[i] = new double [N];
- }
- double *Y = new double [N];
- int *indeksy = new int [N];
- b=new double[N];
- z=new double[N];
- for(j=0;j<N;j++)
- A[0][j]=0;
- for(i=0;i<M;i++){
- A[i][0]=0;
- A[i][N-1]=0;
- }
- for( i = 0; i < N; i++ )
- for( j =0; j< N; j++ ){
- if( j == i-1 || j == i+1){
- WSP[i][j] = -lam_pos/2.0;
- }else if( j == i ){
- WSP[i][j] = 1 + lam_pos;
- }else{
- WSP[i][j] = 0;
- }
- }
- init_kolumny(indeksy,N);
- rozklad_LU(LU,WSP,indeksy);
- for( i=1;i<M;i++){
- b[0]=b[N-1]=0;
- for(j=1;j<N-1;j++)
- //b[j]=(lambda/2.0)*A[i][j-1]+(1-lambda)*A[i][j]+(lambda/2.0)*A[i][j+1];
- b[j]=((lam_pos/2.0)*A[i-1][j-1]+(1.0-lam_pos)*A[i-1][j]+(lam_pos/2.0)*A[i-1][j+1]);
- for(k=0;k<N-1;k++)
- b[k]=b[k]+(D*dt*PI*PI*sin(PI*k*h));
- obliczanie_y( LU,Y, b,indeksy,N);
- obliczanie_x( LU,z, Y,indeksy,N);
- for(j=1;j<N-1;j++)
- A[i][j]=z[j];
- }
- zapisz(A, file4);
- }x`x
- //**************************************************************************************************
- //gauss Seidel
- void GS(double **A, double *x, double *b){
- double tmp;
- double *xp = new double[N];
- double *tmp_w = new double[N];
- for( int i = 0; i < N; i++ )
- xp[i] =x[i] + 1;
- do{
- for( int i = 0; i < N; i++ ){
- tmp = 0;
- for( int j = 0; j < i-1; j++ ){
- tmp += A[j][i] * x[j];
- }
- for( int j = i+1; j < N; j++ ){
- tmp += A[j][i] * x[j];
- }
- x[i] = b[i] - tmp;
- x[i] /= A[i][i];
- }
- for( int i = 0; i < N; i++ ){
- tmp_w[i] = fabs(x[i] - xp[i]);
- xp[i] = x[i];
- }
- }while( norma(tmp_w) > 0.00000000001 );
- }
- void crank_nicolson_GS(double **A, double *x, double *t, double lam_pos){
- register int i,j,k;
- double *b,*z;
- double **LU = new double * [N];
- double **WSP = new double * [N];
- for( i = 0 ; i < N ; i++ ){
- LU[i] = new double [N];
- WSP[i] = new double [N];
- }
- b=new double[N];
- z=new double[N];
- for(j=0;j<N;j++){
- A[0][j]=0;
- z[j] = 1;
- }
- for(i=0;i<M;i++){
- A[i][0]=0;
- A[i][N-1]=0;
- }
- for( i = 0; i < N; i++ )
- for( j =0; j< N; j++ ){
- if( j == i-1 || j == i+1){
- WSP[i][j] = -lam_pos/2.0;
- }else if( j == i ){
- WSP[i][j] = 1 + lam_pos;
- }else{
- WSP[i][j] = 0;
- }
- }
- for( i=1;i<M-1;i++){
- b[0]=b[N-1]=0;
- for(j=1;j<N-1;j++)
- //b[j]=(lambda/2.0)*A[i][j-1]+(1-lambda)*A[i][j]+(lambda/2.0)*A[i][j+1];
- b[j]=((lam_pos/2.0)*A[i-1][j-1]+(1.0-lam_pos)*A[i-1][j]+(lam_pos/2.0)*A[i-1][j+1]);
- for(k=0;k<N-1;k++)
- b[k]=b[k]+(D*dt*PI*PI*sin(PI*k*h));
- GS(WSP,z,b);
- for(j=1;j<N-1;j++)
- A[i][j]=z[j];
- }
- zapisz(A, file6);
- }
- //**************************************************************************************************
- int main(){
- file.open("Analit.txt",ios::out);
- file6.open("CN-GS.txt",ios::out);
- file5.open("Bledy_CN-LU.txt",ios::out);
- file4.open("CN-LU.txt",ios::out);
- file7.open("Bledy_CN-GS.txt",ios::out);
- double *x = new double[N];
- double *t = new double[M];
- for(int i =1; i<N; i++)
- x[i] = x_min + (i*h);
- for(int i =1; i<M; i++)
- t[i] = t_min + (i*dt);
- double **A = macierz(M, N); //macierz do rozwiązania analitycznego
- double **KMB = macierz(M,N);
- double *BLAD = new double[M];
- double **CN_LU = macierz(M,N);
- double **CN_GS = macierz(M,N);
- analitycznie(A,x,t);
- crank_nicolson_GS(CN_GS,x,t,lambda_bez);
- blad(BLAD,CN_GS,A,file7);
- crank_nicolson_LU(CN_LU,x,t,lambda_bez);
- blad(BLAD,CN_LU,A,file5);
- return 0;
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement