The Gibbs free energy of a reacting mixture and the equilibrium composition

| categories: optimization | tags: reaction engineering, thermodynamics | View Comments

Table of Contents

Matlab post

In this post we derive the equations needed to find the equilibrium composition of a reacting mixture. We use the method of direct minimization of the Gibbs free energy of the reacting mixture.

The Gibbs free energy of a mixture is defined as \(G = \sum\limits_j \mu_j n_j\) where \(\mu_j\) is the chemical potential of species \(j\), and it is temperature and pressure dependent, and \(n_j\) is the number of moles of species \(j\).

We define the chemical potential as \(\mu_j = G_j^\circ + RT\ln a_j\), where \(G_j^\circ\) is the Gibbs energy in a standard state, and \(a_j\) is the activity of species \(j\) if the pressure and temperature are not at standard state conditions.

If a reaction is occurring, then the number of moles of each species are related to each other through the reaction extent \(\epsilon\) and stoichiometric coefficients: \(n_j = n_{j0} + \nu_j \epsilon\). Note that the reaction extent has units of moles.

Combining these three equations and expanding the terms leads to:

$$G = \sum\limits_j n_{j0}G_j^\circ +\sum\limits_j \nu_j G_j^\circ \epsilon +RT\sum\limits_j(n_{j0} + \nu_j\epsilon)\ln a_j $$

The first term is simply the initial Gibbs free energy that is present before any reaction begins, and it is a constant. It is difficult to evaluate, so we will move it to the left side of the equation in the next step, because it does not matter what its value is since it is a constant. The second term is related to the Gibbs free energy of reaction: \(\Delta_rG = \sum\limits_j \nu_j G_j^\circ\). With these observations we rewrite the equation as:

$$G - \sum\limits_j n_{j0}G_j^\circ = \Delta_rG \epsilon +RT\sum\limits_j(n_{j0} + \nu_j\epsilon)\ln a_j $$

Now, we have an equation that allows us to compute the change in Gibbs free energy as a function of the reaction extent, initial number of moles of each species, and the activities of each species. This difference in Gibbs free energy has no natural scale, and depends on the size of the system, i.e. on \(n_{j0}\). It is desirable to avoid this, so we now rescale the equation by the total initial moles present, \(n_{T0}\) and define a new variable \(\epsilon' = \epsilon/n_{T0}\), which is dimensionless. This leads to:

$$ \frac{G - \sum\limits_j n_{j0}G_j^\circ}{n_{T0}} = \Delta_rG \epsilon' + RT \sum\limits_j(y_{j0} + \nu_j\epsilon')\ln a_j $$

where \(y_{j0}\) is the initial mole fraction of species \(j\) present. The mole fractions are intensive properties that do not depend on the system size. Finally, we need to address \(a_j\). For an ideal gas, we know that \(A_j = \frac{y_j P}{P^\circ}\), where the numerator is the partial pressure of species \(j\) computed from the mole fraction of species \(j\) times the total pressure. To get the mole fraction we note:

$$y_j = \frac{n_j}{n_T} = \frac{n_{j0} + \nu_j \epsilon}{n_{T0} + \epsilon \sum\limits_j \nu_j} = \frac{y_{j0} + \nu_j \epsilon'}{1 + \epsilon'\sum\limits_j \nu_j} $$

This finally leads us to an equation that we can evaluate as a function of reaction extent:

$$ \frac{G - \sum\limits_j n_{j0}G_j^\circ}{n_{T0}} = \widetilde{\widetilde{G}} = \Delta_rG \epsilon' + RT\sum\limits_j(y_{j0} + \nu_j\epsilon') \ln\left(\frac{y_{j0}+\nu_j\epsilon'}{1+\epsilon'\sum\limits_j\nu_j} \frac{P}{P^\circ}\right) $$

we use a double tilde notation to distinguish this quantity from the quantity derived by Rawlings and Ekerdt which is further normalized by a factor of \(RT\). This additional scaling makes the quantities dimensionless, and makes the quantity have a magnitude of order unity, but otherwise has no effect on the shape of the graph.

