11.27

As described in Sec. PT3.1.2, linear algebraic equations can arise in the solution of differential equations. For example, the following differential equation results from a steady state mass balance for a chemical in a one-dimensional canal.

0 = D * d^c/dx^2 - U * dc/dx - kc

 where c = concentration, t = time, x = distance, D = diffusion coefficient, U = fluid velocity, and k = a first-order decay rate.Convert this differential equation to an equivalent system of simultaneous algebraic equations. Given D = 2, U = 1, k = .2, c(0) = 80 and c(10) = 20, solve these equations from x = 0 to 10 with delta_x = 2, and develop a plot of concentration versus distance.

derivatives

dc_dx = (c(x + dx) - c(x)) ./ dx ;

d2c_dx2 = (c(x + dx) - 2 * c(x) +c(x - dx)) ./ dx^2;

delta x = 1

dc_dx = (c(x + dx) - c(x)) ;

d2c_dx2 = (c(x + dx) - 2 * c(x) +c(x - dx));

0 = D * c(x + dx) - 2 * c(x) + c(x - dx) - U * (c(x + dx) - c(x)

c1 0 = D*c(2) - 2*D*c(1)+ D*c(0) - U*c(2) + U * c(1)

c2 0 = D*c(3) - 2*D*c(2)+ D*c(1) - U*c(3) + U * c(2)

c3 0 = D*c(4) - 2*D*c(3)+ D*c(2) - U*c(4) + U * c(3)

c4 0 = D*c(5) - 2*D*c(4)+ D*c(3) - U*c(5) + U * c(4) 

c5 0 = D*c(6) - 2*D*c(5)+ D*c(4) - U*c(6) + U * c(5)

c6 0 = D*c(7) - 2*D*c(6)+ D*c(5) - U*c(7) + U * c(6)

c7 0 = D*c(8) - 2*D*c(7)+ D*c(6) - U*c(8) + U * c(7)

c8 0 = D*c(9) - 2*D*c(8)+ D*c(7) - U*c(9) + U * c(8)

c9 0 = D*c(10) - 2*D*c(9)+ D*c(8) - U*c(10) + U * c(9)

c1 0 = D*c(2) - U*c(2) + U * c(1) - 2*D*c(1) + D*c(0)

c2 0 = D*c(3) - U*c(3) + U * c(2) - 2*D*c(2) + D*c(1)

c3 0 = D*c(4) - U*c(4) + U * c(3) - 2*D*c(3) + D*c(2)

c4 0 = D*c(5) - U*c(5) + U * c(4) - 2*D*c(4) + D*c(3)

c5 0 = D*c(6) - U*c(6) + U * c(5) - 2*D*c(5) + D*c(4)

c6 0 = D*c(7) - U*c(7) + U * c(6) - 2*D*c(6) + D*c(5)

c7 0 = D*c(8) - U*c(8) + U * c(7) - 2*D*c(7) + D*c(6)

c8 0 = D*c(9) - U*c(9) + U * c(8) - 2*D*c(8) + D*c(7)

c9 0 = D*c(10)- U*c(10)+ U * c(9) - 2*D*c(9) + D*c(8)

Matlab code

D = 2;

U = 1;

k = .2;

%c(0) = 80;

%c(10) = 20;

x = 0:10 ;

dx = 1;% using a delta 1

%        1       2      3     4     5    6    7     8      9

A = [ (U-2*D) (D-U)     0     0     0    0    0     0      0;...

         D   (U-2*D)  (D-U)   0     0    0    0     0      0; ...

         0      D    (U-2*D) (D-U)  0    0    0     0      0; ...

         0      0      D   (U-2*D)  (D-U) 0   0     0      0; ...

         0      0      0     D   (U-2*D)  (D-U) 0   0      0; ...

         0      0      0     0      D   (U-2*D)  (D-U) 0   0; ... 

         0      0      0     0      0   D  (U-2*D) (D-U)   0; ... 

         0      0      0     0      0   0   D  (U-2*D) (D-U); ...

         0      0      0     0      0   0   0     D  (U-2*D);];

B = [-D * 80; 0; 0; 0; 0;0 ; 0 ;0 ; -D*20 + U*20;];

C = A\B;

display(A)

display(B)

display(C)

plot(C)

Results

A =

    -3     1     0     0     0     0     0     0     0

     2    -3     1     0     0     0     0     0     0

     0     2    -3     1     0     0     0     0     0

     0     0     2    -3     1     0     0     0     0

     0     0     0     2    -3     1     0     0     0

     0     0     0     0     2    -3     1     0     0

     0     0     0     0     0     2    -3     1     0

     0     0     0     0     0     0     2    -3     1

     0     0     0     0     0     0     0     2    -3

B =

  -160

     0

     0

     0

     0

     0

     0

     0

   -20

Concentration values =

   79.9413

   79.8240

   79.5894

   79.1202

   78.1818

   76.3050

   72.5513

   65.0440

   50.0293