Learning R using a Chemical Reaction Engineering Book: Part 2

January 26, 2013

(This article was first published on Notes of a Dabbler » R, and kindly contributed to R-bloggers)

In case you missed part 1, you can view it here. In this part, I tried to recreate the examples in section A.2.2 of the computational appendix in the reaction engineering book by Rawlings and Ekerdt.

Solving a nonlinear system of equations

This example involves determining reaction equilibrium conditions by solving the following system of nonlinear equations.
$\begin{aligned} PK_1y_Iy_B-y_{P1}&=&0, \\ PK_2y_Iy_B-y_{P2}&=&0 \end{aligned}$

The relation between the variables $y_I,y_B,y_{P1},y_{P2}$ and extent of reactions $x_1,x_2$ are:
$\begin{aligned} y_I&=&\frac{y_{I0}-x_1-x_2}{1-x_1-x_2} \\ y_B&=&\frac{y_{B0}-x_1-x_2}{1-x_1-x_2} \\ y_{P1}&=&\frac{y_{p10}+x_1}{1-x_1-x_2} \\ y_{P2}&=&\frac{y_{p20}+x_2}{1-x_1-x_2} \end{aligned}$

Here I have used R package rootSolve for solving the above set of equations to determine $x_1$ and $x_2$ . The library is loaded and the functions to be solved are defined in the R function fns.

# load library rootSolve
library(rootSolve)

# function defining F(x)=0
fns=function(x){
K1=108; K2=284; P=2.5
yI0=0.5; yB0=0.5; yP10=0; yP20=0;
d=1-x[1]-x[2]
yI=(yI0-x[1]-x[2])/d
yB=(yB0-x[1]-x[2])/d
yP1=(yP10+x[1])/d
yP2=(yP20+x[2])/d
F1=P*K1*yI*yB-yP1
F2=P*K2*yI*yB-yP2
c(F1=F1,F2=F2)
}

Next, an initial guess of (0.2,0.2) is set for the variables and the equations are solved using the function multiroot (from package rootSolve)

# initial guess for x
xinit=c(0.2,0.2)

# solve the equations
xans=multiroot(f=fns,start=xinit)

# object returned by multiroot
> xans
$root
[1] 0.1333569 0.3506793

$f.root
F1           F2
6.161738e-15 1.620926e-14

$iter
[1] 7

$estim.precis
[1] 1.11855e-14

# solution to the equations
> xans$root
[1] 0.1333569 0.3506793
>

The solution to the equations is accessed from the variable xans$root which in this case is (0.1334,0.3507)

MATLAB/Octave functions for solving nonlinear equations (fsolve) have been used in Chemical Engineering computations for a long time and are robust. R has traditionally not been used in this domain. So it is hard to say how the functions I have used in this blog will perform across the range of problems encountered in Reaction Engineering.