Finally, if we know the initial mole fractions, the initial total pressure, the Gibbs energy of reaction, and the stoichiometric coefficients, we can plot the scaled reacting mixture energy as a function of reaction extent. At equilibrium, this energy will be a minimum. We consider the example in Rawlings and Ekerdt where isobutane (I) reacts with 1-butene (B) to form 2,2,3-trimethylpentane (P). The reaction occurs at a total pressure of 2.5 atm at 400K, with equal molar amounts of I and B. The standard Gibbs free energy of reaction at 400K is -3.72 kcal/mol. Compute the equilibrium composition.

import numpy as np

R = 8.314
P = 250000  # Pa
P0 = 100000 # Pa, approximately 1 atm
T = 400 # K

Grxn = -15564.0 #J/mol
yi0 = 0.5; yb0 = 0.5; yp0 = 0.0; # initial mole fractions

yj0 = np.array([yi0, yb0, yp0])
nu_j = np.array([-1.0, -1.0, 1.0])   # stoichiometric coefficients

def Gwigglewiggle(extentp):
    diffg = Grxn * extentp
    sum_nu_j = np.sum(nu_j)
    for i,y in enumerate(yj0):
        x1 = yj0[i] + nu_j[i] * extentp
        x2 = x1 / (1.0 + extentp*sum_nu_j)
        diffg += R * T * x1 * np.log(x2 * P / P0)
    return diffg

There are bounds on how large \(\epsilon'\) can be. Recall that \(n_j = n_{j0} + \nu_j \epsilon\), and that \(n_j \ge 0\). Thus, \(\epsilon_{max} = -n_{j0}/\nu_j\), and the maximum value that \(\epsilon'\) can have is therefore \(-y_{j0}/\nu_j\) where \(y_{j0}>0\). When there are multiple species, you need the smallest \(epsilon'_{max}\) to avoid getting negative mole numbers.

epsilonp_max = min(-yj0[yj0 > 0] / nu_j[yj0 > 0])
epsilonp = np.linspace(1e-6, epsilonp_max, 1000);

import matplotlib.pyplot as plt

plt.plot(epsilonp,Gwigglewiggle(epsilonp))
plt.xlabel('$\epsilon$')
plt.ylabel('Gwigglewiggle')
plt.savefig('images/gibbs-minim-1.png')
>>> >>> >>> __main__:7: RuntimeWarning: divide by zero encountered in log
__main__:7: RuntimeWarning: invalid value encountered in multiply
[<matplotlib.lines.Line2D object at 0x10b1c7710>]
<matplotlib.text.Text object at 0x10b1c3d10>
<matplotlib.text.Text object at 0x10b1c9b90>

Now we simply minimize our Gwigglewiggle function. Based on the figure above, the miminum is near 0.45.

from scipy.optimize import fminbound

epsilonp_eq = fminbound(Gwigglewiggle, 0.4, 0.5)
print epsilonp_eq

plt.plot([epsilonp_eq], [Gwigglewiggle(epsilonp_eq)], 'ro')
plt.savefig('images/gibbs-minim-2.png')
>>> >>> 0.46959618249
>>> [<matplotlib.lines.Line2D object at 0x10d4d3e50>]

To compute equilibrium mole fractions we do this:

yi = (yi0 + nu_j[0]*epsilonp_eq) / (1.0 + epsilonp_eq*np.sum(nu_j))
yb = (yb0 + nu_j[1]*epsilonp_eq) / (1.0 + epsilonp_eq*np.sum(nu_j))
yp = (yp0 + nu_j[2]*epsilonp_eq) / (1.0 + epsilonp_eq*np.sum(nu_j))

print yi, yb, yp

# or this
y_j = (yj0 + np.dot(nu_j, epsilonp_eq)) / (1.0 + epsilonp_eq*np.sum(nu_j))
print y_j
>>> >>> >>> 0.0573220186324 0.0573220186324 0.885355962735
>>> ... >>> [ 0.05732202  0.05732202  0.88535596]

\(K = \frac{a_P}{a_I a_B} = \frac{y_p P/P^\circ}{y_i P/P^\circ y_b P/P^\circ} = \frac{y_P}{y_i y_b}\frac{P^\circ}{P}\).

We can express the equilibrium constant like this :\(K = \prod\limits_j a_j^{\nu_j}\), and compute it with a single line of code.

K = np.exp(-Grxn/R/T)
print 'K from delta G ',K
print 'K as ratio of mole fractions ',yp / (yi * yb) * P0 / P
print 'compact notation: ',np.prod((y_j * P / P0)**nu_j)
K from delta G  107.776294742
K as ratio of mole fractions  107.779200065
compact notation:  107.779200065

These results are very close, and only disagree because of the default tolerance used in identifying the minimum of our function. You could tighten the tolerances by setting options to the fminbnd function.

1 Summary

In this post we derived an equation for the Gibbs free energy of a reacting mixture and used it to find the equilibrium composition. In future posts we will examine some alternate forms of the equations that may be more useful in some circumstances.

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

org-mode source

Org-mode version = 8.2.7c

Read and Post Comments

Find the minimum distance from a point to a curve.

| categories: optimization | tags: | View Comments

A problem that can be cast as a constrained minimization problem is to find the minimum distance from a point to a curve. Suppose we have \(f(x) = x^2\), and the point (0.5, 2). what is the minimum distance from that point to \(f(x)\)?

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import fmin_cobyla

P = (0.5, 2)

def f(x):
    return x**2

def objective(X):
    x,y = X
    return np.sqrt((x - P[0])**2 + (y - P[1])**2)

def c1(X):
    x,y = X
    return f(x) - y

X = fmin_cobyla(objective, x0=[0.5,0.5], cons=[c1])

print 'The minimum distance is {0:1.2f}'.format(objective(X))

# Verify the vector to this point is normal to the tangent of the curve
# position vector from curve to point
v1 = np.array(P) - np.array(X)
# position vector
v2 = np.array([1, 2.0 * X[0]])
print 'dot(v1, v2) = ',np.dot(v1, v2)

x = np.linspace(-2, 2, 100)

plt.plot(x, f(x), 'r-', label='f(x)')
plt.plot(P[0], P[1], 'bo', label='point')
plt.plot([P[0], X[0]], [P[1], X[1]], 'b-', label='shortest distance')
plt.plot([X[0], X[0] + 1], [X[1], X[1] + 2.0 * X[0]], 'g-', label='tangent')
plt.axis('equal')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc='best')
plt.savefig('images/min-dist-p-func.png')
The minimum distance is 0.86
dot(v1, v2) =  0.000336477214214

   Normal return from subroutine COBYLA

   NFVALS =   44   F = 8.579598E-01    MAXCV = 0.000000E+00
   X = 1.300793E+00   1.692061E+00

