Methode iterative

Résolu
Bonjour,
je voudrais programmer sous matlab avec la methode iterative ce qui suit

U (k+1) (j,i)=U(k)(j,i)+(w/4)*(U(k)(j,i+1)+U(k+1)(j,i-1)+U(k)(j,i)+U(k+1)(j-1,i)-4*U(k)(j,i));

(les k et les k+1 sont les iterations, j'ai pas su vous les ecrire comme exposant :s)

avec w appartenant à IR
svp
Configuration: Windows XP
Internet Explorer 6.0

23 réponses

  1. Version avec produit matrice-vecteur par MatLab
    k=0.13;
    c=0.11;
    p=7.8;
    dx=0.25;
    cfl=0.49;
    dt=cfl*(p*c*dx^2)/k;
    a=0;
    b=2;
    T=100;
    nx=(b-a)/dx+1;
    nt=1000;
    x=(a:dx:b);
    U0=100*[(0:dx:b/2-dx).' ; 1 ; fliplr(0:dx:b/2-dx).'];
    U=U0;
    plot(x,U);
    axis([0 2 0 100]);
    mat = -2*r*eye(nx-2,nx-2)...
        +r*diag(ones(nx-3,1),1) ...
        + r*diag(ones(nx-3,1),-1);
    mat = eye(nx-2,nx-2)+mat;
    for n=1:nt
    U(2:nx-1) = mat*U(2:nx-1);
    U(1)=0;
    U(nx)=0;
    pause(0.5);
    plot(x,U);
    axis([0 2 0 100]);
    end
    Version avec produit matrice-vecteur à la main
    k=0.13;
    c=0.11;
    p=7.8;
    dx=0.25;
    cfl=0.49;
    dt=cfl*(p*c*dx^2)/k;
    a=0;
    b=2;
    T=100;
    nx=(b-a)/dx+1;
    nt=1000;
    x=(a:dx:b);
    U0=100*[(0:dx:b/2-dx)  1  fliplr(0:dx:b/2-dx)];
    U=U0;
    plot(x,U);
    axis([0 2 0 100]);
    for n=1:nt
    Up = circshift(U,[0 -1]);Up(1)=0;Up(end)=0;
    Um = circshift(U,[0 1]);Um(1)=0;Um(end)=0;
    U = U + r*(Up+Um-2*U);
    pause(0.5);
    plot(x,U);
    axis([0 2 0 100]);
    end
    1. :s:s désolée pour le triple poste mais je trouve un autre prbleme :s
      comment je fait pour enlever les 0 des bord du U0 que tu m'as passé?
      si je veux juste que je commence par le 0,25...?
      1. désolée j'ai pas vu ce que t'avais poster, j'integre où ce U0 dans mon prgramme (qui est faux :( )?
        1. clear;
          clc;
          %les donnees
          k=0.13;
          c=0.11;
          p=7.8;
          dx=0.25;
          dt=0.5;
          a=0;
          b=2;
          nx=(b-a)/dx;
          nt=(b-a)/dx;
          r=(k*dt)/(p*c*dx^2);
          x=0.25
          %la matrice A
          E=zeros(nx-1);
          for i=2:nx-1
          j=i-1;
          E(i,j)=r;
          end;
          C=E';
          D=(1-2*r)*eye(nx-1);
          A=E+C+D

          %le vecteur U
          V(1)=25;
          U=A*V
          for j=1:nt-1
          U=A*U
          end
          U
          x=1:nx-1;
          y=1:nt-1;
          mesh(x,y,U)

          mais c'est pas ca vu que je n'ai pas programmer le U(x,0) qui est une fonction dependante du x :( j'ai pas su le faire
          ... :(
          1. j'ai oublier un truc :s comment je programme mon x que j'ai donner dans l'ennoncé :s parceque j'en aurai besoin pour le U(x,0)?
            stp
            1. De rien ;-)
              Tu me montreras ton programme quand tu l'auras écrit ?
              Et oublie pas de virer U à chaque extrémité du segment dans ton équation, c'est pas un ddl, il est toujours égal à zéro.
              1. voila j'ai pris des cas partiulier, comme quoi la matrice A est E dans le programme qui suit et le vecteur U est P
                j'ai preciser les limites et j'ai donné au départ le vecteur initial U comme a

                a(1)=25;
                U=zeros(8,1);
                A=zeros(8);
                for i=2:8
                j=i-1;
                A(i,j)=1;
                end;
                C=A';
                D=-4*eye(8);
                E=A+C+D
                P=E*a
                for j=1:8
                P=E*P
                end
                P

                alors est ce que mon programme doit etre dans les environs de ce qui se trouve en dessus
                stp !!
                1. C'est un produit matriciel, il y a rien à coder, MatLab le fait.
                  Cela dit, vu que la matrice A est creuse, c'est vrai que faire un produit matrice-vecteur est un peu idiot, vaut mieux coder à la main.
                  L'équation de la chaleur que tu veux coder, c'est pour une condition de Dirichlet homogène au bord (la température vaut 0 aux extrémités), donc tes degrés de libertés sont seulement les valeurs de la température en tous les points de discrétisations du segment, sauf le premier et le dernier.
                  Si t'as U au temps j, en faisant :
                  Up = circshift(U,[0 -1]);Up(1)=0;Up(end)=0;
                  Um = circshift(U,[0 1]);Um(1)=0;Um(end)=0;
                  U = U + r*(Up+Um-2*U);
                  T'as U au temps j+1.
                  Si tu veux stocker tous les résultats, (temps et espace), tu peux le faire sous forme de matrice. Chaque ligne de la matrice pour un temps donné et chaque colonne pour une position donnée.
                  1. resalut ^^ !
                    j'ai reussi a programmer mon U(i) initial (d'apres les donnees) et j'ai implementé mon A(i,j)
                    donc maintenant il ne reste plus que programmer leur produit, mais je sais toujours pas où intervient les iteration :s et comment les programmer :s
                    j'ai lu le sujet que tu m'as passé mais il m'aide pas trop :s (je comprend mais je sais pas faire :s )
                    alors si ca te derrage pas de m'expliquer de la maniere la plus bete ^^"
                    1. OK, c'est bien ce que je me disais au début, c'était donc un schéma numérique pour l'équation de la chaleur.
                      Mais vu que le schéma ressemblait à rien du tout, je m'étais dit que c'était peut-être pour autre chose.

                      Si tu l'avais dit depuis le début, je te l'aurais dit directement. L'équation de ton tout premier message, si tu l'as bien recopiée, est complètement fausse.

                      Pour ton autre question, rien de compliqué, c'est un schéma explicite, noté vectoriellement tu as
                      U(j+1)=A*U(j)
                      T'as juste à implémenter A.

                      Tu peux jeter un œil ici si tu veux
                      http://www.commentcamarche.net/forum/affich 6135901 resolution de l equation de chaleur
                      1. l'equation de l'iteration est celle de laplace :s et en effet y a des conditions qui manquent :s :s donc l’erreur vient peut etre du prof
                        Là je suis encore bloquée sur un autre programme, si quelqu’un peut m’aider svp
                        Il s’agit de resoudre l’equation de la chaleur suivante
                        du/dt=(k/p*c)*d²u/dx²
                        K=0.13 c=0.11 p=7.8 dx=0.25 et dt est donné par dt<((dx)²*p*c)/(2*k)
                        U(x,o)=100*x si x appartient à [0,1]
                        U(x,o)=100*(2-x) si x appartient à [1,2]
                        U(0,t)=0
                        U(2,t)=0
                        Faudrait resoudre cette equation en utilisant la methode iterative suivante :
                        (les j et les j+1 sont les iterations)
                        U(j+1)(i)=r*(U((j)(i+1)+U(j)(i-1))+(1-2*r)U(j)(i)
                        r=(k/p*c)*(dt/(dx)²)
                        1. D'accord !
                          Bon, pour le reste, ça va être dur de t'aider vu que je sais pas ce que tu veux mettre exactement comme initialisation pour U dans ta méthode itérative, mais pour ton problème d'itération, il faut que tu fasses passer tous les termes de l'itération k+1 à gauche.
                          U(k+1) (j,i)-w/4*U(k+1)(j,i-1)-w/4*+U­(k+1)(j-1,i)=U(k)(j,i)+(w/4)*(U(k)(j,i+1)+U(k)(j,i)-4*U(k)(j,i))
                          Ensuite changer l'ordre des lignes ou des colonnes d'une matrice, ça revient à multiplier par une matrice de permutation.

                          Cela dit, il manque une info dans ton énoncé, pour les première, dernière lignes et pour les première, dernière colonnes, quelle est l'équation pour l'itération ?
                          1. ben justement c'est là où se présume mon problème, je sais pas comment faire ma programmation alors qu'il y a des k+1 à droite :s
                            1. clc
                              clear
                              ep=10^(-10000);%ep=epsilon
                              h=2.5;
                              a=0;
                              b=20;
                              c=0;
                              d=10;
                              nx=(b-a)/h;
                              ny=(d-c)/h;
                              n=(nx+1)*(ny+1);
                              %t=cos(pi/nx)+cos(pi/ny);
                              %w=4/(2+sqrt(4-t^2))
                              w=1
                              M=ones(n,1);

                              %la déclaration et l'initialisation de U
                              sol=zeros(n,1);
                              for i=(nx+3):(n-(nx+1))
                              if (mod(i,(nx+1))~=0 & mod(i,(nx+1))~=1 )
                              sol(i)=i;
                              elseif mod(i,(nx+1))==0
                              sol(i)=100;
                              end
                              end
                              sol;
                              %la transformation du vecteur U en une matrice "sol" de dimension (ny+1,nx+1)

                              for j=1:(ny+1)
                              for i=1:(nx+1)
                              U(j,i)=sol(i+(nx+1)*(j-1));
                              end
                              end
                              U

                              % l est le nombre d'itérations
                              L=0;

                              while(abs(norm(U)-norm(M))>ep)
                              M=U;
                              for i=2:nx-1
                              for j=2:ny-1
                              U(j,i)=U(j,i)+(w/4)*(U(j,i+1)+U(j,i-1)+U(j,i)+U(j-1,i)-4*U(j,i));
                              end
                              end
                              L=L+1;
                              end

                              U
                              L

                              % la trace de la matrice "sol"
                              x=1:nx+1;
                              y=1:ny+1;
                              mesh(x,y,U)

                              voilà le programme que j'ai fait, mais il me donne pas la solution auquelle je m'attendais, je crois qu'à la ligne où je mets le if (mod(i,(nx+1))~=0 & mod(i,(nx+1))~=1 )
                              sol(i)=i;
                              est fausse, et je trouve un prbleme aussi là où j'essaie de programmer mon U(i,j) que je vous ai donné en premier :s
                              alors si vous pouviez m'aider
                              • 1
                              • 2