Résolution de l'equation de chaleur

Résolu/Fermé
mirinda - 27 avril 2008 à 02:33
 dhekra88 - 19 déc. 2012 à 11:39
Bonjour,je suis un nouveau membre ,j'espère que vous m'accepterez parmi vous .j'ai un projet en analyse numérique qui porte sur la résolution de l'equation de chaleur et sa programmation en scilab en utilisant la méthode de différence finie pour l'eqution en dimension 1 puis 2 ,mais malheureusement je n'ai pas compris cette méthode, c'est pour cela que je demande votre aide pour résoudre cette equation en traitant tous les cas particuliers et générals ,j'attend vos réponce ,merci d'avance.

53 réponses

mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
2 mai 2008 à 15:06
salut ,les liens ce sont utile mais malheureusement on n'a pas encore fait des cours concernant ceci et moi je cherche des trucs trés détaillé.voilà le programme en scilab que j'ai fait mais je n'ai pas pu le terminer :

lg = 10. ; // lg = demie longueur du domaine
dx = 0.05 ; // dx = pas d'espace
nx = lg/dx ; // nx = demi nombre de mailles
cfl = 0.4 ; // cfl
dt = dx*dx*cfl ; // dt = pas de temps
nt = 500 ; // nt = nombre de pas de temps effectues
//
// initialisation
//
x=zeros(1,2*nx+1) ;
u0=zeros(1,2*nx+1) ;
for i=1:2*nx+1
x(i) = (i-nx-1)*dx ;
u0(i) = max(0.,1.-x(i)**2) ;
end
u=u0 ;
up=u0 ;
um=u0 ;
uexacte=u0 ;

tics=[4,10,4,10];
plotframe([-lg,-0.1,lg,1.1],tics);
plot2d(x,u0,1,"000")
xtitle ('donnee initiale' ,' ',' ') ;


halt() ;
////////////////////////////////////////////////////////////////
// schema explicite: cfl=0.4
////////////////////////////////////////////////////////////////
//
for n=1:nt
//
up = sh('+1',u) ;
um = sh('-1',u) ;

u=u + dt/(dx*dx)*(up+um-2*u) ;

plotframe([-lg,-0.1,lg,1.1],tics);
plot2d(x,u,[1,1],"100","schema explicite")
plot2d(x,u0,[2,2],"100","donnee initiale")
xtitle('schema explicite, cfl=0.4',' ',' ');
xpause(1000) ;
end