In the code above, we demonstrate that the point we find on the curve that minimizes the distance satisfies the property that a vector from that point to our other point is normal to the tangent of the curve at that point. This is shown by the fact that the dot product of the two vectors is very close to zero. It is not zero because of the accuracy criteria that is used to stop the minimization is not high enough.

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

org-mode source

Read and Post Comments

Using constrained optimization to find the amount of each phase present

| categories: optimization | tags: | View Comments

The problem we solve here is that we have several compounds containing Ni and Al, and a bulk mixture of a particular composition of Ni and Al. We want to know which mixture of phases will minimize the total energy. The tricky part is that the optimization is constrained because the mixture of phases must have the overall stoichiometry we want. We formulate the problem like this.

Basically, we want to minimize the function \(E = \sum w_i E_i\), where \(w_i\) is the mass of phase \(i\), and \(E_i\) is the energy per unit mass of phase \(i\). There are some constraints to ensure conservation of mass. Let us consider the following compounds: Al, NiAl, Ni3Al, and Ni, and consider a case where the bulk composition of our alloy is 93.8% Ni and balance Al. We want to know which phases are present, and in what proportions. There are some subtleties in considering the formula and molecular weight of an alloy. We consider the formula with each species amount normalized so the fractions all add up to one. For example, Ni_3Al is represented as Ni_{0.75}Al_{0.25}, and the molecular weight is computed as 0.75*MW_{Ni} + 0.25*MW_{Al}.

We use scipy.optimize.fmin_slsqp to solve this problem, and define two equality constraint functions, and the bounds on each weight fraction.

Note: the energies in this example were computed by density functional theory at 0K.

import numpy as np
from scipy.optimize import fmin_slsqp

# these are atomic masses of each species
Ni = 58.693
Al = 26.982

COMPOSITIONS = ['Al', 'NiAl',              'Ni3Al',  'Ni']
MW = np.array(  [Al,  (Ni + Al)/2.0, (3*Ni + Al)/4.0, Ni])

xNi = np.array([0.0, 0.5, 0.75, 1.0])  # mole fraction of nickel in each compd
WNi = xNi*Ni / MW                      # weight fraction of Ni in each cmpd

