%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function a=laplaceD(n) // // calcule la matrice du laplacien 1D // discrétisé par différences finies centrées // conditions de Dirichlet a=zeros(n,n) for i=1:n-1 a(i,i)=2;a(i,i+1)=-1;a(i+1,i)=-1; end; a(n,n)=2; a=a*(n+1)^2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function a=laplace(n) // // calcule la matrice du laplacien 1D // discrétisé par différences finies centrées // conditions de Dirichlet a=zeros(n,n) for i=1:n-1 a(i,i)=2;a(i,i+1)=-1;a(i+1,i)=-1; end; a(n,n)=2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function a=laplaceP(n) // // calcule la matrice du laplacien 1D // discrétisé par différences finies centrées // conditions Periodiques a=laplaceD(n);a(1,n)=-1;a(n,1)=-1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [x,sol]=solvelaplace(a,b,ua,ub,c,f,n) // //resout le pobleme -u''+cu=f //avec CL de Dirichlet //discretisation de [a,b] x=linspace(a,b,n+2)'; h=(a-b)/(n+2); // le x utile x(n+2)=[]; x(1)=[]; //calcul du second membre bb=h*h*f(x); bb(1)=bb(1)+ua; bb(n)=bb(n)+ub; //calcul de la matrice A=laplaceD(n); A=A+h*h*diag(c(x)); //calcul de la solution sol=A\bb; x=[a;x;b]; sol=[ua;sol;ub] %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A=laplace2dD(n) // // Calul de la matrice du laplacien en 2D // posé sur le carré unité // conditions aux limites de Dirichlet // n = nombre de points internes en x // = nombre de points internes en x I0=-eye(n,n); B0=2*eye(n,n)-diag(ones(n-1,1),1);B0=B0+B0'; O=zeros(n,n); A=[B0 I0;I0 B0]; X=O; for j=3:n A=[A [X I0]';X I0 B0]; X=[X O]; end; A=A*(n+1)*(n+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A=laplace2dDBis(n) // A=4*eye(n*n,n*n); for i=1:n*n-1 A(i,i+1)=-1; A(i+1,i)=-1; end; for i=1:n-1 A(n*i,n*i+1)=0; A(n*i+1,n*i)=0; end; for i=1:(n-1)*(n-1) A(i,i+n)=-1; A(i+n,i)=-1; end; A=A*(n+1)*(n+1); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function A=laplace2dDSparse(n) // // Calul de la matrice sparse du laplacien en 2D // posé sur le carré unité // conditions aux limites de Dirichlet // n = nombre de points internes en x // = nombre de points internes en x // bloc 1 i=1; ii=[1,1,1];jj=[i,i+1,i+n];uu=[4,-1,-1]; for i=2:n-1 ii=[ii,i,i,i,i];jj=[jj,i-1,i,i+1,i+n]; uu=[uu,-1,4,-1,-1]; end; i=n; ii=[ii,i,i,i];jj=[jj,i-1,i,i+n];uu=[uu,-1,4,-1]; // blocs 2 à n-1 for k=2:n-1 I=(k-1)*n; i=I+1; ii=[ii,i,i,i,i];jj=[jj,i-n,i,i+1,i+n]; uu=[uu, -1,4,-1,-1]; for i=I+2:I+n-1 ii=[ii,i,i,i,i,i];jj=[jj,i-n,i-1,i,i+1,i+n]; uu=[uu,-1,-1,4,-1,-1]; end; i=I+n; ii=[ii,i,i,i,i];jj=[jj,i-n,i-1,i,i+n]; uu=[uu,-1,-1,4,-1]; end; // bloc n i=(n-1)*n+1; ii=[ii,i,i,i];jj=[jj,i-n,i,i+1]; uu=[uu,-1,4,-1]; for i=n*(n-1)+2:n*n-1 ii=[ii,i,i,i,i];jj=[jj,i-n,i-1,i,i+1]; uu=[uu,-1,-1,4,-1]; end; i=n*n; ii=[ii,i,i,i];jj=[jj,i-n,i-1,i];uu=[uu,-1,-1,4]; uu=uu*(n+1)*(n+1); A=sparse([ii;jj]',uu); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function b=Smlaplace2dD(n,f) // // Calul le second membre du laplacien en 2D grillex=(1:n)/(n+1);grilley=(1:n)/(n+1); xx=grillex'*ones(1,n); yy=ones(n,1)*grilley; b=f(xx,yy);b=b(:); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sol=SolveLaplace2d(n,f) // // résout le laplacien en 2D // posé sur le carré unité // conditions aux limites homogènes // f = fonction donnant le second membre A=laplace2dD(n); b=Smlaplace2dD(n,f); sol=A\b; // sol est un vecteur sol=matrix(sol,n,n);// sol est une matrice %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sol=SolveLaplace2dSparse(n,f) // // résout le laplacien en 2D // posé sur le carré unité // conditions aux limites homogènes // f = fonction donnant le second membre A=laplace2dDSparse(n); b=Smlaplace2dD(n,f); sol=A\b; // sol est un vecteur sol=matrix(sol,n,n);// sol est une matrice