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 their blog: Notes of a Dabbler » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, 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.

Sponsors

Mango solutions



RStudio homepage



Zero Inflated Models and Generalized Linear Mixed Models with R

Dommino data lab

Quantide: statistical consulting and training



http://www.eoda.de







ODSC

ODSC

CRC R books series





Six Sigma Online Training





Contact us if you wish to help support R-bloggers, and place your banner here.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)