Equation de la chaleur en 2D. Problème avec mon logiciel.

Marsh Posté le 30-05-2007 à 13:14:46

Voilà dans le cadre d'un projet avec mon école, j'ai à programmer une méthode de résoluton numérique de l'équation de la chaleur en 2D. VOilà la sujet :

Voilà mon code :

Code :

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int N;
//---------------------------------------------
// Structure et fonctions pour gérer des matrices
//---------------------------------------------
typedef struct _Matrix
{
float * data;
int line;
int col;
} Matrix;
void AllocMatrix(Matrix * matrix, int l, int c)
{
matrix->line = l;
matrix->col = c;
matrix->data = (float *)malloc(sizeof(float) * l * c);
}
void FreeMatrix(Matrix * matrix)
{
free(matrix->data);
}
void MatrixSet(Matrix * matrix, int i, int j, float value)
{
matrix->data[j * matrix->line + i] = value;
}
float MatrixGet(Matrix * matrix, int i, int j)
{
return matrix->data[j * matrix->line + i];
}
//-----------------------------------------------------
//Fonction Choleski
//-----------------------------------------------------
Choleski (Matrix A, int N, Matrix H, Matrix Out)
{
int I=1,J=1,K,L,M;
Matrix B;
float S;
float X[N+1], Y[N+1];
int Bool=1;
AllocMatrix(&B, N, N);
//Affichage de Aj dans Choleski
/*printf("Aj dans Choleski \n" );
for (I=1 ; I<=N ; I++)
{
for (K=1 ; K<=N ; K++)
{
printf("%f ",MatrixGet(&A, I, K));
}
printf("\n" );
}
//*/
for (I=0 ; I<=N ; I++)
{
S=0;
for (K=0 ; K<=I-1 ; K++)
{
S=(S+( MatrixGet(&B, I, K)*MatrixGet(&B, I, K) ) );
}
if ( (MatrixGet(&A, I, I) - S )> 0 )
{
MatrixSet(&B, I, I, sqrt(MatrixGet(&A, I, I)-S));
for (J=I+1 ; J<=N ; J++)
{
S=0;
for (K=1 ; K<=I-1 ; K++)
{
S=S+(MatrixGet(&B, I, K)*MatrixGet(&B, J, K));
}
MatrixSet(&B, J, I, (MatrixGet(&A, I, J)-S)/(MatrixGet(&B, I, I)));
}
for (L=1 ; L<=N-1 ; L++)
{
for (M=I+1 ; M<=N ; M++)
{
MatrixSet(&B ,L, M, 0);
}
}
}
else
{
Bool=0;
}
}
if (Bool==0)
{
printf("\n\nResolution impossible\n" );
}
else
printf("\n\nResolution de By=h" );
Y[1]=MatrixGet(&H, 1, 1)/(MatrixGet(&B, 1, 1));
for (I=2 ; I<=N ; I++)
{
S=0;
for (K=1 ; K<=I-1 ; K++)
{
S=S+(MatrixGet(&B, K, I))*Y[K];
}
Y[I]=(MatrixGet(&H, I,1)-S)/MatrixGet(&B, I, I);
}
printf("\n\nResolution de B trans x = y" );
X[N]=Y[N]/(MatrixGet(&B ,N ,N));
for (I=N-1 ; I>= 1 ; I--)
{
S=0;
for (K=I+1 ; K <= N ; K++)
{
S=S+(MatrixGet(&B, K, I)*X[K]);
}
X[I]=(Y[I]-S)/(MatrixGet(&B, I, I));
}
printf("\nX=\n" );
for (I=1 ; I<=N ; I++)
{
printf(" %f\n",X[I]);
MatrixSet(&Out, I, 1, X[I]);
}
FreeMatrix(&B);
}
//------------------------------------------------------------------------
//Fonction principale
//------------------------------------------------------------------------
int main(void)
{
//Déclaration des variables
int n;
//Choix du nombre de colonnes du tableau et donc du pas
printf ("Bienvenue dans le programme de resolution de l'équation de la chaleur \n" );
printf (" en deux dimension d'espace. \n" );
printf("\n" );
printf("Veuillez entrer le nombre de colonnes pour la matrice initiale \n" );
printf(" On rapelle que le pas est egal a 1/(n+1) \n \n" );
printf("n = " );
scanf("%d",&n);
printf("\n" );
//
//Suite de la déclaration des variables
float T=10;
int i,j,l,k,c=15,t_etoile=4,p=200;
float dy,dt,dx,x=0,y=0,t,f;
Matrix u_etoile;
Matrix u_avant;
Matrix u_apres;
Matrix u0, uc1;
Matrix Aj;
Matrix ucl;
Matrix bj;
Matrix Xc;
Matrix Xc2;
Matrix bi;
dx=1.0/(n+1);
dy=1.0/(n+1);
dt= T/p;
float n1= 0.5*dt/(dx*dx);
float n2= 0.5*dt/(dy*dy);
printf("Le pas sera donc de %f \n \n",dx);
printf("Appuyez sur une touche pour continuer \n" );
getchar();
getchar();
// Allocation des divers tableaux
AllocMatrix(&u_etoile, n+1, n+1);
AllocMatrix(&u_avant, n+1, n+1);
AllocMatrix(&u0, n+1, n+1);
AllocMatrix(&ucl, n+1, n+1);
AllocMatrix(&Aj, n+1, n+1);
AllocMatrix(&bj, n+1, 1);
AllocMatrix(&Xc, n+1, 1);
AllocMatrix(&Xc2, n+1, 1);
AllocMatrix(&bi, n+1, 1);
//Fin allocation variable
//-----------------------------------------------------------------
//1ere ETAPE : PASSAGE de u_k à u*
//-----------------------------------------------------------------
// Remplissage de Aj
for (i=1 ; i<n+1 ; i++)
{
for (j=1 ; j<n+1 ; j++)
{
MatrixSet(&Aj, i, j,0);
MatrixSet(&Aj, i, i,1+2*n1);
if (j==(i+1))
MatrixSet(&Aj, i, j,-n1);
if (i==(j+1))
MatrixSet(&Aj, i, j,-n1);
}
}
//
//Début de la boucle principal qui effectuera les divers calculs
for (k=0 ; k<1 ; k=k+1)
{
for (y=dy ; y<1 ; y=y+dy)
{
for (x=dx ; x<1 ; x=x+dx)
{
i = 1 / dx * x;
j = 1 / dy * y;
// Conditions aux limites
MatrixSet(&ucl, 0, 0, 1);
MatrixSet(&ucl, 0, j, 1);
MatrixSet(&ucl, i, 0, 1);
MatrixSet(&ucl, n+1, j, exp(-y));
MatrixSet(&ucl, i, n+1, exp(-x));
MatrixSet(&ucl, n+1, n+1, exp(-1));
MatrixSet(&u0, i, 0, 1);
MatrixSet(&u0, i, n+1, exp(-x));
// Conditions initiales
if (x+y>1)
MatrixSet(&u0, i, j, exp(1-x-y));
else
MatrixSet(&u0, i, j, 1);
// Calcul des valeurs de la fonction f
t=k*(dt);
f=c*exp(-((x-0.8)*(x-0.8)+(y-0.8)*(y-0.8)))*exp(-(t-t_etoile)*(t-t_etoile));
//
//
if(k==0)
{
MatrixSet(&u_avant, i, j, MatrixGet(&u0, i, j));
MatrixSet(&u_avant, i, j-1, MatrixGet(&u0, i, j));
MatrixSet(&u_avant, i, j+1, MatrixGet(&u0, i, j));
}
//
//
MatrixSet(&u_avant, i, 0, MatrixGet(&ucl, i, 0));
MatrixSet(&u_avant, i, n+1, MatrixGet(&ucl, i, n+1));
//
// Remplissage de bj
MatrixSet(&bj, i,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1) - 2*MatrixGet(&u_avant, i,j) ) + 0.5*dt*f );
MatrixSet(&bj, 0,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1))+ (0.5*(dt)*f) + n1*MatrixGet(&ucl, 0, j) );
MatrixSet(&bj, n,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1) - 2*MatrixGet(&u_avant, i,j) ) + (0.5*(dt)*f) + n1*MatrixGet(&ucl, n+1, j) );
//
}
//Appel de la fonction Choleski
Choleski (Aj, n, bj, Xc);
printf("\n" );
//Remplissage de u* avec les résultats de Choleski
for (l=1 ; l<n+1 ; l++)
{
MatrixSet(&u_etoile, l, j, MatrixGet(&Xc, l, 1));
}
}
//////////////////////////////Fin étape 1//////////////////////////
//-----------------------------------------------------------------
//2nde ETAPE : PASSAGE de u* à u_k+1
//-----------------------------------------------------------------
for (y=dy ; y<1 ; y=y+dy)
{
for (x=dx ; x<1 ; x=x+dx)
{
i = 1 / dx * x;
j = 1 / dy * y;
if(k==0)
{
MatrixSet(&u_avant, i, j, MatrixGet(&u_etoile, i, j));
MatrixSet(&u_avant, i, j-1, MatrixGet(&u_etoile, i, j));
MatrixSet(&u_avant, i, j+1, MatrixGet(&u_etoile, i, j));
}
//Calcul des valeurs de la fonction f
t=k*(dt);
f=c*exp(-((x-0.8)*(x-0.8)+(y-0.8)*(y-0.8)))*exp(-(t-t_etoile)*(t-t_etoile));
//
// Remplissage de bi
MatrixSet(&bi, i,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1) - 2*MatrixGet(&u_avant, i,j) ) + 0.5*dt*f );
MatrixSet(&bi, 0,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1))+ (0.5*dt*f) + n1*MatrixGet(&ucl, 0, j) );
MatrixSet(&bi, n,1, MatrixGet(&u_avant, i, j) + n2*( MatrixGet(&u_avant, i, j-1) + MatrixGet(&u_avant, i, j+1) - 2*MatrixGet(&u_avant, i,j) ) + (0.5*(dt)*f) + n1*MatrixGet(&ucl, n+1, j) );
//
//Appel de Choleski pour déterminer u_k+1
//Choleski (Aj, n, bi, Xc2);
//
/* //Remplissage de u_k+1 avec les résultats de Choleski
for (l=1 ; l<n+1 ; l++)
{
MatrixSet(&u_apres, i, l, MatrixGet(&Xc2, l, 1));
}
*/
}
}
/* for (y=dy ; y<1 ; y=y+dy)
{
for (x=dx ; x<1 ; x=x+dx)
{
i = 1 / dx * x;
j = 1 / dy * y;
MatrixSet(&u_avant, i, j, MatrixGet(&u_apres, i, j));
}
}*/
}
//-----------------------------------------------------------------
//Affichage divers pour tests des valeurs
//-----------------------------------------------------------------
//Test de u0
printf("\nAffichage de u0 \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
for(y=dy ; y<1 ; y=y+dy)
{
i= 1 / dx * x;
j= 1 / dy * y;
printf("%f ", MatrixGet(&u0, i, j));
}
printf("\n" );
}
//Test de Ai
printf("\nAffichage de Aj \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
for(y=dy ; y<1 ; y=y+dy)
{
i= 1 / dx * x;
j= 1 / dy * y;
printf("%f ", MatrixGet(&Aj, i, j));
}
printf("\n" );
}
printf("\n" );
//Test de bj
printf("\nAffichage de bj \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
int i= 1 / dx * x;
printf("%f \n", MatrixGet(&bj, i, 1));
}
//Test de Xc
printf("\nAffichage de Xc \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
int i= 1 / dx * x;
printf("%f \n", MatrixGet(&Xc, i, 1));
}
//Test de u*
printf("\nAffichage de u* \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
for(y=dy ; y<1 ; y=y+dy)
{
i= 1 / dx * x;
j= 1 / dy * y;
printf("%f ", MatrixGet(&u_etoile, i, j));
}
printf("\n" );
}
//Test de bi
printf("\nAffichage de bi \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
int i= 1 / dx * x;
printf("%f \n", MatrixGet(&bi, i, 1));
}
//Test de Xc2
printf("\nAffichage de Xc2 \n \n" );
for(x=dx ; x<1 ; x=x+dx)
{
int i= 1 / dx * x;
printf("%f \n", MatrixGet(&Xc2, i, 1));
}
////////////////////////////Fin affichage des tests///////////////////////
getchar();
getchar();
// Libération des divers tableaux
FreeMatrix(&bj);
FreeMatrix(&bi);
FreeMatrix(&Aj);
FreeMatrix(&ucl);
FreeMatrix(&u0);
FreeMatrix(&u_avant);
FreeMatrix(&u_etoile);
FreeMatrix(&Xc);
FreeMatrix(&Xc2);
}
//Fin programme

