Solving an ode for a specific solution value

August 31, 2011 at 08:01 PM | categories: odes | View Comments

Solving an ode for a specific solution value

Solving an ode for a specific solution value

John Kitchin 8/31/2011

The analytical solution to an ODE is a function, which can be solved to get a particular value, e.g. if the solution to an ODE is y(x) = exp(x), you can solve the solution to find the value of x that makes $y(x)=2$. In a numerical solution to an ODE we get a vector of independent variable values, and the corresponding function values at those values. To solve for a particular function value we need a different approach. This post will show one way to do that in Matlab.

Contents

The problem

given that the concentration of a species A in a constant volume, batch reactor obeys this differential equation $\frac{dC_A}{dt}=- k C_A^2$ with the initial condition $C_A(t=0) = 2.3$ mol/L and $k = 0.23$ L/mol/s, compute the time it takes for $C_A$ to be reduced to 1 mol/L.

Solution

we will use Matlab to integrate the ODE to obtain the solution, and then use the deval command in conjuction with fsolve to get the numerical solution we are after.

clear all; close all; clc

k = 0.23;
Cao = 2.3;
dCadt = @(t,Ca) -k*Ca^2;
tspan = [0 10];
sol = ode45(dCadt,tspan,Cao)
% sol here is a struct with fields corresponding to the solution and other
% outputs of the solver.

sol = 

     solver: 'ode45'
    extdata: [1x1 struct]
          x: [1x12 double]
          y: [1x12 double]
      stats: [1x1 struct]
      idata: [1x1 struct]

We need to make sure the solution vector from the ode solver contains the answer we need

plot(sol.x,sol.y)
xlabel('Time (s)')
ylabel('C_A (mol/L)')
% the concentration of A is at approximately t=3

We can make a function Ca(t) using the solution from the ode solver and the deval command, which evaluates the ODE solution at a new value of t.

Ca = @(t) deval(sol,t);

% now we create a function for fsolve to work on. remember this function is
% equal to zero at the solution.
tguess = 3;
f = @(t) 1 - Ca(t) ;
t_solution = fsolve(f,tguess)

Ca(t_solution) % this should be 1.0

sprintf('C_A = 1 mol/L at a time of %1.2f seconds.',t_solution)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




t_solution =

    2.4598


ans =

    1.0000


ans =

C_A = 1 mol/L at a time of 2.46 seconds.

A related post is Post 954 on solving integral equations.

'done'

% categories: ODEs
% tags: reaction engineering
% post_id = 968; %delete this line to force new post;

ans =

done


Published with MATLAB® 7.12

Read and Post Comments

Solving integral equations

August 31, 2011 at 12:29 AM | categories: nonlinear algebra | View Comments

Solving integral equations

Solving integral equations

John Kitchin

Contents

Occasionally we have integral equations we need to solve in engineering problems, for example, the volume of plug flow reactor can be defined by this equation: $V = \int_{Fa(V=0)}^{Fa} \frac{dFa}{r_a}$ where $r_a$ is the rate law. Suppose we know the reactor volume is 100 L, the inlet molar flow of A is 1 mol/L, the volumetric flow is 10 L/min, and $r_a = -k Ca$, with $k=0.23$ 1/min. What is the exit molar flow rate? We need to solve the following equation:

$$100 = \int_{Fa(V=0)}^{Fa} \frac{dFa}{-k Fa/\nu}$$

k = 0.23;
nu = 10;
Fao = 1;

We start by creating a function handle that describes the integrand. We can use this function in the quad command to evaluate the integral.

f1 = @(Fa) -1./(k*Fa/nu);

Now we create a second function handle that defines the equation to solve. We use the integrand function handle above in the quad command to evaluate the integral

f2 = @(Fa) 100 - quad(f1,Fao,Fa);

plot f2 to get an idea of where the solution is.

the maximum that Fa can be is Fao, and the minimum is zero. From the graph you can see the solution is at approximately Fa = 0.1.

ezplot(f2,0,Fao)

Warning: Function failed to evaluate on array inputs; vectorizing the function may
speed up its evaluation and avoid the need to loop over array elements.

Get the solution

Fa_guess = 0.1;
Fa_exit = fsolve(f2, Fa_guess);

Ca_exit = Fa_exit/nu;
sprintf('The exit concentration is %1.2f mol/L',Ca_exit)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.




ans =

The exit concentration is 0.01 mol/L


% categories: nonlinear algebra
% tags: reaction engineering
% post_id = 954; %delete this line to force new post;


Published with MATLAB® 7.12

