Bonjour,
J'ai besoin de résoudre une équation différentielle homogène autonome de second ordre :
a d²u/dx²+b du/dx+c u=0 pour x ϵ [0 , L], les conditions de bord sont u(0) = U0 et u(L) = UL.
Ceci revient à résoudre l'équation matricielle : A * U = B (voir fichier ci-joint schéma-numérique.docx )
Voici mon code :
Pour le traçage des courbes, j'utilise Gnuplot (solution analytique en vert, solution numérique en violet).Code:#include <stdio.h> #include <stdlib.h> #include <math.h> /* Triangularisation d'un système tridiagonale par la méthode du pivot de Gauss */ void Triangularisation (double **A, double *B, int N) { int i; double p,q; for( i = 0 ; i<N-1 ; i++ ) { p= A[i][i]; q= A[i+1][i]; A[i+1][i]=0; A[i+1][i+1]= A[i+1][i+1]- (q/p)* A[i][i+1]; B[i+1]=B[i+1]-(q/p)*B[i]; } } /* Résolution d'un système triangulaire supérieure par la méthode de remontée */ double* SolutionSyst (double **A, double *B,int N) { int i; double *U; U = calloc( N, sizeof(double)); if( U == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } U[N-1]= B[N-1]/A[N-1][N-1]; for( i = N-2 ; i>=0 ; i-- ) { U[i]= (B[i]-A[i][i+1]*U[i+1])/A[i][i]; } return (U); } /* Résolution de l'équation différentielle par la méthode des différences finies */ double* MethodesFinies (double a, double b, double c, double L, double U0, double UL, int N, double h) { int i; double *B,*U; double **A; double e,f,g; /* Calcul de e, f et g */ e = a/pow(h,2)+b/h; f = c - (2*a)/pow(h,2)-b/h; g = a/pow(h,2); /* Allocation de la matrice A et des vecteurs U et B */ A = malloc( N * sizeof(double*)); if( A == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } for( i = 0 ; i < N ; i++ ) { A[i] = calloc (N, sizeof(double)); if( A[i] == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } } U = calloc( N, sizeof(double)); if( U == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } B = calloc( N, sizeof(double)); if( B == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } /* Remplissage de la matrice tridiagonale A de taille N*N et du vecteur colonne B de taille N */ A[0][0]= f; A[0][1]= e; A[N-1][N-2]= g; A[N-1][N-1]=f; for( i = 1 ; i<N-1 ; i++ ) { A[i][i-1] = g; A[i][i] = f; A[i][i+1]= e; } B[0]= -g*U0; B[N-1]= -e*UL; /* Triangularisation du système */ Triangularisation (A,B,N); /* Résolution du système */ U = SolutionSyst (A, B,N); return (U); } /* Résolution analytique de l'équation différentielle avec un Delta > 0 */ double SolAnal1 (double x, double a, double b, double L, double U0, double UL,double D) { double r1, r2, A, B, y; r1 = (-b+sqrt(D))/(2*a); r2 = (-b-sqrt(D))/(2*a); B = (UL-U0*exp(r1*L))/(exp(r2*L)-exp(r1*L)); A = U0-B; y= A*exp(r1*x)+B*exp(r2*x); return (y); } /* Résolution analytique de l'équation différentielle avec un Delta = 0 */ double SolAnal2 (double x, double a, double b, double L, double U0, double UL) { double r, A, B, y; r = -b /(2*a); B = U0; A = ((UL/exp (r*L))- U0)/L; y = (A*x+B)*exp(r*x); return (y); } int main(int argc, char const *argv[]) { int N,i; double a,b,c,h,L,U0,UL,D,y; double *U; FILE* f = NULL; /* Sasie de a, b, c, L, U0, UL et N par l'utilisateur */ printf("Saisissez la valeur de a \n"); scanf("%lf",&a); printf("Saisissez la valeur de b \n"); scanf("%lf",&b); printf("Saisissez la valeur de c \n"); scanf("%lf",&c); printf("Saisissez la valeur de L \n"); scanf("%lf",&L); printf("Saisissez la valeur de U0 \n"); scanf("%lf",&U0); printf("Saisissez la valeur de UL \n"); scanf("%lf",&UL); printf("Saisissez le nombre des points de maillage \n"); scanf("%d", &N); /* Calcul du pas */ h= L/(N+1); /* Allocation et initialisation du vecteur U */ U = calloc( N, sizeof(double)); if( U == NULL ) { printf("Allocation impossible"); exit(EXIT_FAILURE); } /* Résolution de l'équation avec la méthode des différences finies */ U = MethodesFinies (a, b, c, L, U0, UL, N,h); /* Résolution analytique de l'équation et écriture des images numériques, exactes et de l'erreur pour chaque noeud dans le fichier "fichier.txt" */ f = fopen ("fichier.txt","w"); if (f!= NULL) { fprintf(f,"0 \t\t %lf \t%lf \t 0 \n", U0, U0); D = pow(b,2)-4*a*c ; if (D>0) for(i=1 ; i<N+1 ; i++) { fprintf(f,"%lf \t %lf \t",i*h, U[i-1]); y = SolAnal1 (i*h, a, b, L, U0, UL, D); fprintf (f,"%lf \t %lf\n",y,fabs(y-U[i-1])); } else for(i=1 ; i<N+1 ; i++) { fprintf(f,"%lf \t %lf \t",i*h, U[i-1]); y = SolAnal2 (i*h, a, b, L, U0, UL); fprintf (f,"%lf \t %lf\n",y,fabs(y-U[i-1])); } fprintf(f,"%lf \t %lf \t%lf \t 0 \n", L, UL, UL); fclose (f); } return 0; }
Pour le cas de a = -1, b = 0.1, c = 1, L = 1, U0 = 2, UL = 1, N = 10000, j'obtiens 2 courbes non confondues pour les solutions numérique et analytique (sachant que pour N = 10, 100, 1000 j'obtenais 2 courbes confondues).
1.png
Je pense que l'écart observé provient des erreurs d'arrondi et je me demandais s'il y'aurait un moyen de pallier à ces erreurs.
Merci de votre aide
-----