/*  TP1 Differences finies, programme en scilab a traduire en C */

#include 
#include 
#define N 50
#define MAX N*N

void produit(float x[MAX],float y[MAX],float alfa)
{float L=4;
int i,j,k,n=N/3;
y[0]=4*x[0]-2*x[1]-2*x[N];
y[N-1]=4*x[N-1]-2*x[N-2]-2*x[2*N-1];
y[N*N-N]=4*x[N*N-N]-2*x[N*N-2*N]-2*x[N*N-N+1];
y[N*N-1]=4*x[N*N-1]-2*x[N*N-2]-2*x[N*N-N-1];

for (k=N*N-N+1;N*N-2;k++)
   {y[k]=4*x[k]-2*x[k-N]-x[k+1]-x[k-1];}

for(k=1;k=N-2;k++)
   {y[k*N]=4*x[k*N]-2*x[k*N+1]-x[k*N-N]-x[k*N+N];}

for(k=n;2*n;k++)
{y[k]=x[k];
 y[k*N]=[4+2*alfa*h]*x[k*N]-2*x[k*N-1]-x[k*N+N]-x[k*N-N];}



for (k=1;k=n-1;k++)
  {y[k]=4*x[k]-2*x[k+N]-x[k-1]-x[k+1];}

for(k=2*n+1;k=N-2)
  {y[k]=4*x[k]-2*x[k+N]-x[k-1]-x[k+1];}

for(k=2;k=n-1;k++)
{y[k*N]=4*x[k*N]-2*x[k*N-1]-x[k*N+N]-x[k*N-N];}

for(k=2*n+1;k=N-2;k++)
     {y[k*N]=4*x[k*N]-2*x[k*N-1]-x[k*N+N]-x[k*N-N];}

for(k=1;k=N-2,k++)
{for(j=k*N+1;j=(k+1)*N-2;j++)
  {y[j]=4*x[j]-x[j-N]-x[j+N]-x[j+1]-x[j-1];}
}

void transproduit(float x[MAX], float y[MAX], float alfa)
{y[0]=4*x[0]-x[1]-x[N];
y[1]=4*x[1]-2*x[0]-x[2]-x[N+1];
y[n-1]=4*x[n-1]-x[n-2]-x[n-1+N];
y[n]=x[n]-x[n-1]-x[n+N];
y[2*n]=x[2*n]-x[2*n+1]-x[2*n+N];
y[2*n+1]=4*x[2*n+1]-x[2*n+2]-x[2*n+N+1];
y[N-2]=4*x[N-2]-2*x[N-1]-x[N-3]-x[2*N-2];
y[N-1]=4*x[N-1]-x[N-2]-x[2*N-1];
y[N]=4*x[N]-2*x[0]-x[N+1]-x[2*N];
y[N+1]=-2*x[N]+4*x[N+1]-x[N+2]-2*x[1]-x[2*N+1];
y[2*N-2]=4*x[2*N-2]-2*x[N-2]-2*x[2*N-1]-x[2*N-3]-x[3*N-2];
y[2*N-1]=4*x[2*N-1]-2*x[N-1]-x[2*N-2]-x[3*N-1];

 int p=N*N-2*N;/*les 2 premieres lignes de l'avant dernier bloc*/
y[p]=4*x[p]-x[p+1]-2*x[p+N]-x[p-N];
y[p+1]=-2*x[p]+4*x[p+1]-x[p+2]-2*x[p+1+N]-x[p+1-N];
  
 int q=N*N-N-2; /*les 2 dernieres lignes de l'avant-dernier bloc*/
y[q]=4*x[q]-2*x[q+1]-x[q-1]-2*x[q+N]-x[q-N];
y[q+1]=-x[q]+4*x[q+1]-x[q+1-N]-2*x[q+1+N];
  
 int r=N*N-N;/* les 2 premieres lignes du dernier bloc*/
y[r]=4*x[r]-x[r+1]-x[r-N];
y[r+1]=-2*x[r]+4*x[r+1]-x[r+2]-x[r+1-N];
  
 int s=N*N-1; /*les 2 dernieres lignes du dernier bloc */
y[s-1]=4*x[s-1]-2*x[s]-x[s-2]-x[s-1-N];
y[s]=4*x[s]-x[s-1]-x[s-N];


  
// les 2 premiers blocs

for(k=2;k=n-2;k++)
  {y[k]=4*x[k]-x[k-1]-x[k+1]-x[k+N];}

 for(k=2*n+2;k=N-3;k++)
   {y[k]=4*x[k]-x[k-1]-x[k+1]-x[k+N];}

 for(k=n+1;k=2*n-1;k++)
       {y[k]=x[k]-x[k+N];}

for(k=N+2;k=N+n-1;k++)
y[k]=4*x[k]-2*x[k-N]-x[k-1]-x[k+1]-x[k+N];
end

for(k=N+2*n+1;k=2*N-3;k++)
y[k]=4*x[k]-2*x[k-N]-x[k-1]-x[k+1]-x[k+N];
end

  for(k=N+n;k=N+2*n;k++)
y[k]=4*x[k]-x[k-1]-x[k+1]-x[k+N];
end

// les 2 derniers blocs

for(k=N*N-2*N+2;k=N*N-N-3;k++)
  {y[k]=4*x[k]-x[k-1]-x[k+1]-x[k-N]-2*x[k+N];}
for(k=N*N-N+2;k=N*N-3;k++)
   {y[k]=4*x[k]-x[k-1]-x[k+1]-x[k-N];}

// les autres blocs

for(k=2;k=N-3;k++)
{y[k*N]=4*x[k*N]-x[k*N+1]-x[k*N-N]-x[k*N+N];
 y[k*N+1]=4*x[k*N+1]-2*x[k*N]-x[k*N+2]-x[k*N+1-N]-x[k*N+1+N];
 y[k*N+N-2]=4*x[k*N+N-2]-2*x[k*N+N-1]-x[k*N-2]-x[k*N+2*N-2]-x[k*N+N-3];
 y[k*N+N-1]=4*x[k*N+N-1]-x[k*N+N-2]-x[k*N-1]-x[k*N+1*N];}

for(j=3;k=N-2;k++)
{y[k*N+j]=4*x[k*N+j]-x[k*N+j+1]-x[k*N+j-1]-x[k*N+j-N]-x[k*N+j+N];}


// sans oublier les points de la fenetre

for(k=n;k=2*n;k++)
{y[k*N]=y[k*N]+2*h*alfa*x[k*N];}
}


void vect(float b[MAX],float lambda, float Uext, float alfa)
{int k,N=MAX,n=N/3;
 float M=MAX,h=1/M;
for(k=n;k=2*n;k++)
  {b[k]=b[k]+lambda;} /* radiateur */

 for(k=n;2*n;k++)
{b[k*N]=b[k*N]+2*h*alfa*Uext;} /* fenetre */
}

float norm(float x[MAX])
{float r=0;
 int k,N=MAX;
 for (k=0;k=N-1;k++)
   {r=r+x[k]*x[k];}
 return r;}

void resoud(float lambda,float Uext,float alfa)
{float b[MAX],bt[MAX],x[MAX];
 int i,j,k,N=MAX;
b=vect[f,lambda,Uext,alfa,h];
transproduit(b,bt,alfa);

/* resolution de Au=b par la methode du gradient conjugue
 On multiplie le systeme a gauche par la transposee de A
 Ce qui nous donne un systeme symetrique Nu=c avec N=tA*A et c=tA*b
 On part du vecteur [1 1 1 1 1 ...... ] */


 for(k=1:

T=ones[N*N,1];
produit(T,alfa,h);
transproduit(T,alfa,h);
p=c-z1;
r=p;
while norm(r)>0.01,
  produit(p,alfa,h);
  transproduit(z,alfa,h);
  a=r'*r;
  d=z1'*p;
  pas=a/d; 
  T=T+pas*p;
  r1=r-pas*z1;
  e=r1'*r1; 
  g=r'*r;
  beta=e/g;
  p=r1+beta*p;
  r=r1;
  norm[r),
  
end

/*-----------------------------------------------------------------*/
/* TP2 Differences finies - version sans interrogation.           */                             
/*  a rediriger vers scilab                                      */
/*---------------------------------------------------------------*/

#include 
#include 
#define MAX 70

void nextstep(float A[MAX][MAX], float B[MAX][MAX], float
f[MAX][MAX], float alfa,float beta, float Uext)
{int i,j,k,N=MAX;
 float M=MAX;
 float h=4/M,t=h*h/4;
 float c=2*h*alfa;
 int m=N/3;

B[0][0]=A[0][0]+t*f[0][0]+beta*(2*A[0][1]+2*A[1][0]-4*A[0][0]);
B[0][N-1]=A[0][N-1]+t*f[0][N-1]+beta*(2*A[1][N-1]+2*A[0][N-2]-4*A[0][N-1]);
B[N-1][0]=A[N-1][0]+t*f[N-1][0]+beta*(2*A[N-1][1]+2*A[N-2][1]-4*A[N-1][0]);
B[N-1][N-1]=A[N-1][N-1]+t*f[N-1][N-1]+beta*(2*A[N-2][N-1]+2*A[N-1][N-2]-4*A[N-1][N-1]);


for(k=1;k<=N-2;k++)
{B[0][k]=A[0][k]+t*f[0][k]-beta*(4*A[0][k]-2*A[1][k]-A[0][k+1]-A[0][k-1]);

B[N-1][k]=A[N-1][k]+t*f[N-1][k]-beta*(4*A[N-1][k]-2*A[N-2][k]-A[N-1][k+1]-A[N-1][k-1]);}

for(k=1;k<=m-1;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}


for(k=2*m+1;k<=N-2;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}

for (k=m;k<=2*m;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+c*Uext+A[k+1][0]+A[k-1][0]-(4+c)*A[k][0]);
 B[k][N-1]=A[k][N-1];}

for(i=1;i<=N-2;i++)
{for(j=1;j<=N-2;j++)

{B[i][j]=A[i][j]+t*f[i][j]+beta*(A[i-1][j]+A[i+1][j]+A[i][j+1]+A[i][j-1]-4*A[i][j]);}
}}


void copier(float A[MAX][MAX],float B[MAX][MAX])
{int i,j,N=MAX;

for(i=0;i<=N-1;i++)
 {for(j=0;j<=N-1;j++)
 {B[i][j]=A[i][j];}
 }
}


void afficher(float A[MAX][MAX])
{int i,j,N=MAX;
 printf("T=[");
for(i=0;i<=N-2;i++)
{for(j=0;j<=N-1;j++)
{printf("%3.1f ",A[i][j]);}
 printf(";\n");}
for(j=0;j<=N-1;j++)
{printf("%3.1f ",A[i][j]);}
printf ("];\n");
printf("x=linspace(0,4,70);\ny=x;\nplot3d(x,y,T)");}

main()
{float A[MAX][MAX];
float B[MAX][MAX];
float f[MAX][MAX];
float C[MAX][MAX];
float alfa=0.1;

int i,j,k,P,N=MAX,m=N/3;
 float M=MAX;
float h=4/M,t=h*h/4;
P=50/t;

for(i=0;i<=N-1;i++)
{for(j=0;j<=N-1;j++)
 {A[i][j]=15;

 f[i][j]=0;}}

 for(k=m;k<2*m;k++)
   {A[k][N-1]=28;}


for(i=1;i<=P;i++)
{nextstep(A,B,f,alfa,0.25,-10);
 copier(B,A); }
 afficher(A);
 printf("\n");}



/* TP2 differences finies en C. Autre version avec interrogation */


#include 
#include 
#define MAX 28

void nextstep(float A[MAX][MAX], float B[MAX][MAX], float
f[MAX][MAX], float alfa,float beta, float Uext)
{int i,j,k,N=MAX;
 float M=MAX;
 float h=4/M,t=h*h/4;
 float c=2*h*alfa;
 int m=N/3;

B[0][0]=A[0][0]+t*f[0][0]+beta*(2*A[0][1]+2*A[1][0]-4*A[0][0]);
B[0][N-1]=A[0][N-1]+t*f[0][N-1]+beta*(2*A[1][N-1]+2*A[0][N-2]-4*A[0][N-1]);
B[N-1][0]=A[N-1][0]+t*f[N-1][0]+beta*(2*A[N-1][1]+2*A[N-2][1]-4*A[N-1][0]);
B[N-1][N-1]=A[N-1][N-1]+t*f[N-1][N-1]+beta*(2*A[N-2][N-1]+2*A[N-1][N-2]-4*A[N-1][N-1]);


for(k=1;k<=N-2;k++)
{B[0][k]=A[0][k]+t*f[0][k]-beta*(4*A[0][k]-2*A[1][k]-A[0][k+1]-A[0][k-1]);

B[N-1][k]=A[N-1][k]+t*f[N-1][k]-beta*(4*A[N-1][k]-2*A[N-2][k]-A[N-1][k+1]-A[N-1][k-1]);}

for(k=1;k<=m-1;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}


for(k=2*m+1;k<=N-2;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+A[k-1][0]+A[k+1][0]-4*A[k][0]);

B[k][N-1]=A[k][N-1]+t*f[k][N-1]+beta*(2*A[k][N-2]+A[k+1][N-1]+A[k-1][N-1]-4*A[k][N-1]);}

for (k=m;k<=2*m;k++)
{B[k][0]=A[k][0]+t*f[k][0]+beta*(2*A[k][1]+c*Uext+A[k+1][0]+A[k-1][0]-(4+c)*A[k][0]);
 B[k][N-1]=A[k][N-1];}

for(i=1;i<=N-2;i++)
{for(j=1;j<=N-2;j++)

{B[i][j]=A[i][j]+t*f[i][j]+beta*(A[i-1][j]+A[i+1][j]+A[i][j+1]+A[i][j-1]-4*A[i][j]);}
}}


void copier(float A[MAX][MAX],float B[MAX][MAX])
{int i,j,N=MAX;

for(i=0;i<=N-1;i++)
 {for(j=0;j<=N-1;j++)
 {B[i][j]=A[i][j];}
 }
}


void afficher(float A[MAX][MAX])
{int i,j,N=MAX;
for(i=0;i<=N-1;i++)
{for(j=0;j<=N-1;j++)
{printf("%3.1f ",A[i][j]);}
 printf("\n");}
}


main()
{float A[MAX][MAX];
float B[MAX][MAX];
float f[MAX][MAX];
float C[MAX][MAX];
float alfa=0.1;

int i,j,k,P,N=MAX,m=N/3;
 float M=MAX;
float h=4/M,t=h*h/4;
 float heure,lambda,Uext;

printf("Vous voulez la temperature a quelle heure: ");
scanf("%f",&heure);
printf("Rentrez la temperature du radiateur: ");
scanf("%f",&lambda);
printf("Rentrez la temperature exterieure: ");
scanf("%f",&Uext);
 P=heure/t;

for(i=0;i<=N-1;i++)
{for(j=0;j<=N-1;j++)
 {A[i][j]=15;

 f[i][j]=0;}}

 for(k=m;k<2*m;k++)
   {A[k][N-1]=lambda;}


for(i=1;i<=P;i++)
{nextstep(A,B,f,alfa,0.25,Uext);
 copier(B,A); }
 afficher(A);
 printf("\n");}

/*--------------------------------------------------------------*/
/* 			TP1 d'AVP  				*/
/*--------------------------------------------------------------*/

function y=solution(N,g0,g1,teta,lambda,f)
N=100; l=40;
h=1/N;
x=linspace(0,1,N);
f=x.*(1-x); f=f';  f1=f(1:l-1); f2=f(l:N);
x1=ones(l-2,1);
x2=ones(N-l,1); x3=x2;
x3(1)=2;
A1=(2*eye(l-1,l-1)-diag(x1,1)-diag(x1,-1))/h;
A2=(2*eye(N-l+1,N-l+1)-diag(x3,1)-diag(x2,-1))/h;

for k=1:500,

b1(1)=g0/h; b1(l-1)=lambda/h;
b1=b1+h*f1;
u1=A1\b1;
b2(1)=2*(u1(l-1)-lambda)/h; b2(N-l+1)=g1/h;
b2=b2+h*f2;
u2=A2\b2;
u=[u1;u2];
lambda=teta*u(l)+(1-teta)*lambda;
u(l)=lambda
plot(x,u)

end

/* TP1 differences finies version Scilab */

function [y]=produit(x,alfa,h) // fonction retournant le produit y=Ax
l=4; L=4; // On a une piece de 3m de largeur et de 4 m de longueur
N=ceil(l/h);    n=ceil(N/3);
M=ceil(L/h);    m=ceil(M/3);
y(1)=4*x(1)-2*x(2)-2*x(N+1);
y(N)=4*x(N)-2*x(N-1)-2*x(2*N);
y(N*M-N+1)=4*x(N*M-N+1)-2*x(N*M-2*N+1)-2*x(N*M-N+2);
y(N*M)=4*x(N*M)-2*x(N*M-1)-2*x(N*M-N); // points du coin
  
for k=N*M-N+2:N*M-1,
y(k)=4*x(k)-2*x(k-N)-x(k+1)-x(k-1);
end  // mur droit

for k=1:M-2,
y(k*N+1)=4*x(k*N+1)-2*x(k*N+2)-x(k*N+1-N)-x(k*N+1+N);
end // mur bas

for k=n:2*n,
y(k)=x(k);  // radiateur
end

for k=m:2*m,
y(k*N)=(4+2*alfa*h)*x(k*N)-2*x(k*N-1)-x(k*N+N)-x(k*N-N); // fenetre
end

for k=[2:n-1 2*n+1:N-1],
y(k)=4*x(k)-2*x(k+N)-x(k-1)-x(k+1);
end
for k=[2:m-1 2*m+1:N-1],
y(k*N)=4*x(k*N)-2*x(k*N-1)-x(k*N+N)-x(k*N-N);
end     // points du mur situes de part et d'autre de la fenetre et
       // du radiateur

// restent les points de l'interieur

for k=1:M-2,
for j=k*N+2:(k+1)*N-1,
y(j)=4*x(j)-x(j-N)-x(j+N)-x(j+1)-x(j-1);
 end
end


y(1)=4*x(1)-2*x(2)-2*x(N+1);
y(N)=4*x(N)-2*x(N-1)-2*x(2*N);
y(N*M-N+1)=4*x(N*M-N+1)-2*x(N*M-2*N+1)-2*x(N*M-N+2);
y(N*M)=4*x(N*M)-2*x(N*M-1)-2*x(N*M-N);

// De meme avec la transposee, dont la forme est bien plus compliquee.

function [y]=transproduit(x,alfa,h)
l=4; L=4;
N=ceil(l/h);    n=ceil(N/3);
M=ceil(L/h);    m=ceil(M/3);
 // illisible, mais malheureusement necessaire
y(1)=4*x(1)-x(2)-x(N+1);
y(2)=4*x(2)-2*x(1)-x(3)-x(N+2);
y(n-1)=4*x(n-1)-x(n-2)-x(n-1+N);
y(n)=x(n)-x(n-1)-x(n+N);
y(2*n)=x(2*n)-x(2*n+1)-x(2*n+N);
y(2*n+1)=4*x(2*n+1)-x(2*n+2)-x(2*n+N+1);
y(N-1)=4*x(N-1)-2*x(N)-x(N-2)-x(2*N-1);
y(N)=4*x(N)-x(N-1)-x(2*N);   // premier bloc
y(N+1)=4*x(N+1)-2*x(1)-x(N+2)-x(2*N+1);
y(N+2)=-2*x(N+1)+4*x(N+2)-x(N+3)-2*x(2)-x(2*N+2);
y(2*N-1)=4*x(2*N-1)-2*x(N-1)-2*x(2*N)-x(2*N-2)-x(3*N-1);
y(2*N)=4*x(2*N)-2*x(N)-x(2*N-1)-x(3*N);  // second bloc

  p=N*M-2*N+1; // les 2 premieres lignes de l'avant dernier bloc
y(p)=4*x(p)-x(p+1)-2*x(p+N)-x(p-N);
y(p+1)=-2*x(p)+4*x(p+1)-x(p+2)-2*x(p+1+N)-x(p+1-N);
  
  q=N*M-N-1; // les 2 dernieres lignes de l'avant-dernier bloc
y(q)=4*x(q)-2*x(q+1)-x(q-1)-2*x(q+N)-x(q-N);
y(q+1)=-x(q)+4*x(q+1)-x(q+1-N)-2*x(q+1+N);
  
  r=N*M-N+1;// les 2 premieres lignes du dernier bloc
y(r)=4*x(r)-x(r+1)-x(r-N);
y(r+1)=-2*x(r)+4*x(r+1)-x(r+2)-x(r+1-N);
  
  s=N*M; // les 2 dernieres lignes du dernier bloc
y(s-1)=4*x(s-1)-2*x(s)-x(s-2)-x(s-1-N);
y(s)=4*x(s)-x(s-1)-x(s-N);


  
// les 2 premiers blocs

for k=[3:n-2 2*n+2:N-2],
y(k)=4*x(k)-x(k-1)-x(k+1)-x(k+N);
end

for k=n+1:2*n-1,
y(k)=x(k)-x(k+N);
end

for k=[N+3:N+n-1 N+2*n+1:2*N-2],
y(k)=4*x(k)-2*x(k-N)-x(k-1)-x(k+1)-x(k+N);
end

for k=N+n:N+2*n,
y(k)=4*x(k)-x(k-1)-x(k+1)-x(k+N);
end

// les 2 derniers blocs

for k=N*M-2*N+3:N*M-N-2,
y(k)=4*x(k)-x(k-1)-x(k+1)-x(k-N)-2*x(k+N);
end

for k=N*M-N+3:N*M-2,
y(k)=4*x(k)-x(k-1)-x(k+1)-x(k-N);
end

// les autres blocs

for k=2:M-3,
y(k*N+1)=4*x(k*N+1)-x(k*N+2)-x(k*N+1-N)-x(k*N+1+N);
y(k*N+2)=4*x(k*N+2)-2*x(k*N+1)-x(k*N+3)-x(k*N+2-N)-x(k*N+2+N);
y(k*N+N-1)=4*x(k*N+N-1)-2*x(k*N+N)-x(k*N-1)-x(k*N+2*N-1)-x(k*N+N-2);
y(k*N+N)=4*x(k*N+N)-x(k*N+N-1)-x(k*N)-x(k*N+2*N);
for j=3:N-2,
y(k*N+j)=4*x(k*N+j)-x(k*N+j+1)-x(k*N+j-1)-x(k*N+j-N)-x(k*N+j+N);
  end
end

// sans oublier les points de la fenetre

for k=m:2*m,
y(k*N)=y(k*N)+2*h*alfa*x(k*N);
end


function [b]=vect(f,lambda,Uext,alfa,h)
l=4; L=4; 
N=ceil(l/h);    n=ceil(N/3);
M=ceil(L/h);    m=ceil(M/3);
for k=1:N*M,      // on pourrait essayer avec d'autres fonctions
  b(k)=h*h*f(k);
end

for k=n:2*n,
b(k)=b(k)+lambda; // radiateur
end

for k=m:2*m
b(k*N)=b(k*N)+2*h*alfa*Uext; // fenetre
end


function [T]=resoud(f,lambda,Uext,alfa)
l=4; L=4;
h=0.01;
N=ceil(l/h);    n=ceil(N/3);
M=ceil(L/h);    m=ceil(M/3);

// On pourra changer le pas de discretisation a notre guise

b=vect(f,lambda,Uext,alfa,h);
c=transproduit(b,alfa,h);

// resolution de Au=b par la methode du gradient conjugue
// On multiplie le systeme a gauche par la transposee de A
// Ce qui nous donne un systeme symetrique Mu=c avec M=tA*A et c=tA*b
// On part du vecteur [1 1 1 1 1 ...... ]

T=ones(N*M,1);
z=produit(T,alfa,h);
z1=transproduit(T,alfa,h);
p=c-z1; // Po=Ro=c-MXo
r=p;
while norm(r)>0.01,
  z=produit(p,alfa,h);
  z1=transproduit(z,alfa,h);
  a=r'*r;
  d=z1'*p;
  pas=a/d; 
  T=T+pas*p;
  r1=r-pas*z1;
  e=r1'*r1; 
  g=r'*r;
  beta=e/g;
  p=r1+beta*p;
  r=r1;
  norm(r),
  
end

function dessine(lambda,Uext,alfa,N)
f=zeros(N*N,1);
T=resoud(f,lambda,Uext,alfa,N);
x=linspace(0,3,N);
y=linspace(0,4,M);
for i=0:M-1,
for j=1:N,
  R(j,i+1)=T(N*i+j);
end
end
 
// On trace la nappe qui represente le profil de la temperature
// selon la position dans la piece. On retrouve bien les hautes
// temperatures pres du radiateur
plot3d(x,y,R)

function [y]=lambdaopt(f,Uext,alfa,Zd,N)
// en principe, on prend une temperature constante sur toute la
// piece. Zd peut etre soit une constante, soit une matrice contenant
// le profil de la temperature desire.
U0=resoud(f,0,Uext,alfa,N);
U1=resoud(zeros(N*N,1);1,0,alfa,N);
P0=resoud(U0-Zd,0,0,0,N);
P1=resoud(U1,0,0,0,N);

// l'integrale sur les bords de du/dn
n=ceil(N/3);
for k=1:n,
  p=k+n-1;  // indice variant de n a 2n. On procede ainsi pour eviter
  // de stocker des zeros inutilement dans les vecteurs D1, D2, D3, D4
  D1(k)=P0(p)-P0(p+N);
  D2(k)=U1(p*N)*P0(p*N)-U1(p*N-1)*P0(p*N-1);
  D3(k)=P0(p)-P0(p+N);
  D4(k)=U1(p*N)*P1(p*N)-U1(p*N-1)*P1(p*N-1);
end

I1=sum(D1); I2=sum(D2);
I3=sum(D3); I4=sum(D3);
y=-(I1+alfa*I2)/(I3+alfa*I4);

function y=cout(T,Uext,alfa,N)
f=zeros(N*N,1);
x=lambdaopt(f,Uext,alfa,N);
U=resoud(f,x,Uext,alfa,N);

n=ceil(N/3);
for k=1:n,
  p=k+n-1;
  D1(k)=U(p)-U(p+N);
end

y=sum(D1);

/* TP2 differences finies, version Scilab */

function B=pasuivant(A,beta,N)
Uext=-10; // On travaille sur une piece carree de 4m de cote
alfa=0.1;  // Le nombre t*mu/(h*h) doit rester petit, si on ne
m=ceil(N/3); // veut pas que le schema explicite explose
B(1,1)=A(1,1)-beta*(4*A(1,1)-2*A(1,2)-2*A(2,1));
B(1,N)=A(1,N)-beta*(4*A(1,N)-2*A(2,N)-2*A(1,N-1));
B(N,1)=A(N,1)-beta*(4*A(N,1)-2*A(N,2)-2*A(N-1,1));
B(N,N)=A(N,N)-beta*(4*A(N,N)-2*A(N,N-1)-2*A(N-1,N));

for k=2:N-1,
B(1,k)=A(1,k)-beta*(4*A(1,k)-2*A(2,k)-A(1,k+1)-A(1,k-1));
B(N,k)=A(N,k)-beta*(4*A(N,k)-2*A(N-1,k)-A(N,k+1)-A(N,k-1));
end

for k=[2:m-1 2*m+1:N-1],
B(k,1)=A(k,1)-beta*(4*A(k,1)-2*A(k,2)-A(k-1,1)-A(k+1,1));
B(k,N)=A(k,N)-beta*(4*A(k,N)-2*A(k,2)-A(k-1,N)-A(k+1,N));
end

for k=m:2*m,
B(k,1)=A(k,1)-beta*((4+2*alfa*h)*A(k,1)-2*A(k,2)-2*h*alfa*Uext-A(k+1,1)-A(k-1,1));
B(k,N)=A(k,N);
end


for i=2:N-1,
for j=2:N-1,
B(i,j)=A(i,j)-beta*(4*A(i,j)-A(i-1,j)-A(i+1,j)-A(i,j+1)-A(i,j-1));
 end
end

function T=explicit(beta,N,P)
h=4/N;
m=ceil(N/3);
t=h*h/2;
T=15*ones(N,N);
T(m:2*m,N)=25;
x=linspace(0,4,N);
y=x;

for k=1:P,
T=pasuivant(T,beta,N);
end




/* TP de recherche operationnelle */

function [y,g,ind]=costf(q,ind)
p=ones(31,1);
eps=0.001;
a=p'*q;
b=(A*q-b)'(A*q-b);
y=a+b/eps;
g=p+2*A'*(A*q-b)/eps;

/* TP2 differences finies, autres version Scilab */
function [s]=gradient(A,b,x)
M=A'*A;
p=A'*b-M*x;
r=p;

while norm(r)>0.01,
  z=M*p;
  a=r'*r;
  d=z'*p;
  alfa=a/d;  
  x=x+alfa*p;  
  r1=r-alfa*z;
  e=r1'*r1;
  f=r'*r;   
  beta=e/f;
  p=r1+beta*p;
  r=r1;
  norm(r),
  
end

function B=nextstep(A,f,N)
h=1/N;
m=ceil(N/3);
B(1,1)=4*A(1,1)-2*A(1,2)-2*A(2,1);
B(1,N)=4*A(1,N)-2*A(2,N)-2*A(1,N-1);
B(N,1)=4*A(N,1)-2*A(N,2)-2*A(N-1,1);
B(N,N)=4*A(N,N)-2*A(N,N-1)-2*A(N-1,N);

for k=2:N-1,
B(1,k)=4*A(1,k)-2*A(2,k)-A(1,k+1)-A(1,k-1);
B(N,k)=4*A(N,k)-2*A(N-1,k)-A(N,k+1)-A(N,k-1);
end

for k=[2:m-1 2*m+1:N-1],
B(k,1)=4*A(k,1)-2*A(k,2)-A(k-1,1)-A(k+1,1);
B(k,N)=4*A(k,N)-2*A(k,2)-A(k-1,N)-A(k+1,N);
end

for k=m:2*m,
B(k,1)=4*A(k,1)-2*A(k,2)-A(k+1,1)-A(k-1,1);
B(k,N)=A(k,N);
end



for i=2:N-1,
for j=2:N-1,
B(i,j)=4*A(i,j)-A(i-1,j)-A(i+1,j)-A(i,j+1)-A(i,j-1);
 end
end

/* TP1 differences finies, trafiqué (résolution directe) */

function [A]=matrice(alfa,N);
h=1/N;
NN=N*N;
m=ceil(N/3);

A(1,1)=4;
A(1,2)=-2;
A(1,N+1)=-2;
A(N,N)=4;
A(N,N-1)=-2;
A(N,2*N)=-2;
A(NN-N+1,NN-N+1)=4;
A(NN-N+1,NN-N+2)=-2;
A(NN-N+1,NN-2*N+1)=-2
A(NN,NN-1)=-2;
A(NN,NN-N)=-2;
A(NN,NN)=4;

for k=[2:m-1 2*m+1:N-1],
  A(k,k)=4;
  A(k,k-1)=-1;
  A(k,k+1)=-1;
  A(k,k+N)=-2;
end

for k=1:N-2;
  for j=2:N-1;
  A(k*N+j,k*N+j)=4;
  A(k*N+j,k*N+j-1)=-1;
  A(k*N+j,k*N+j+1)=-1
  A(k*N+j,k*N+j-N)=-1;
  A(k*N+j,k*N+j+N)=-1;
end
A(k*N+1,k*N+1)=4;
A(k*N+1,k*N+2)=-2;
A(k*N+1,k*N+1+N)=-1;
A(k*N+1,k*N+1-N)=-1;
A(k*N+N,k*N+N)=4;
A(k*N+N,k*N+N-1)=-2;
A(k*N+N,k*N)=-1;
A(k*N+N,k*N+2*N)=-1;
end

for k=NN-N+2:NN-1;
  A(k,k)=4;
  A(k,k+1)=-1;
  A(k,k-1)=-1;
  A(k,k-N)=-2;
end



for k=m:2*m,
  A(k,k)=1;
  A(k*N,k*N)=4+2*h*alfa;
end

function b=vect(f,lambda,Uext,alfa,N)
h=1/N;
n=ceil(N/3);
for k=1:N*N,      // on pourrait essayer avec d'autres fonctions
  b(k)=h*h*f(k);
end

for k=n:2*n,
b(k)=b(k)+lambda; // radiateur
b(k*N)=b(k*N)+2*h*alfa*Uext; // fenetre
end


function [T]=resoud(f,lambda,Uext,alfa,N)
A=matrice(alfa,N);
b=vect(f,lambda,Uext,alfa,N);
T=A\b;


function dessine(lambda,Uext,alfa,N)
f=zeros(N*N,1);
T=resoud(f,lambda,Uext,alfa,N);
x=linspace(0,3,N);   x=x';
y=linspace(0,3,N);   y=y';
for i=0:N-1,
for j=1:N,
  Z(j,i+1)=T(N*i+j);
end
end
plot3d(x,y,Z)

function [y]=lambdaopt(f,Zd,Uext,alfa,N)
// en principe, on prend une temperature constante sur toute la
// piece. Zd peut etre soit une constante, soit une matrice contenant
// le profil de la temperature desire.
U0=resoud(f,0,Uext,alfa,N);
U1=resoud(zeros(N*N,1),1,0,alfa,N);
P0=resoud(U0-Zd,0,0,0,N);
P1=resoud(U1,0,0,0,N);

// l'integrale sur les bords de du/dn
n=ceil(N/3);
for k=1:n+1,
  p=k+n-1;  // indice variant de n a 2n. On procede ainsi pour eviter
  // de stocker des zeros inutilement dans les vecteurs D1, D2, D3, D4
  D1(k)=P0(p)-P0(p+N);
  D2(k)=U1(p*N)*P0(p*N)-U1(p*N-1)*P0(p*N-1);
  D3(k)=P1(p)-P1(p+N);
  D4(k)=U1(p*N)*P1(p*N)-U1(p*N-1)*P1(p*N-1);
end

I1=sum(D1); I2=sum(D2);
I3=sum(D3); I4=sum(D3);
y=-(I1+alfa*I2)/(I3+alfa*I4);

function y=cout(T,Uext,alfa,N)
f=zeros(N*N,1);
x=lambdaopt(f,T,Uext,alfa,N);
U=resoud(f,x,Uext,alfa,N);

n=ceil(N/3);
for k=1:n+1,
  p=k+n-1;
  D1(k)=U(p)-U(p+N);
 end

y=sum(D1);

// Farid IKKEN
// DESS IMOI
// Projet d'elements finis pour M. Conrad

function A=matriceA(N,m)
h=1/N;
h2=h*h;
// voici le bloc qu'on va repeter N-2 fois dans la matrice
// on procede par concatenation successive
x=[9/70,13*h/420,26/35,0,9/70,-13*h/420;
  -13*h/420,-h2/140,0,2*h2/105,13*h/420,-h2/140];
y=[x(1:2,3:6) zeros(2,2);x];
// premier bloc, forme des 4 dernieres colonnes de X, et de 4 zeros
// et de X.
for i=1:N-3,
  y=[y zeros(2*i+2,2); zeros(2,2*i) x];
end
y(2*N-1,2*N-3)=9/70;
y(2*N-1,2*N-2)=13*h/420;
y(2*N-1,2*N-1)=13/35;
y(2*N-1,2*N)=-11*h/210;

y(2*N,2*N-3)=-13*h/420;
y(2*N,2*N-2)=-h2/140;
y(2*N,2*N-1)=-11*h/210;
y(2*N,2*N)=h2/105;
A=h*y;
A(2*N-1,2*N-1)=A(2*N-1,2*N-1)+m;

function C=matriceC(N)
h=1/N;
h2=h*h;
h3=h*h*h;
// on procede comme en A
x=[-12,-6*h,24,0,-12,6*h;6*h,2*h2,0,8*h2,-6*h,2*h2];
y=[x(1:2,3:6) zeros(2,2);x];
for i=1:N-3,
  y=[y zeros(2*i+2,2); zeros(2,2*i) x];
end
y(2*N-1,2*N-3)=-12;
y(2*N-1,2*N-2)=-6*h;
y(2*N-1,2*N-1)=12;
y(2*N-1,2*N)=-6*h;

y(2*N,2*N-3)=6*h;
y(2*N,2*N-2)=2*h2
y(2*N,2*N-1)=-6*h;
y(2*N,2*N)=4*h2;
C=y/h3;

function B=matriceB(N,alfa)
h=1/N;
B(2*N,2*N)=alfa;

function [AA,BB]= rondes(N,alfa,m)
Z=ones(2*N,2*N);
I=eye(2*N,2*N);
A=matriceA(N,m);
B=matriceB(N,alfa);
C=matriceC(N);
AA=[-B -C;I Z];
BB=[A Z; Z I];

function y=valprp(N,alfa,m)
[a b]=rondes(N,alfa,m);
x=spec(inv(b)*a);
y=[]; // vecteur vide qu'on va remplir avec des elements de x
for i=1:4*N,    
  im=imag(x(i));
  re=real(x(i));
  // Toutes les valeurs propres situees a gauche
  // de l'axe imaginaire. De toute facon, elles sont toutes
  // a gauches. Toutes les valeurs propres etant
  // conjuguees, on ne prend que celle qui sont de partie
  // imaginaire positive, ce qui nous fait 4*N/2=2*N v.p.
if (im>0 & re <0) then
y=[y;x(i)];
end
end
j=size(y); j=j(1);
//n=floor(1.1*N);
//p=floor(1.2*N);
y=y(N:j);


function [x1,x2]=spectre(N,alfa,m)
x=valprp(N,alfa,m);
j=size(x); j=j(1);
s=-1/alfa;
x1=x(2:2:j);
x2=x(3:2:j);

k=size(x1); k=k(1);
z1=x1(k:-1:1);

square(-5,0,0,1e5)
plot([real(x1)' real(z1)' real(x2)'],[imag(x1)' imag(z1)' imag(x2)'])


function courbe(m)
z=1:9;
alfa=[1e-5*z,1e-4*z,1e-3*z,1e-2*z]
for i=1:36,
  x=valprp(50,alfa(i),m);
  y(i)=max(real(x));
end
plot(y,alfa)

// Farid IKKEN
// DESS IMOI
// TP1 de Volume finis


// 1/ Cas de la discretisation uniforme
// ====================================

function [y] = lambda(T)
y=0.25*(T^0.81) + 0.0004*T + 0.43;

function [a]=coeff(T,N)
for i=1:N+1,
a(i)=2*(lambda(T(i))*lambda(T(i+1)))/(lambda(T(i))+lambda(T(i+1)));
end



function uniforme(N)
T=linspace(400,0,N+2);
s=T;
f=linspace(-20,0,N);
x=linspace(0,0.1,N+2);
h=0.1/N;
b=h*h*f';

// On fait passer le -Ao*To de la premiere equation du cote
// du vecteur image, sachant que To=400 et Ao=2*lambda(To)*lambda(T1)/.....

// ====================================================================
// Cette region doit etre copiee, puis collee dans scilab. 
// Pour voir chaque iteration, cliquer sur le bouton du milieu
// de la souris. La fonction qui devait remplir la matrice
// ne marchait pas, et remplissait une matrice nulle,
// avec le message d'erreur 'singular matrix'

r=s;
a=coeff(s,N);
b(1)=h*h*f(1)+400*a(1);

for i=1:N,
A(i,i)=a(i)+a(i+1);
end

for i=1:N-1,
A(i,i+1)=-a(i+1);
end

for i=2:N,
A(i,i-1)=-a(i);
end

T1=A\b;
s=[400 T1' 0];
plot (x,s);
sum((r-s)^2)

// ====================================================================

// On rappelle que le schema aux volumes finis
// pour la discretisation uniforme s'ecrit
// -A_i(T_i+1-T_i)+A_i-1(T_i-T_i-1)=h2*f_i
// ou f est la fonction S(t)=200(t-20)

// La matrice de ce systeme lineaire s'ecrit:
// 
//
//     | A1+A2      -A2       |
//     |-A2   A2+A3   -A3     | 
//  A= |   \    \        \    |
//     |    \    \        \   |
//     |  0  \    \        \  |
//     |      \    \       -An|
//     |       -An   An + An+1| 
//
// On remplit alors cette matrice



// A chaque iteration k, on trace la courbe de X ----> Tk


// 2/ Cas de la discretisation non uniforme
//==========================================
//
// On genere des pas aleatoires avec la fonction rand de scilab


function [a]=coeff2(T,N)
for i=1:N+1,
a(i)=2*(lambda(T(i))*lambda(T(i+1)))/(lambda(T(i))*h(2*i-1)+lambda(T(i+1))*h(2*i));
end

function not_uniforme(N)
T=linspace(400,0,N+2);
f=linspace(-20,0,N);
x=linspace(0,0.1,N+2);

h=rand(1,2*N+2,'uniform');
// On doit avoir 2N pas et non N a cause des h(i+1/2)+ et h(i+1/2)-

h=0.1*h/sum(h);
// les h(i) doivent etre les pas d'une subdivision de [0,0.1]
// donc la somme doit faire 0.1= longueur de l'intervalle


// ===============================================================
// La encore cette region doit etre copiee

for i=1:N,
  b(i)=(h(2*i)+h(2*i-1))*f(i);
end

r=T;
a=coeff2(T,N);
b(1)= (h(1)+h(2))*f(1)+400*a(1);

for i=1:N,
A(i,i)=a(i)+a(i+1);
end

for i=1:N-1,
A(i,i+1)=-a(i+1);
end

for i=2:N,
A(i,i-1)=-a(i);
end

T1=A\b;
T=[400 T1' 0];
plot (x,T);
sum((r-T)^2)

//==============================================================


function B=nextstep(A,f,Uext,alfa,N)
h=4/N;     // On travaille sur une piece carree de 4m de cote
t=h*h;   // Le nombre t*mu/(h*h) doit rester petit, si on ne
m=ceil(N/3); // veut pas que le schema explicite explose
mu=0.1;
beta=t*mu/(h*h);
B(1,1)=A(1,1)+t*f(1,1)-beta*(4*A(1,1)-2*A(1,2)-2*A(2,1));
B(1,N)=A(1,N)+t*f(1,N)-beta*(4*A(1,N)-2*A(2,N)-2*A(1,N-1));
B(N,1)=A(N,1)+t*f(N,1)-beta*(4*A(N,1)-2*A(N,2)-2*A(N-1,1));
B(N,N)=A(N,N)+t*f(N,N)-beta*(4*A(N,N)-2*A(N,N-1)-2*A(N-1,N));

for k=2:N-1,
B(1,k)=A(1,k)+t*f(1,k)-beta*(4*A(1,k)-2*A(2,k)-A(1,k+1)-A(1,k-1));
B(N,k)=A(N,k)+t*f(N,k)-beta*(4*A(N,k)-2*A(N-1,k)-A(N,k+1)-A(N,k-1));
end

for k=[2:m-1 2*m+1:N-1],
B(k,1)=A(k,1)+t*f(k,1)-beta*(4*A(k,1)-2*A(k,2)-A(k-1,1)-A(k+1,1));
B(k,N)=A(k,N)+t*f(k,N)-beta*(4*A(k,N)-2*A(k,2)-A(k-1,N)-A(k+1,N));
end

for k=m:2*m,
B(k,1)=A(k,1)+t*f(k,1)-beta*((4+2*alfa*h)*A(k,1)-2*A(k,2)-2*h*alfa*Uext-A(k+1,1)-A(k-1,1));
B(k,N)=A(k,N);
end



for i=2:N-1,
for j=2:N-1,
B(i,j)=A(i,j)+t*f(i,j)-beta*(4*A(i,j)-A(i-1,j)-A(i+1,j)-A(i,j+1)-A(i,j-1));
 end
end

function T=tresoud(lambda,Uext,alfa,heure,N)
h=4/N;
m=ceil(N/3);
t=h*h/2;
mu=0.1;
P=ceil(heure/t);
f=zeros(N,N);
T=15*ones(N,N);
T(m:2*m,N)=lambda;

for k=1:P,
T=nextstep(T,f,Uext,alfa,N);
end

x=linspace(0,4,N);
y=x;
plot3d(x,y,T)