Figure 1.4:

Two quadratic functions and their sum.

Code for Figure 1.4

Text of the GNU GPL.

main.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
# Example sum of two quadratic functions.

import numpy as np
import matplotlib.pyplot as plt
from scipy.io import savemat
from misc import gnuplotsave
from misc import ellipse
import os

# def ellipse(A, level, npts, center):
#     """Generate points on an ellipse"""
#     theta = np.linspace(0, 2*np.pi, npts)
#     eigvals, eigvecs = np.linalg.eig(A)

#     # Get points on unit circle
#     circle_pts = np.vstack([np.cos(theta), np.sin(theta)])

#     # Transform circle to ellipse
#     scale = np.sqrt(level/eigvals)[:, np.newaxis]
#     transform = eigvecs @ np.diag(scale.flatten())
#     ellipse_pts = transform @ circle_pts

#     # Shift by center point
#     x = ellipse_pts[0,:] + center[0]
#     y = ellipse_pts[1,:] + center[1]

#     return x, y

# Example sum of two quadratic functions
a = np.array([-1, 0]).reshape(-1,1)
b = np.array([1, 1]).reshape(-1,1)
A = np.array([[1.25000, 0.75000], [0.75000, 1.25000]])
B = np.array([[1.50000, -0.50000], [-0.50000, 1.50000]])

# add the two quadratics
H = A + B
h = np.linalg.solve(H, A@a + B@b)
d = -(A@a + B@b).T@h + a.T@A@a + b.T@B@b

# plot resulting ellipse for combining two ellipses
# V(x) = V_1(x) + V_2(x)
level = 0.5
npts = 50
xA, yA, *_ = ellipse(A, level, npts, a.flatten())

xB, yB, *_ = ellipse(B, level, npts, b.flatten())

level = (2 - d/2)*2
xH, yH, *_ = ellipse(H, level, npts, h.flatten())

if not os.getenv('OMIT_PLOTS') == 'true':
    plt.plot(xA, yA, '-r', a[0], a[1], 'or',
             xB, yB, '-g', b[0], b[1], 'og',
             xH, yH, '-b', h[0], h[1], 'ob')
    plt.legend(['A', 'a', 'B', 'b', 'H = A + B', 'h'],
              loc='upper left')

# Save data
data = {}
data['contours'] = np.column_stack([xA, yA, xB, yB, xH, yH])
data['centers'] = np.column_stack([a.T, b.T, h.T])

# Save to dat file
gnuplotsave('nestedV.dat', data);