avec sh une fct qui décale un tableau d'element de +1 ou -1 .
ici j'ai traité le cas de shéma explicite , est ce que vous pouvez m'aider à le terminer ?merci et by.
17
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
30 avril 2008 à 01:13
salut ,merci beaucoup ,mais pour l'equation de chaleur ca dépend de t et x l'equation est Ut+aUx,x=0 (ux,x est la dérivés seconde par rapport à x et a>0)avec u(x,0)=u°(x) et u(0,t)-u(l,t)=0 pour l appartient à l'inetrvalle ]0,1[
voila ce que j'ai fait
(Un+1j-Ujn)/Δt – a Uj+1n -2Ujn +Uj-1n/Δx2 =0 ;
U0n+1=UM+1n+1=0 (n+1 c'est la puissace et M+1 est l' indice ...)
D’où Ujn+1=Uj n+ a Δt/Δx2(Uj+1n – 2Ujn + Uj-1n) (2 c'est la puissance)
D’où on obtient le système Un+1=BUn

C’est ici que je n’ai pas compris comment obtenir le schéma explicite et implicite et aussi comment étudier et mettre en evidence la convergence et la stabilité de chaque schéma
5
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
2 juin 2008 à 20:58
salut Sacabouffe
merci pour la remarque , mais ça ne marche pas meme si j'ai fait
D=inv(A);
W=D*B;
2
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
29 avril 2008 à 14:47
merci beaucoup ,j'ai compris un petit peu la méthode mais je n'ai pas compris comment arrivé au système linéaire qu'il faut résoudre et aussi la différence entre la cas explicite et implicite ,c'est pas bien expliqué dans les sites .concernant la programmation ,vous aver des propositions?
1

Vous n’avez pas trouvé la réponse que vous recherchez ?

Posez votre question
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
3 mai 2008 à 22:19
salut ,merci beaucoup pour ton aide ,le problème c'est que dans le cas du shéma explicite j'ai programmé juste la solution initiale ,il fallait que je programme aussi la solution exacte en prenant en considération la convergence et la stabilité et aussi les erreurs mais je n'ai pas pu le faire.
Concernant le shéma implicite ,voilà ce que j'ai fait

////////////////////////////////////////////////////////////////
// schema implicite: cfl=2.
////////////////////////////////////////////////////////////////
cfl = 2. ;
dt = dx*dx*cfl ;
u = u0 ;
nt=200 ;
mat=zeros(2*nx+1,2*nx+1) ;

for i=2:2*nx
mat(i,i) = 1. + 2*dt/(dx*dx) ;
mat(i,i+1) = -dt/(dx*dx) ;
mat(i,i-1) = -dt/(dx*dx) ;
end
mat(1,1) = 1. + 2*dt/(dx*dx) ;
mat(1,2) = -dt/(dx*dx) ;
mat(2*nx+1,2*nx) = -dt/(dx*dx);
mat(2*nx+1,2*nx+1) = 1. + 2*dt/(dx*dx) ;
smat = sparse(mat) ;
spcho = chfact(smat) ; // factorisation de Cholesky

//
for n=1:nt
//
u = chsolve(spcho,u) ; // resolution du systeme lineaire

plotframe([-lg,-0.8,lg,1.6],tics);
plot2d(x,u,[1,1],"100","schema implicite 8")
plot2d(x,u0,[2,2],"100","donnee initiale 8")
xtitle('schema implicite 8, cfl=2.',' ',' ');

end
//
j'ai utilisé la factorisation de cholesky ,mais dans ce cas aussi j'ai traité que la solution initiale .
Dans les deux programmes ca ne marche pas bien ,je ne sais pas pourquoi .j'attend tes propositions avec patience,merci.
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
7 mai 2008 à 20:10
salut Sacabouffe
merci beaucoup , ah oui je me trompe toujours j'écris solution initiale au lieu de solution numérique mais tu as compris maintenant ce que j'ai voulu dire .
Puisqu' elle me reste que la solution au dernier pas de temps il faut que je fait la comparaison juste au temps final n'est ce pas ?tu m'a dit "Par exemple , avec un Dirac de température en (0,0) la solution exacte de l'edp est très simple."
mais je n'ai pas étudié " la distribution de Dirac" ni "Transformée de Fourier ....." donc je ne peux pas l'utiliser.
je n'ai pas compris ce que tu veux dire par "sauve la température dans une matrice. Les lignes pour l'espace et les colonnes pour le temps ".merci
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
10 mai 2008 à 00:55
salut sacabouffe
Merci beaucoup pour tous ces informations , oui je suis d'accord avec toi c'est maintenant que je veux modifier le programme pour qu'il utilise u(0,t)=u(1,t)=0 càd les conditions aux limites de direchlet et puis je veux modifier le programme en utilisant les conditions aux limites de Neuman .je veux que ca sera sur l'intervalle [0,1] ,et je veux savoir qu'elles sont les modifications que je doit faire pour qe ca marche bien en faisant apparaitre la stabilité quand on change la condition cfl, la convergence et la consistence .J'attend votre réponce ,merci .
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
12 mai 2008 à 15:08
salut
Voila le nouveau problème détaillé :

Ut – a uxx =0 pour (x,t) dans ]0,1[*R+*
U(x,0)=u0(x) pour x dans ]0,1[
U(0,t)=u(1,t)=0 pour t dans R+* condition de Dirichlet

on pose Uj0=U0(xj ) ; condition initiale
xj =j dx ;
tn =n dt ;
dx=1/(M+1) ; avec M le nombre de nœud de maillage sur ]0,1[

le schéma explicite est le suivant :
Ujn+1 = Ujn + a dt/dx2 (Uj+1n – 2 Ujn + Uj-1n)
U0n+1 = UM+1n+1=0 conditions aux bords
merci .
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
14 mai 2008 à 02:08
Salut
J'ai un peu de mal à comprendre, un coup ta longueur c'est l, ton programme est fait pour x Є [-l/2,l/2], et un coup c'est 1.

Pour le cas périodique, il y a quand même un petit truc à corriger pour le cas explicite et le cas implicite, j'avais pas fait gaffe...
Il faut prendre qu'un seule des extrémités comme ddl, vu que les températures sont toujours égales.
Du coup, pour les deux schémas, j'écrirais plutôt ça:
x=zeros(1,2*nx) ;
u0=zeros(1,2*nx) ;
for i=1:2*nx
x(i) = (i-nx-1)*dx ;
u0(i) = max(0.,1.-x(i)**2) ;
end
u=u0 ;
up=u0 ;
um=u0 ;
uexacte=u0 ;

Et pour le schéma implicite, ta matrice est la même sauf qu'elle est de taille 2*nx x 2*nx et que comme je te l'avais fait remarquer, faut pas oublier les valeurs non nulles de M(2*nx,1) et M(1,2*nx).
Et après, tout pareil. Sauf qu'à la fin, si tu veux aussi ta température à l'extrémité droite (qui est la même qu'à l'extrémité gauche), tu complètes le vecteur u en un vecteur uper tel que uper(1:2*nx)=u(1:2*nx) et uper(2*nx+1)=u(1).

Pour le problème de Dircihlet, c'est presque pareil. Cette fois, les deux extrémités ne sont pas des ddl. Tu travailles donc avec un vecteur de taille 2nx -1 pour la température et t'appliques tes schémas. Pour le schéma implicite, la matrice est quasiment la même, c'est en fait celle que t'avais écrite au départ, elle est de taille 2*nx-1 x 2*nx-1 et cette fois donc, les valeursM(2*nx-1,1) et M(1,2*nx-1) sont bien nulles.

Pour le problème de Neumann, de nouveau, il y a pas grand chose qui change. Les deux extrémités ne sont pas des ddl. Vu que la dérivée est nulle, la température aux extrémités peut être approchée par la température de la barre au premier pas d'espace pour l'extrémité gauche et au dernier pas d'espace pour l'extrémité droite.
Du coup, pour le schéma implicite, seules les valeurs M(1,1) et M(2*nx-1,2*nx-1) changeront.
A plus
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
17 mai 2008 à 19:45
salut Sacabouffe

merci pour tous ces informations ,le problème maintenant c'est que le programme concernant le schéma explicite est lent ,il prend beaucoup de temps pour donner le graphe final et je ne sais pas pourquoi ,contrairement au schéma implicite qui donne le graphe rapidement.Malgré ca je peux dire maintenant que les programmes marchent bien
il me reste maintenant à comparer la solution numérique obtenue avec la solution exacte premièrement pour le schéma explicite,je veux savoir quelles sont les modification que je dois faire pour faire ca , je veux savoir aussi est ce qu'il ya d'autre manière pour programmer le schéma explicite ,merci d'avance.
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
18 mai 2008 à 01:31
Salut
Pour une condition donnée (Dirichlet, Neumann ou périodique) et à pas de temps et pas d'espace égaux, le schéma explicite et le schéma implicite mettent le même temps.
Si t'as fixé le même pas d'espace pour le schéma explicite et le schéma implicite, le pas de temps est limité par la condition CFL pour le schéma explicite. Si t'as mis un pas de temps énorme pour le schéma implicite, c'est normal que tu vois ta solution évoluer plus rapidement vers la solution stationnaire. Je sais pas si c'était vraiment ça ta question...
Pour la solution exacte, il te faut la calculer d'abord.
A plus
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
22 mai 2008 à 14:49
salut Sacabouffe
merci beaucoup pour ton grand aide .
Grace à ton aide j'ai pu términer la première partie de mon travail, maintenant je dois passer à la deuxième partie qui concerne la résolution de l'équation de chaleur en dimension 2 c'est pour cela je te demande si tu peux me donner des informations sur cette dernière et sur sa programmation .merci
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
26 mai 2008 à 02:30
Salut
C'est la même chose en un peu plus compliqué. T'auras 3 indices, un pour la discrétisation selon l'axe (Ox), un pour la discrétisation selon l'axe (Oy) et un pour la discrétisation temporelle. Si on note i l'indice de discrétisation selon (Ox), j l'indice de discrétisation selon (Oy), k l'indice de discrétisation temporelle ainsi que Δx, Δy, Δt, les pas respectifs de discrétisation, pour le schéma explicite t'as un truc du genre:

(u(i,j,k+1)-u(i,j,k))/(Δt)=D[(u(i+1,j,k)-2u(i,j,k)+u(i-1,j,k))/(2Δx)+(u(i,j+1,k)-2u(i,j,k)+u(i,j-1,k))/(2Δy)]

D est la diffusivité thermique.
A plus
0
bonjour mes amis.
s'il vous plait,vous pouvez me donner le programme de la solution de l'equation de chaleur que vous avez trouvé? dans le cas implicite + explicite.s'il vous plait c urgent.
merci d'avance.
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
30 mai 2008 à 10:55
Salut
Apprends à lire. Et quand t'auras appris, lis ça aussi:
http://www.commentcamarche.net/faq/sujet 10925 demander de l aide pour vos exercices sur ccm
A plus
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
31 mai 2008 à 01:55
salut Sacabouffe
merci pour tout, je revient à la pemière partie du travail càd l'équation de la chaleur en dimension 1
je doit la programmer en dimension 1 par un schéma de Crank Nicolson (teta schéma implicite) qui s'écrit:

uj,n+1 - µ*teta[uj+1,n+1 - 2*uj,n+1 + uj-1,n+1]=µ*(1-teta)[uj+1,n - 2*uj,n + uj-1,n] + uj,n

avec µ=a*(dt/dx²) et teta dans [0,1]
càd
A* Un+1=B*Un
je veux savoir comment je peux résoudre ce système pour obtenir la solution u sur un intervalle [0,1], dans le cas de l'implicite j'avais A*Un+1=Un pour le résoudre j'ai utilisé la factorisation de cholesky mais maintenant j'ai deux matrices,comment faire.
voilà ce que j'ai essayé mais je n'ai pas pu avancer :

A=zeros(2*nx+1,2*nx+1) ;
B=zeros(2*nx+1,2*nx+1) ;

for i=2:2*nx
A(i,i) = 1. + dt/(dx*dx) ;
A(i,i+1) = -1/2*dt/(dx*dx) ;
A(i,i-1) = -1/2*dt/(dx*dx) ;

end
A(1,1) = 1. + dt/(dx*dx) ;
A(1,2) = -dt/(dx*dx) ;
A(2*nx+1,2*nx) = -1/2*dt/(dx*dx);
A(2*nx+1,2*nx+1) = 1. + dt/(dx*dx) ;

for i=2:2*nx

B(i,i) = -dt/(dx*dx) ;
B(i,i+1) = 1/2*dt/(dx*dx) ;
B(i,i-1) = 1/2*dt/(dx*dx) ;
end

B(1,1) = -dt/(dx*dx) ;
B(1,2) = 1/2*dt/(dx*dx) ;
B(2*nx+1,2*nx) = 1/2*dt/(dx*dx);
B(2*nx+1,2*nx+1) = -dt/(dx*dx) ;
D=inv(B);
W=A*D;
Est ce que tu peux m'aider à terminer ce programme?(avec toujours les memes conditions qu'avant et conditions de direchlet)
Merci d'avance.
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
31 mai 2008 à 19:14
salut
Dans la programmation j'ai posé teta=1/2 .
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
2 juin 2008 à 18:09
Salut mirinda
J'ai pas vérifié entièrement ton code mais il y a un truc bizarre à la fin:
D=inv(B);
W=A*D;

Si tu dois résoudre A* U(n+1)=B*U(n) c'est plutôt un truc comme ça que tu dois écrire:
D=inv(A);
W=D*B;

Ensuite tu calcules U(n+1) en faisant U(n+1)=W*U(n).
Dis-moi si ça marche comme ça. Si c'est pas le cas, je regarderai plus attentivement ce que t'as écrit.
A plus
0
Sacabouffe Messages postés 9427 Date d'inscription dimanche 19 août 2007 Statut Membre Dernière intervention 29 mai 2009 1 832
2 juin 2008 à 21:01
Ah... bon ben je vais regarder plus en détails ton programme alors.
Je te réponds ce soir ou demain.
0
mirinda Messages postés 51 Date d'inscription dimanche 27 avril 2008 Statut Membre Dernière intervention 4 mai 2009 38
3 juin 2008 à 01:41
D'accord, merci .
0