Read and Post Comments

Constrained minimization to find equilibrium compositions

August 12, 2011 at 02:03 PM | categories: optimization | View Comments

Constrained minimization to find equilibrium compositions

Constrained minimization to find equilibrium compositions

adapated from Chemical Reactor analysis and design fundamentals, Rawlings and Ekerdt, appendix A.2.3.

Contents

The equilibrium composition of a reaction is the one that minimizes the total Gibbs free energy. The Gibbs free energy of a reacting ideal gas mixture depends on the mole fractions of each species, which are determined by the initial mole fractions of each species, the extent of reactions that convert each species, and the equilibrium constants.

Reaction 1: $I + B \rightleftharpoons P1$

Reaction 2: $I + B \rightleftharpoons P2$

The equilibrium constants for these reactions are known, and we seek to find the equilibrium reaction extents so we can determine equilibrium compositions.

function main

clc, close all, clear all

constraints

we have the following constraints, written in standard less than or equal to form:

$-\epsilon_1 \le 0$

$-\epsilon_2 \le 0$

$\epsilon_1 + \epsilon_2 \le 0.5$

We express this in matrix form as Ax=b where $A = \left[ \begin{array}{cc} -1 & 0 \\ 0 & -1 \\ 1 & 1 \end{array} \right]$ and $b = \left[ \begin{array}{c} 0 \\ 0 \\ 0.5\end{array} \right]$

A = [[-1 0];[0 -1];[1 1]];
b = [0;0;0.5];

% this is our initial guess
X0 = [0.1 0.1];

This does not work!

X = fmincon(@gibbs,X0,A,b)

Warning: Trust-region-reflective algorithm does not solve this type of
problem, using active-set algorithm. You could also try the interior-point or
sqp algorithms: set the Algorithm option to 'interior-point' or 'sqp' and
rerun. For more help, see Choosing the Algorithm in the documentation. 

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunEvals = 200 (the default value).


X =

      NaN +    NaNi      NaN +    NaNi

the error message suggests we try another solver. Here is how you do that.

options = optimset('Algorithm','interior-point');
X = fmincon(@gibbs,X0,A,b,[],[],[],[],[],options)
% note the options go at the end, so we have to put a lot of empty [] as
% placeholders for other inputs to fmincon.

% we can see what the gibbs energy is at our solution
gibbs(X)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the function tolerance,
and constraints were satisfied to within the default value of the constraint tolerance.




X =

    0.1334    0.3507


ans =

   -2.5594

getting the equilibrium compositions

these are simply mole balances to account for the consumption and production of each species

e1 = X(1);
e2 = X(2);

yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0; %initial mole fractions

d = 1 - e1 - e2;
yI = (yI0 - e1 - e2)/d;
yB = (yB0 - e1 - e2)/d;
yP1 = (yP10 + e1)/d;
yP2 = (yP20 + e2)/d;

sprintf('y_I = %1.3f y_B = %1.3f y_P1 = %1.3f y_P2 = %1.3f',yI,yB,yP1,yP2)

ans =

y_I = 0.031 y_B = 0.031 y_P1 = 0.258 y_P2 = 0.680

Verifying our solution

One way we can verify our solution is to plot the gibbs function and see where the minimum is, and whether there is more than one minimum. We start by making grids over the range of 0 to 0.5. Note we actually start slightly above zero because at zero there are some numerical imaginary elements of the gibbs function or it is numerically not defined since there are logs of zero there.

[E1, E2] = meshgrid(linspace(0.001,0.5),linspace(0.001,0.5));

% we have to remember that when e1 + e2 > 0.5, there is no solution to our
% gibbs functions, so we find the indices of the matrix where the sum is
% greater than 0.5, then we set those indices to 0 so no calculation gets
% done for them.
sumE = E1 + E2;
ind = sumE > 0.5;
E1(ind) = 0; E2(ind) = 0;

% we are going to use the arrayfun function to compute gibbs at each e1 and
% e2, but our gibbs function takes a vector input not two scalar inputs. So
% we create a function handle to call our gibbs function with two scalar
% values as a vector.
f = @(x,y) gibbs([x y]);

% finally, we compute gibbs at each reaction extent
G = arrayfun(f, E1, E2);

% and make a contour plot at the 100 levels
figure; hold all;
contour(E1,E2,G,100)
xlabel('\epsilon_1')
ylabel('\epsilon_2')

% let's add our solved solution to the contour plot
plot([X(1)],[X(2)],'ko','markerfacecolor','k')
colorbar
% you can see from this contour plot that there is only one minimum.

