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 :

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 traçage des courbes, j'utilise Gnuplot (solution analytique en vert, solution numérique en violet).

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