Learning R using a Chemical Reaction Engineering Book: Part 2

[This article was first published on Notes of a Dabbler » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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 about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

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)