'done'

ans =

done


function retval = gibbs(X)
% we do not derive the form of this function here. See chapter 3.5 in the
% referenced text book.
e1 = X(1);
e2 = X(2);

K1 = 108; K2 = 284; P = 2.5;
yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0;

d = 1 - e1 - e2;
yI = (yI0 - e1 - e2)/d;
yB = (yB0 - e1 - e2)/d;
yP1 = (yP10 + e1)/d;
yP2 = (yP20 + e2)/d;

retval = -(e1*log(K1) + e2*log(K2)) + ...
    d*log(P) + yI*d*log(yI) + ...
    yB*d*log(yB) + yP1*d*log(yP1) + yP2*d*log(yP2);

% categories: Optimization
% tags: thermodynamics, reaction engineering


Published with MATLAB® 7.10

Read and Post Comments

first order reversible reaction in batch reactor

August 07, 2011 at 06:56 PM | categories: odes | View Comments

first order reversible reaction in batch reactor

first order reversible reaction in batch reactor

John Kitchin

Contents

Problem statement

$A \rightleftharpoons B$

forward rate law: $-r_a = k_1 C_A$

backward rate law: $-r_b = k_{-1} C_B$

this example illustrates a set of coupled first order ODES

function main

clc; clear all; close all;

u = cmu.units; %note that using units makes this m-file run pretty slowly
tspan = [0 5]*u.min;
init = [1 0]*u.mol/u.L;

[t,C] = ode45(@myode,tspan,init);

plot the solution

plot(t/u.min,C/(u.mol/u.L))
xlabel('Time (min)')
ylabel('Concentration (mol/L)')
legend C_A C_B

you need this to make the plot appear in the right place above when published. It appears you need a final cell with an output to avoid having the plot showup at the end of the published document.

disp('end')

end

function dCdt = myode(t,C)
% ra = -k1*Ca
% rb = -k_1*Cb
% net rate for production of A:  ra - rb
% net rate for production of B: -ra + rb
u = cmu.units;
k1 = 1/u.min;
k_1 = 0.5/u.min;

Ca = C(1);
Cb = C(2);

ra = -k1*Ca; rb = -k_1*Cb;

dCadt = ra - rb;
dCbdt = -ra + rb;

dCdt = [dCadt; dCbdt];

% categories: ODEs
% tags: reaction engineering


% post_id = 706; %delete this line to force new post;


Published with MATLAB® 7.10

Read and Post Comments

Determining linear independence of a set of vectors

August 02, 2011 at 02:48 PM | categories: linear algebra | View Comments

Determining linear independence of a set of vectors

Determining linear independence of a set of vectors

Occasionally we have a set of vectors and we need to determine whether the vectors are linearly independent of each other. This may be necessary to determine if the vectors form a basis, or to determine how many independent equations there are, or to determine how many independent reactions there are.

Reference: Kreysig, Advanced Engineering Mathematics, sec. 7.4

Contents

Problem set 7.4 - #15

v1 = [6 0 3 1 4 2];
v2 = [0 -1 2 7 0 5];
v3 = [12 3 0 -19 8 -11];

R = rank([v1; v2; v3])
[rows, ~] = size([v1; v2; v3])
% note the ~ in the output indicates we do not care what the value is.

