/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 15.08.2 ] */

/* [wxMaxima: comment start ]
Fct utiles
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
DitCeQueCest([s]):=apply(print,append(s,[%]))$DCQ:DitCeQueCest$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Advection
   [wxMaxima: comment end   ] */

/* [wxMaxima: section start ]
Diffusion
   [wxMaxima: section end   ] */

/* [wxMaxima: comment start ]
Equation différentielle ordinaire en F
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
Vu:u=x/(2*sqrt(D*t))$DCQ("u vs x et t : ")$
depends(F,u)$DCQ("F depend de u")$
gradef(u,x,diff(ev(u,Vu),x))$DCQ("du/dx")$
gradef(u,t,diff(ev(u,Vu),t))$DCQ("du/dt")$
rho:F$DCQ("rho=F")$
rhoxx:diff(rho,x,2)$DCQ("d^2rho/dx^2 = ")$
rhot:diff(rho,t)$DCQ("drho/dt = ")$
doe:factor(ev(D*rhoxx-rhot,solve(Vu,x)))$DCQ("EDO")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Résolution
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
soldoe:ode2(doe,F,u)$DCQ("Resolution EDO")$
eqc:[ev(F,eval,soldoe,u=0)=rhos,limit(ev(F,soldoe),u,inf)=0]$DCQ("conditions : ")$
solc:solve(eqc,[%k1,%k2]);DCQ("solution")$
rho:sol:factor(ev(ev(F,soldoe),solc,eval,Vu))$DCQ("rho = ")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Exploitation
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
plot2d(ev(rho/rhos,D=1.0e-9,t=3600),[x,0,0.1])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
plot2d([ev(rho/rhos,D=1.0e-9,x=0.1,t=24*3600*t),30/360],[t,0,25])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Régime stationnaire
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
rhop:%k1*x+%k2$DCQ("forme de base : ")$
eq:[ev(rhop,x=0)=rhos,ev(-D*diff(rhop,x)-alpha*rhop,eval,x=L)]$DCQ("Eqs : ")$
sol:solve(eq,[%k1,%k2])$DCQ("Sols : ")$
rhop:ev(rhop,sol)$DCQ("rhop = ")$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: section start ]
Propagation
   [wxMaxima: section end   ] */

/* [wxMaxima: input   start ] */
u0(t):=if ((t<0) or (t>T)) then 0 else (t*(T-t))*4/T^2$ 
plot2d(ev(u0(t),T=1),[t,-3,3],[y,0,1.5])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
u0(t):=if ((t<0) or (t>T)) then 0 else 
 if (t<T/4) then (t*(T/4-t))*4/(T/4)^2 
 else if (t>3/4*T) then -((t-3*T/4)*(T-t))*4/(T/4)^2 else 0$
plot2d(ev(u0(t),T=1),[t,-3,3],[y,-1.5,1.5])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
u0(t):=if t<0 then 0 else sin(2*%pi*t/T)$
plot2d(ev(u0(t),T=1),[t,-3,3],[y,-1.5,1.5])$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
u(t,x=L)=0
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
block([T:0.3,L:1,c:1,N:20,ts:0.1],
    for i:0 thru N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(u0(t-x/c),[x,0,1],[y,-1.5,1.5],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
block([T:0.1,L:1,c:1,N:20,ts:0.1],
    for i:0 thru 2*N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(
             u0(t-x/c)-u0(t+x/c-2*L/c)
            ,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
block([T:0.1,L:1,c:1,N:100,ts:0.1],
    for i:0 thru 3*N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(
             u0(t-x/c)-u0(t+x/c-2*L/c)+u0(t-x/c-2*L/c)
            ,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
block([T:0.3,L:1,c:1,N:100,ts:0.01],
    for i:0 thru 4*N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(
             u0(t-x/c)-u0(t+x/c-2*L/c)+u0(t-x/c-2*L/c)-u0(t+x/c-4*L/c)
            ,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
block([T:0.1,L:1,c:1,N:100,ts:0.01,P:10],
 ut:apply("+",makelist('u0(t-x/c-2*i*L/c)-'u0(t+x/c-2*(i+1)*L/c),i,0,P)),
 for i:0 thru 2*(P+1)*N+1 do
   block([t:i/N*L],?sleep(ts),plot2d(ut,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
makelist('u0(t-x/c-2*i*L/c)-'u0(t+x/c-2*(i+1)*L/c),i,0,0);
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
block([T:0.3,L:1,c:1,N:100,ts:0.01,P:10],
 ut:flatten(makelist(['u0(t-x/c-2*i*L/c),-'u0(t+x/c-2*(i+1)*L/c)],i,0,P)),
 for i:0 thru 2*(P+1)*N+1 do
   block([t:i/N*L],?sleep(ts),plot2d(ut,[x,0,1],[y,-2.1,2.1],[ylabel,"u"],[legend,false])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
du/dx(t,x=L)=0
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
block([T:0.4,L:1,c:1,N:100,ts:0.1],
    for i:0 thru 4*N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(
             u0(t-x/c)+u0(t+x/c-2*L/c)-u0(t-x/c-2*L/c)-u0(t+x/c-4*L/c)
            ,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
du/dx(t,x=L)= -c du/dx(t,x=L)
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
block([T:0.4,L:1,c:1,N:100,ts:0.1],
    for i:0 thru 4*N+1 do
      block([t:i/N*L],?sleep(ts),plot2d(
             u0(t-x/c)
            ,[x,0,1],[y,-2.1,2.1],[ylabel,"u"])))$
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$