Solving integral equations with fsolve

| categories: nonlinear algebra | tags: reaction engineering

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.

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)

We will need an initial guess, so we make a plot of our function to get an idea.

import matplotlib.pyplot as plt

f = np.linspace(0.01, 1)
plt.plot(f, vfunc(f))
plt.xlabel('Molar flow rate')
>>> >>> [<matplotlib.lines.Line2D object at 0x964a910>]
<matplotlib.text.Text object at 0x961fe50>

Now we can see a zero is near Fa = 0.1, so we proceed to solve the equation.

Fa_guess = 0.1
Fa_exit, = fsolve(vfunc, Fa_guess)
print 'The exit concentration is {0:1.2f} mol/L'.format(Fa_exit / nu)
>>> The exit concentration is 0.01 mol/L

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

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source

Discuss on Twitter