Diffusion and Reaction in a One Dimensional Slab
> restart:
This example involves the concentration in plug flow reactor with axial dispersion. The governing ODE is
> ODE := diff(C[A](z),z,z)=k/D[A,B]*C[A](z): ODE;
where is the diffusion coefficient and is the concentration. The boundary conditions for this example are written in the following form.
> BC0 := C[A](0)=C[A,0]: BC0;
> BC1:= D(C[A])(L)=0: BC1;
Method 1
This problem has an analytical solution that is easily obtained with Maple. (It's easy with R4. If you are using R5 then you need to execute the follwoing read statement to load Joe Riel's patch to dsolve, which, otherwise, cannot handle indexed variables in this context!!)
> read `e:/maple/utils/dsolve.mtx`:
> dsolve({ODE, BC0,BC1},C[A](z));
We can express this in the more familir form involving hyperbolic functions as follows
> convert(%,trig);
Method 2
We shall also attempt a numerical solution. The problem is not straightforward because the ODE to be solved is second order and we will need an estimate of the first derivative in order to use the shooting method as described in the Polymath solution.
We define as the first derivative term.
> ODE3:=diff(C[A](z),z)=y(z): ODE3;
The second derivative is obtained on differentiation.
> diff(ODE3,z);
Substitution of the original ODE gives the following expression for the ODE
> ODE4:=subs(ODE,%): ODE4;
The parameter values are specified
> Params :={D[A,B]=1.2e-9,k=1e-3};
and substituted into the set of differential equations.
> Deqns := subs(Params,{ODE3,ODE4});
The initial conditions are the known initial concentration and a guess for the unknown first derivative term.
> IC := {C[A](0)=0.2,y(0)=-130};
We combine the sets of ODEs and initial conditions.
> Alleqns := Deqns union IC;
dsolve is invoked to make a procedure for numerical computation of the result.
> g:=dsolve(Alleqns,{C[A](z),y(z)}, type=numeric, method=rkf45, output=procedurelist);
We compute the "solution" (subject to the guessed initial value of y) as follows
> g(1e-3);
where the value of 0.001 in the call to g is the final value of .
We need to find an initial value of y that leads to (more or less). We make a procedure to automate the preceding computations:
> f := proc(y0) local IC, Deqns, Params, Alleqns,g, err; Params :={D[A,B]=1.2e-9,k=1e-3, y[f,0]=0}; Deqns := {833333.3333*C[A](z) = diff(y(z),z), diff(C[A](z),z) = y(z)}; IC := {C[A](0)=0.2,y(0)=y0}; Alleqns := Deqns union IC; g:=dsolve(Alleqns,{C[A](z),y(z)}, type=numeric, method=rkf45, output=procedurelist); err := subs(Params,rhs(g(1e-3)[2]) - y[f,0]); end:
and test it for some differing initial values
> f(-130);
> f(-140);
The different signs obtained above is an indication that there is a solution to this problem at a value of somewhere between -130 and -140. In Maple Releave 4 we could ask fsolve to compute a numerical solution to this equation using the command shown below.
> # result:=fsolve('f'(y0),y0=-140..-130);
However, Release 5 is unable to solve this problem quite so easily, and we shall resort to using our own implementation of Newton's method. First we must rewrite the above proc so that it returns a list,
rather than a single number.
> f := proc(y0) local IC, Deqns, Params, Alleqns,g, err; Params :={D[A,B]=1.2e-9,k=1e-3, y[f,0]=0}; Deqns := {833333.3333*C[A](z) = diff(y(z),z), diff(C[A](z),z) = y(z)}; IC := {C[A](0)=0.2,y(0)=y0}; Alleqns := Deqns union IC; g:=dsolve(Alleqns,{C[A](z),y(z)}, type=numeric, method=rkf45, output=procedurelist); err := [subs(Params,rhs(g(1e-3)[2]) - y[f,0])]; end:
Now we can invoke Newton's method.
> read `e:/maple/numerics/nonlin/newton.mpl`:
> result:=Newton(f,y0=-135): result;
The value obtained above is the "correct" value of the first derivative at the origin. We need to compute all the other variables at the end using this value for the first derivative. To this end we repeat our initial computations, first remaking the initial condition.
> IC := {C[A](0)=0.2,y(0)=op(2,result)};
then reforming the full set of differential equations and initial conditions,
and asking dsolve to generate a procedure from which we may compute a numerical solution
Finally, we invoke the procedure to find the desired concentration.
Method 3
Here, we are going to use a finite difference method of solving the ODE. Our first step is to stop using functional notation in the original equation
> ODE2:=subs(diff=Diff,C[A](z)=C[A],ODE): ODE2;
as well as in the boundary conditions
> BC1:= Diff(C[A](L),z)=0: BC1;
The governing equation is assumed to hold for the concentration at certain selected points (nodes), here indicated by the index i.
> i:='i': ODE3:=subs(C[A]=C[A,i],ODE2): ODE3;
We set the number of nodes
> n:=11;
We use a central difference approximation to the second derivative in the above equation
> fddiff:=Diff(Diff(C[A,i],z),z) = (C[A,i+1]-2*C[A,i]+C[A,i-1])/h^2: fddiff;
We can now generate finite difference approximations for each interior node (the boundaries are dealt with separately below).
> for i from 2 to n-1 do eqn[i] :=subs(fddiff,ODE3); print(eqn[i]); od:
The boundary condition at is
> BC0a:=subs(C[A](0)=C[A,1],BC0): BC0a;
where is the initial concentration. The derivative boundary condition at the other boundary must be replaced using a 2 term backward difference approximation
> fddiffb2:=Diff(C[A,i],z)=1/2*(3*C[A,i]-4*C[A,i-1]+C[A,i-2])/h: fddiffb2;
The second boundary condition now becomes
> BC1a:=subs(C[A](L)=C[A,i],fddiffb2,i=n,BC1): BC1a;
We place all the equations into a list.
> eqnlist :=[BC0a,seq(eqn[i],i=2..n-1),BC1a];
In order to obtain a numerical solution we must assign values to the diffusion coefficient, reaction rate coefficient and so on. The values used here are as follows.
> Params := {D[A,B]=1.2e-9,k=1e-3,h=1e-3/(n-1),C[A,0]=0.2};
The step size, , is determined by the fact that the distance over which diffusion and reaction take place is just m.
We substitute these variables into the finite difference equations.
> eqnlist2:=subs(Params,eqnlist);
The variables appearing in these equations are
> vars:=indets(eqnlist2,name);
and they number
> nops(vars);
The number of equations is
> nops(eqnlist2);
So there appear to be no further degrees of freedom. Before proceeding it might help to determine the structure of the equations. Here is a procedure for looking at the structure of equation sytems.
> Eqnanalysis := proc (Eqns::{set,list,vector},x::{set,list,vector}) local neqns, Imat, i, j; # Determine number of equations neqns := nops(Eqns); # Now to construct the incidence matrix Imat := linalg[matrix](neqns,neqns,0); for i from 1 to neqns do for j from 1 to neqns do if has (Eqns[i], x[j]) then Imat[i,j] := 1; fi; od; od; plots[sparsematrixplot](Imat); end:
> Eqnanalysis(eqnlist2,[seq(C[A,i],i=1..n)]);
Here we see that the system of equations has an almost tridiagonal form. We could, had we a mind to do so, take advantage of this structure in devising a suitable computational method for calculating the concentrations. In this case there is no need to do so since Maple can solve these equations directly.
> result:=solve({op(eqnlist2)},vars): result;
Maple can also solve the purely symbolic problem in which numerical values have not been given to the model parameters.
> solve({op(eqnlist)},vars);
We can check that this gives the same result as before on substituting the parameter values
> subs(Params,%);
As an alternative - not needed in this case but necessary in others - we could use Newton's method to solve the equations numerically.
A list of initial guesses is created.
> Initial:=[seq(C[A,i]=0.2,i=1..n)];
We invoke Newton's method as follows:
> result:=Newton(eqnlist2,Initial,tolerance=1e-6,output={variables});
An advantage of this solution is that the concentrations at each node are listed in the original order.
The concentration profile can be plotted as shown in the next command.
> plot([seq([k,rhs(result[k])],k=1..n)],labels=['n','C']);
Note the profile flattens out on the right hand side of the plot, as is to be expected from a boundary condition that forces to zero the concentration derivative.
Inspection of the iteration history shows that the solution was obtained in a single step (round off is responsible for the procedure taking the extra step). This means that the equations are linear (we could have deduced this fact if only we had bothered to look).
>
