** Solving integral equations with fsolve
:PROPERTIES:
:categories: Nonlinear algebra
:tags: reaction engineering
:date: 2013/01/23 09:00:00
:updated: 2013/03/06 16:26:42
:END:
[[http://matlab.cheme.cmu.edu/2011/08/30/solving-integral-equations/][Original post in Matlab]]
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{1}{r_a} dFa$ 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{1}{-k Fa/\nu} dFa$$
We start by creating a function handle that describes the integrand. We can use this function in the quad command to evaluate the integral.
#+BEGIN_SRC python :session
import numpy as np
from scipy.integrate import quad
from scipy.optimize import fsolve
k = 0.23
nu = 10.0
Fao = 1.0
def integrand(Fa):
return -1.0 / (k * Fa / nu)
def func(Fa):
integral,err = quad(integrand, Fao, Fa)
return 100.0 - integral
vfunc = np.vectorize(func)
#+END_SRC
#+RESULTS:
We will need an initial guess, so we make a plot of our function to get an idea.
#+BEGIN_SRC python :session
import matplotlib.pyplot as plt
f = np.linspace(0.01, 1)
plt.plot(f, vfunc(f))
plt.xlabel('Molar flow rate')
plt.savefig('images/integral-eqn-guess.png')
plt.show()
#+END_SRC
#+RESULTS:
:
: >>> >>> []
:
[[./images/integral-eqn-guess.png]]
Now we can see a zero is near Fa = 0.1, so we proceed to solve the equation.
#+BEGIN_SRC python :session
Fa_guess = 0.1
Fa_exit, = fsolve(vfunc, Fa_guess)
print 'The exit concentration is {0:1.2f} mol/L'.format(Fa_exit / nu)
#+END_SRC
#+RESULTS:
:
: >>> The exit concentration is 0.01 mol/L
*** Summary notes
This example seemed a little easier in Matlab, where the quad function seemed to get automatically vectorized. Here we had to do it by hand.