ENERGIES = np.array([0.0, -0.7, -0.5, 0.0])

BNi = 0.938

def G(w):
    'function to minimize. w is a vector of weight fractions, ENERGIES is defined above.'
    return np.dot(w, ENERGIES)

def ec1(w):
    'conservation of Ni constraint'
    return BNi - np.dot(w, WNi)

def ec2(w):
    'weight fractions sum to one constraint'
    return 1 - np.sum(w)

w0 = np.array([0.0, 0.0, 0.5, 0.5]) # guess weight fractions

y = fmin_slsqp(G,   
               w0,
               eqcons=[ec1, ec2], 
               bounds=[(0,1)]*len(w0))

for ci, wi in zip(COMPOSITIONS, y):
    print '{0:8s} {1:+8.2%}'.format(ci, wi)
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.233299644373
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
Al         -0.00%
NiAl       +0.00%
Ni3Al     +46.66%
Ni        +53.34%

So, the sample will be about 47% by weight of Ni3Al, and 53% by weight of pure Ni.

It may be convenient to formulate this in terms of moles.

import numpy as np
from scipy.optimize import fmin_slsqp

COMPOSITIONS = ['Al', 'NiAl', 'Ni3Al',  'Ni']
xNi = np.array([0.0, 0.5, 0.75, 1.0])   # define this in mole fractions

ENERGIES = np.array([0.0, -0.7, -0.5, 0.0]) 

xNiB = 0.875  # bulk Ni composition

def G(n):
    'function to minimize'
    return np.dot(n, ENERGIES)

def ec1(n):
    'conservation of Ni'
    Ntot = np.sum(n)
    return (Ntot * xNiB) - np.dot(n,  xNi)

def ec2(n):
    'mole fractions sum to one'
    return 1 - np.sum(n)

n0 = np.array([0.0, 0.0, 0.45, 0.55]) # initial guess of mole fractions

y = fmin_slsqp(G,   
               n0,
               eqcons=[ec1, ec2], 
               bounds=[(0, 1)]*(len(n0)))

for ci, xi in zip(COMPOSITIONS, y):
    print '{0:8s} {1:+8.2%}'.format(ci, xi)
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -0.25
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
Al         +0.00%
NiAl       -0.00%
Ni3Al     +50.00%
Ni        +50.00%

This means we have a 1:1 molar ratio of Ni and Ni_{0.75}Al_{0.25}. That works out to the overall bulk composition in this particular problem.

Let us verify that these two approaches really lead to the same conclusions. On a weight basis we estimate 53.3%wt Ni and 46.7%wt Ni3Al, whereas we predict an equimolar mixture of the two phases. Below we compute the mole fraction of Ni in each case.

# these are atomic masses of each species
Ni = 58.693
Al = 26.982

# Molar case
# 1 mol Ni + 1 mol Ni_{0.75}Al_{0.25}
N1 = 1.0; N2 = 1.0
mol_Ni = 1.0 * N1 + 0.75 * N2
xNi = mol_Ni / (N1 + N2)
print xNi

# Mass case
M1 = 0.533; M2 = 0.467
MW1 = Ni; MW2 = 0.75*Ni + 0.25*Al

xNi2 = (1.0 * M1/MW1 + 0.75 * M2 / MW2) / (M1/MW1 + M2/MW2)
print xNi2
0.875
0.874192746385

You can see the overall mole fraction of Ni is practically the same in each case.

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

org-mode source

Read and Post Comments

Constrained minimization to find equilibrium compositions

| categories: optimization | tags: reaction engineering, thermodynamics | View Comments

adapated from Chemical Reactor analysis and design fundamentals, Rawlings and Ekerdt, appendix A.2.3.

Matlab post

The equilibrium composition of a reaction is the one that minimizes the total Gibbs free energy. The Gibbs free energy of a reacting ideal gas mixture depends on the mole fractions of each species, which are determined by the initial mole fractions of each species, the extent of reactions that convert each species, and the equilibrium constants.

Reaction 1: \(I + B \rightleftharpoons P1\)

Reaction 2: \(I + B \rightleftharpoons P2\)

Here we define the Gibbs free energy of the mixture as a function of the reaction extents.

import numpy as np

