John Cagnol


Paris, France

Introduction au Calcul Scientifique
ESI CS 1402 - Printemps 2010

Schema numérique pour simuler l'équation trouvéee au paragraphe II

> a:=-62.0;
b:=73.60;
-62.0
73.60
> L:=24.0;
T:=1.0;
Nx:=240;
Nt:=2000;
hx:=L/Nx;
ht:=T/Nt;
24.0
1.0
240
2000
.1000000000
0.5000000000e-3
> ConditionInitiale:=proc(x)
if x>17.3 and x<23.55
 then return (sin(x)+1)/8+0.7
 else return 0.7;
fi;
end proc;
proc (x) if `and`(`<`(17.3, x), `<`(x, 23.55)) then return `+`(`*`(`/`(1, 8), `*`(sin(x))), `/`(1, 8), .7) else return .7 end if end proc
proc (x) if `and`(`<`(17.3, x), `<`(x, 23.55)) then return `+`(`*`(`/`(1, 8), `*`(sin(x))), `/`(1, 8), .7) else return .7 end if end proc
proc (x) if `and`(`<`(17.3, x), `<`(x, 23.55)) then return `+`(`*`(`/`(1, 8), `*`(sin(x))), `/`(1, 8), .7) else return .7 end if end proc
> ConditionBordGauche:=proc(t);
if t<0.4
 then return 0.7;
 elif t<0.7 then return 0.85;
 else return 0.7;
fi;
end proc;

ConditionBordDroit:=t->0.1;
proc (t) if `<`(t, .4) then return .7 elif `<`(t, .7) then return .85 else return .7 end if end proc
proc (t) if `<`(t, .4) then return .7 elif `<`(t, .7) then return .85 else return .7 end if end proc
proc (t) if `<`(t, .4) then return .7 elif `<`(t, .7) then return .85 else return .7 end if end proc
proc (t) options operator, arrow; .1 end proc
> for i from 0 to Nx do
 u[i,0] := ConditionInitiale(i*hx);
od:
> with(plots):
> pointplot({seq([i*hx,u[i,0]],i=0..Nx)});
Plot_2d
> for j from 0 to Nt do

 u[0,j+1] :=ConditionBordGauche(j*ht);
 u[Nx+1,j]:=ConditionBordDroit(j*ht);

 for i from 1 to Nx do
  if 2*a*u[i,j]+b>=0
  then
   if i<Nx
   then
    u[i,j+1]:=u[i,j]+ht*(2*a*u[i,j]+b)/hx*(u[i+1,j]-u[i,j]);
   else
    u[i,j+1]:=u[i,j];
   fi;
  else
   u[i,j+1]:=u[i,j]+ht*(2*a*u[i,j]+b)/hx*(u[i,j]-u[i-1,j]);
  fi;
 od;

 if j mod 100=0
  then printf("Calcul fait pour t=%f\n",j*ht);
 fi;

od:
Calcul fait pour t=0.000000
Calcul fait pour t=0.050000
Calcul fait pour t=0.100000
Calcul fait pour t=0.150000
Calcul fait pour t=0.200000
Calcul fait pour t=0.250000
Calcul fait pour t=0.300000
Calcul fait pour t=0.350000
Calcul fait pour t=0.400000
Calcul fait pour t=0.450000
Calcul fait pour t=0.500000
Calcul fait pour t=0.550000
Calcul fait pour t=0.600000
Calcul fait pour t=0.650000
Calcul fait pour t=0.700000
Calcul fait pour t=0.750000
Calcul fait pour t=0.800000
Calcul fait pour t=0.850000
Calcul fait pour t=0.900000
Calcul fait pour t=0.950000
Calcul fait pour t=1.000000
> for j from 0 to Nt do
 P[j]:=pointplot({seq([i*hx,u[i,j]],i=0..Nx)});
od:

display([seq(P[i],i=0..Nt)],insequence=true);
Plot_2d