% the number of rows is greater than the rank, so these vectors are not
% independent. Let's demonstrate that one vector can be defined as a linear
% combination of the other two vectors. Mathematically we represent this
% as:
%
% $x_1 \mathit{v1} + x_2 \mathit{v2} = v3%
%
% or
%
% [x_1 x_2][v1; v2] = v3
%
% this is not the usual linear algebra form of Ax = b. To get there, we
% transpose each side of the equation to get:
%
% [v1.' v2.'][x_1; x_2] = v3'
%
% which is the form Ax = b. In Matlab we solve that as x = A\b (

R =

     2


rows =

     3


x = [v1.' v2.']\v3.'

x(1)*v1 + x(2)*v2
v3
% you can see that v3 = 2v1 - 3v2, illustrating that v3 is not linearly
% independent of v1 and v2.

x =

    2.0000
   -3.0000


ans =

   12.0000    3.0000    0.0000  -19.0000    8.0000  -11.0000


v3 =

    12     3     0   -19     8   -11

Problem set 7.4 #17

v1 = [0.2 1.2 5.3 2.8 1.6];
v2 = [4.3 3.4 0.9 2.0 -4.3];

rank([v1; v2])
% the rank is equal to the number of rows, so these vectors are linearly
% independent. You could also see that by inspection since the signs of the
% last element are different. there is no way to convert v1 to v2 by simple
% scaling.

ans =

     2

Near deficient rank

the rank command roughly works in the following way: the matrix is converted to a reduced row echelon form, and then the number of rows that are not all equal to zero are counted. Matlab uses a tolerance to determine what is equal to zero. If there is uncertainty in the numbers, you may have to define what zero is, e.g. if the absolute value of a number is less than 1e-5, you may consider that close enough to be zero. The default tolerance is usually very small, of order 1e-15.

A = [[1 2 3];
    [0 2 3];
    [0 0 1e-6]];

rank(A)
% matlab considers this matrix to have a full rank of 3 because the default
% tolerance in this case is 2.7e-15.

ans =

     3

rank with tolerance

If we believe that any number less than 1e-5 is practically equivalent to zero, we can use that information to compute the rank like this.

tolerance = 1e-5;
rank(A,tolerance)
% now the A matrix has a rank of only 2.

ans =

     2

Application to independent chemical reactions.

reference: Exercise 2.4 in Chemical Reactor Analysis and Design Fundamentals by Rawlings and Ekerdt.

The following reactions are proposed in the hydrogenation of bromine:

Let this be our species vector: v = [H2 H Br2 Br HBr]^T

the reactions are then defined by M*v where M is a stoichometric matrix in which each row represents a reaction with negative stoichiometric coefficients for reactants, and positive stoichiometric coefficients for products. A stoichiometric coefficient of 0 is used for species not participating in the reaction.

%    [H2  H Br2 Br HBr]
M = [[-1  0 -1  0  2];  % H2 + Br2 == 2HBR
     [ 0  0 -1  2  0];  % Br2 == 2Br
     [-1  1  0 -1  1];  % Br + H2 == HBr + H
     [ 0 -1 -1  1  1];  % H + Br2 == HBr + Br
     [ 1 -1  0  1  -1];  % H + HBr == H2 + Br
     [ 0  0  1 -2  0]]; % 2Br == Br2

rank(M)

% 6 reactions are given, but the rank of the matrix is only 3. so there
% are only four independent reactions. You can see that reaction 6 is just
% the opposite of reaction 2, so it is clearly not independent. Also,
% reactions 3 and 5 are just the reverse of each other, so one of them can
% also be eliminated. finally, reaction 4 is equal to reaction 1 minus
% reaction 3.

ans =

     3

choosing independent reactions

We can identify independent reactions by examining the reduced row echelon form of the matrix where the reactions are in the columns rather than rows. That is simply the transpose of the M matrix above. The columns with leading ones correspond to the reactions that can form a basis, i.e. the independent reactions.

rref(M')
% this shows columns 1,2,3 have leading ones, indicating the corresponding
% reactions are a suitable basis for all the reactions. If that is true,
% then we should be able to find solutions to the following linear
% equation: basis*x_i = reaction_i

ans =

     1     0     0     1     0     0
     0     1     0     0     0    -1
     0     0     1    -1    -1     0
     0     0     0     0     0     0
     0     0     0     0     0     0


basis = transpose(M(1:3,:)); % this is rows 1-3 of the M matrix

reaction_4 = M(4,:).';
reaction_5 = M(5,:).';
reaction_6 = M(6,:).';

x_4 = basis\reaction_4
x_5 = basis\reaction_5
x_6 = basis\reaction_6

% x_4, x_5 and x_6 are all integers, which means reactions 4, 5 and 6 are
% linear combinations of reactions 1, 2 and 3. As pointed out above, x_4
% shows reaction 4 to be reaction 1 - reaction 3; x_5 shows reaction 5 is
% -reaction 3; and x_6 shows reaction 6 is minus reaction 4.

x_4 =

    1.0000
    0.0000
   -1.0000


x_5 =

   -0.0000
    0.0000
   -1.0000


x_6 =

   -0.0000
   -1.0000
    0.0000

lets also demonstrate that reaction 1 is not a linear combination of reactions 2 and 3.

reaction_1 = M(1,:).';
x_1 = basis\reaction_1
% you can see the coefficients for reactions 2 and 3 are equal to zero,
% indicating reaction 1 is linearly independent of reactions 2 and 3.

x_1 =

    1.0000
    0.0000
    0.0000

Related posts

transposition rules

solving linear equations

% categories: Linear algebra
% tags: math, reaction engineering

% post_id = 640; %delete this line to force new post;


Published with MATLAB® 7.10

Read and Post Comments
« Previous Page