Lagrange multipliers example

This is an experiment with org-mode and ob-python that simulates a notebok environment which mix code, text and math (latex). I am going through example presented in John Kitchin work, which is extremely useful for someone planning on start coding in python for science.

I take no credit from the material presented here, it is a copy with a new format. I am particularly interested in this subject due an ongoing project in structural optimization therefore I will put an effort on detailing some of the method derivation.

We seek to maximize the function

$$f(x,y)=x+y$$

subjected to a circle constraint

$$x^2+y^2=1$$

which can be visualized with

copy 1import numpy as np
 2
 3x = np.linspace(-1.5, 1.5)
 4
 5[X, Y] = np.meshgrid(x, x)
 6
 7%matplotlib inline
 8import matplotlib as mpl
 9from mpl_toolkits.mplot3d import Axes3D
10import matplotlib.pyplot as plt
11
12fig = plt.figure()
13ax = fig.gca(projection='3d')
14
15ax.plot_surface(X, Y, X+Y)
16
17theta = np.linspace(0, 2*np.pi)
18R = 1.0
19x1 = np.cos(theta)
20y1 = np.sin(theta)
21
22ax.plot(x1, y1, x1+y1, 'r-')
23plt.tight_layout()

In order to find the maximum we construct the function

$$\mathrm{\Lambda }(x,y;\lambda )=f(x,y)+\lambda g(x,y)$$

where, the constraint function is

$$g(x,y)=x^2+y^2-1=0$$

in code:

copy1def func(X):
2    x = X[0]
3    y = X[1]
4    L = X[2] # the multiplier
5    return x + y + L *(x**2 + y**2 -1)

The maximum and minimum of the augmented function are located where all of the partial derivatives are equal to zero,

$${\displaystyle \frac{\partial \mathrm{\Lambda }}{\partial x}}=0\text{and}{\displaystyle \frac{\partial \mathrm{\Lambda }}{\partial y}}=0.$$

The partial derivatives are usually obtained analytically, which can be hard depending on the function. It is more convenient to compute them numerically, this is done using finite differences

$${\displaystyle \frac{\partial f}{\partial x}}={\displaystyle \frac{f(x+\mathrm{\Delta }x)-f(x-\mathrm{\Delta }x)}{2\mathrm{\Delta }x}}$$

copy1def dfunc(X):
2    dLambda = np.zeros(len(X))
3    h = 1e-3 # step size used in the finite differences
4    for i in range(len(X)):
5        dX = np.zeros(len(X))
6        dX[i] = h
7        dLambda[i] = (func(X + dX) - func(X - dX))/(2*h)
8    return dLambda

Setting the derivative dfunc to zero will enable was to find the maximum or minimum.

copy1from scipy.optimize import fsolve
2
3# the max
4X1 = fsolve(dfunc, [1, 1, 0])
5print(X1, func(X1))
6
7# the min
8X2 = fsolve(dfunc, [-1, -1, 0])
9print(X2, func(X2))

[ 0.70710678  0.70710678 -0.70710678] 1.41421356237
[-0.70710678 -0.70710678  0.70710678] -1.41421356237

The points found before $X_1$ and $X_2$ are the maximum and minimum, respectively, of the function $f(x,y)=x+y$ subjected to the restricted domain defined by $g(x,y)=x^2+y^2-1$.

copy1fig2 = plt.figure() # create an empty figure with no axes
2ax2 = Axes3D(fig2) # create a 3D axes object

copy1fig3 = plt.figure()
2ax3 = Axes3D(fig3)
3ax3.plot_surface(X, Y, X+Y, color='y', rstride=1, cstride=1, alpha=0.5, edgecolor='w')

copy1fig4 = plt.figure()
2ax4 = Axes3D(fig4)
3ax4.plot_surface(X, Y, X+Y, color='y', rstride=1, cstride=1, alpha=0.5, edgecolor='w')
4
5theta = np.linspace(0, 2*np.pi)
6R = 1.0
7x1 = R * np.cos(theta)
8y1 = R * np.sin(theta)
9ax4.plot(x1, y1, x1+y1)

copy 1fig5 = plt.figure()
 2ax5 = Axes3D(fig5)
 3ax5.plot_surface(X, Y, X+Y, color='y', rstride=100, cstride=100, alpha=0.2, edgecolor='w')
 4ax5.plot_surface(X, Y, 0, color='b', alpha=0.05, rstride=100, cstride=100)
 5
 6theta = np.linspace(0, 2*np.pi)
 7R = 1.0
 8x1 = R * np.cos(theta)
 9y1 = R * np.sin(theta)
10ax5.plot(x1, y1, alpha=0.3)
11ax5.plot(x1, y1, x1+y1, c='y', alpha=0.8)
12
13ax5.scatter(X1[0], X1[1], X1[2], s=30, c='r')
14ax5.scatter(X2[0], X2[1], X2[2], s=30, c='g')

We can see from the last plot that the maximum and minimum point are on the edge of the restrained domain, which makes sense, since the function minimized/maximized is linear on both variables.