Learning R using a Chemical Reaction Engineering Book: Part 2

January 26, 2013
By

(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.


To leave a comment for the author, please follow the link and comment on his blog: Notes of a Dabbler » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.