def gibbs(E):
    'function defining Gibbs free energy as a function of reaction extents'
    e1 = E[0]
    e2 = E[1]
    # known equilibrium constants and initial amounts
    K1 = 108; K2 = 284; P = 2.5
    yI0 = 0.5; yB0 = 0.5; yP10 = 0.0; yP20 = 0.0
    # compute mole fractions
    d = 1 - e1 - e2
    yI = (yI0 - e1 - e2) / d
    yB = (yB0 - e1 - e2) / d
    yP1 = (yP10 + e1) / d
    yP2 = (yP20 + e2) / d
    G = (-(e1 * np.log(K1) + e2 * np.log(K2)) +
         d * np.log(P) + yI * d * np.log(yI) +
         yB * d * np.log(yB) + yP1 * d * np.log(yP1) + yP2 * d * np.log(yP2))
    return G

The equilibrium constants for these reactions are known, and we seek to find the equilibrium reaction extents so we can determine equilibrium compositions. The equilibrium reaction extents are those that minimize the Gibbs free energy. We have the following constraints, written in standard less than or equal to form:

\(-\epsilon_1 \le 0\)

\(-\epsilon_2 \le 0\)

\(\epsilon_1 + \epsilon_2 \le 0.5\)

In Matlab we express this in matrix form as Ax=b where

\begin{equation} A = \left[ \begin{array}{cc} -1 & 0 \\ 0 & -1 \\ 1 & 1 \end{array} \right] \end{equation}

and

\begin{equation} b = \left[ \begin{array}{c} 0 \\ 0 \\ 0.5\end{array} \right] \end{equation}

Unlike in Matlab, in python we construct the inequality constraints as functions that are greater than or equal to zero when the constraint is met.

def constraint1(E):
    e1 = E[0]
    return e1


def constraint2(E):
    e2 = E[1]
    return e2


def constraint3(E):
    e1 = E[0]
    e2 = E[1]
    return 0.5 - (e1 + e2)

Now, we minimize.

from scipy.optimize import fmin_slsqp

X0 = [0.2, 0.2]
X = fmin_slsqp(gibbs, X0, ieqcons=[constraint1, constraint2, constraint3],
               bounds=((0.001, 0.499),
                       (0.001, 0.499)))
print(X)

print(gibbs(X))
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -2.55942394906
            Iterations: 7
            Function evaluations: 31
            Gradient evaluations: 7
[ 0.13336503  0.35066486]
-2.55942394906

One way we can verify our solution is to plot the gibbs function and see where the minimum is, and whether there is more than one minimum. We start by making grids over the range of 0 to 0.5. Note we actually start slightly above zero because at zero there are some numerical imaginary elements of the gibbs function or it is numerically not defined since there are logs of zero there. We also set all elements where the sum of the two extents is greater than 0.5 to near zero, since those regions violate the constraints.

import numpy as np
import matplotlib.pyplot as plt

def gibbs(E):
    'function defining Gibbs free energy as a function of reaction extents'
    e1 = E[0]
    e2 = E[1]
    # known equilibrium constants and initial amounts
    K1 = 108; K2 = 284; P = 2.5;
    yI0 = 0.5; yB0 = 0.5; yP10 = 0.0; yP20 = 0.0;
    # compute mole fractions
    d = 1 - e1 - e2;
    yI = (yI0 - e1 - e2)/d;
    yB = (yB0 - e1 - e2)/d;
    yP1 = (yP10 + e1)/d;
    yP2 = (yP20 + e2)/d;
    G = (-(e1 * np.log(K1) + e2 * np.log(K2)) +
         d * np.log(P) + yI * d * np.log(yI) +
         yB * d * np.log(yB) + yP1 * d * np.log(yP1) + yP2 * d * np.log(yP2))
    return G


a = np.linspace(0.001, 0.5, 100)
E1, E2 = np.meshgrid(a,a)

sumE = E1 + E2
E1[sumE >= 0.5] = 0.00001
E2[sumE >= 0.5] = 0.00001

# now evaluate gibbs
G = np.zeros(E1.shape)
m,n = E1.shape

G = gibbs([E1, E2])

CS = plt.contour(E1, E2, G, levels=np.linspace(G.min(),G.max(),100))
plt.xlabel('$\epsilon_1$')
plt.ylabel('$\epsilon_2$')
plt.colorbar()

plt.plot([0.13336503],  [0.35066486], 'ro')