Et enfin, voilà mon problème :
Quand j'appelle Choleski une première fois, tout va bien, il s'exécute. Par contre quand le programme arrive dans la seconde boucle et appelle Choleski une seconde fois, il plante.
Si je n'appelle pas Choleski dans la première et que je le fais dans la seconde, ça fonctionne.

Donc j'aimerais savoir si vous pouviez me dire pourquoi quand j'appelle une seconde fois ma fonction ça plante. Merci d'avance

Reply

Marsh Posté le 30-05-2007 à 13:14:46

Reply

Marsh Posté le 30-05-2007 à 13:18:44

Reflexe : Utiliser le debugger de ton environnement, voir la ligne ou ca part en couille, analyser la valeur des données en jeu à ce moment et hop.

---------------
Töp of the plöp

Reply

Marsh Posté le 30-05-2007 à 13:20:49

Je sais t rès bien d'où vient le problème : à la ligne 328. Quand je rappelle Choleski ça plante. Si je l'appelle pas une seconde fois tout s'exécute très bien. Et si je ne l'appelle pas à la ligne 280, là il s'exécutera bien à la ligne 328.

Je pense à un problème d'allocation mémoire ou un truc qui me dépasse.

Message édité par StarHunter le 30-05-2007 à 13:21:36

Reply

Marsh Posté le 30-05-2007 à 13:24:09

Ben tu rentres dans cette fonction et tu vois exactement où ca plante.

---------------
Töp of the plöp

Reply

Marsh Posté le 30-05-2007 à 13:27:05

Je fais ça comment ?

Reply

Marsh Posté le 30-05-2007 à 13:27:31

T'as quoi comme environnement de dev ?

Si t'as visual, F11

---------------
Töp of the plöp

Reply

Marsh Posté le 30-05-2007 à 13:28:38

DEV C++.

Reply

Marsh Posté le 30-05-2007 à 13:33:36

Je sais pas sous devcpp, check la doc.

---------------
Töp of the plöp

Reply

Equation de la chaleur en 2D. Problème avec mon logiciel.

Sujets relatifs:

Leave a Replay