plt.savefig('images/gibbs-minimization-1.png')
plt.savefig('images/gibbs-minimization-1.svg')
plt.show()

You can see we found the minimum. We can compute the mole fractions pretty easily.

e1 = X[0];
e2 = X[1];

yI0 = 0.5; yB0 = 0.5; yP10 = 0; yP20 = 0; #initial mole fractions

d = 1 - e1 - e2;
yI = (yI0 - e1 - e2) / d
yB = (yB0 - e1 - e2) / d
yP1 = (yP10 + e1) / d
yP2 = (yP20 + e2) / d

print('y_I = {0:1.3f} y_B = {1:1.3f} y_P1 = {2:1.3f} y_P2 = {3:1.3f}'.format(yI,yB,yP1,yP2))
y_I = 0.031 y_B = 0.031 y_P1 = 0.258 y_P2 = 0.680

1 summary

I found setting up the constraints in this example to be more confusing than the Matlab syntax.

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

org-mode source

Org-mode version = 8.2.10

Read and Post Comments

Using Lagrange multipliers in optimization

| categories: optimization | tags: | View Comments

Matlab post (adapted from http://en.wikipedia.org/wiki/Lagrange_multipliers.)

Suppose we seek to maximize the function \(f(x,y)=x+y\) subject to the constraint that \(x^2 + y^2 = 1\). The function we seek to maximize is an unbounded plane, while the constraint is a unit circle. We want the maximum value of the circle, on the plane. We plot these two functions here.

import numpy as np

x = np.linspace(-1.5, 1.5)

[X, Y] = np.meshgrid(x, x)

import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.gca(projection='3d')

ax.plot_surface(X, Y, X + Y)

theta = np.linspace(0,2*np.pi);
R = 1.0
x1 = R * np.cos(theta)
y1 = R * np.sin(theta)

ax.plot(x1, y1, x1 + y1, 'r-')
plt.savefig('images/lagrange-1.png')

1 Construct the Lagrange multiplier augmented function

To find the maximum, we construct the following function: \(\Lambda(x,y; \lambda) = f(x,y)+\lambda g(x,y)\) where \(g(x,y) = x^2 + y^2 - 1 = 0\), which is the constraint function. Since \(g(x,y)=0\), we are not really changing the original function, provided that the constraint is met!

import numpy as np

def func(X):
    x = X[0]
    y = X[1]
    L = X[2] # this is the multiplier. lambda is a reserved keyword in python
    return x + y + L * (x**2 + y**2 - 1)

2 Finding the partial derivatives

The minima/maxima of the augmented function are located where all of the partial derivatives of the augmented function are equal to zero, i.e. \(\partial \Lambda/\partial x = 0\), \(\partial \Lambda/\partial y = 0\), and \(\partial \Lambda/\partial \lambda = 0\). the process for solving this is usually to analytically evaluate the partial derivatives, and then solve the unconstrained resulting equations, which may be nonlinear.

Rather than perform the analytical differentiation, here we develop a way to numerically approximate the partial derivatives.

def dfunc(X):
    dLambda = np.zeros(len(X))
    h = 1e-3 # this is the step size used in the finite difference.
    for i in range(len(X)):
        dX = np.zeros(len(X))
        dX[i] = h
        dLambda[i] = (func(X+dX)-func(X-dX))/(2*h);
    return dLambda

3 Now we solve for the zeros in the partial derivatives

The function we defined above (dfunc) will equal zero at a maximum or minimum. It turns out there are two solutions to this problem, but only one of them is the maximum value. Which solution you get depends on the initial guess provided to the solver. Here we have to use some judgement to identify the maximum.

from scipy.optimize import fsolve

# this is the max
X1 = fsolve(dfunc, [1, 1, 0])
print X1, func(X1)

# this is the min
X2 = fsolve(dfunc, [-1, -1, 0])
print X2, func(X2)
>>> ... >>> [ 0.70710678  0.70710678 -0.70710678] 1.41421356237
>>> ... >>> [-0.70710678 -0.70710678  0.70710678] -1.41421356237

4 Summary

Three dimensional plots in matplotlib are a little more difficult than in Matlab (where the code is almost the same as 2D plots, just different commands, e.g. plot vs plot3). In Matplotlib you have to import additional modules in the right order, and use the object oriented approach to plotting as shown here.

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

org-mode source

Read and Post Comments

« Previous Page -- Next Page »