-
Notifications
You must be signed in to change notification settings - Fork 0
/
jacobian.py
79 lines (60 loc) · 1.91 KB
/
jacobian.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
65
66
67
68
69
70
71
72
73
74
75
76
77
78
import sympy as sp
from sympy import pprint
import numpy as np
# jacobian matricies for multivariate functions
# 2d example
def a1():
x1, x2 = sp.symbols('x1 x2')
f1 = 5*x1*x2
f2 = x1**2 * x2**2 + x1 + 2*x2
f = sp.Matrix([f1, f2])
X = sp.Matrix([x1, x2])
Df = f.jacobian(X)
Df0 = Df.subs([(x1, 1), (x2, 2)])
pprint(Df)
pprint(Df0)
# a1()
# 3d example
def a2():
x1, x2, x3 = sp.symbols('x1 x2 x3')
# assemble vectors f and x
f1 = sp.log(x1**2+x2**2)+x3**2
f2 = sp.exp(x2**2+x3**2)+x1**2
f3 = 1/(x3**2+x1**2)+x2**2
f = sp.Matrix([f1, f2, f3])
X = sp.Matrix([x1, x2, x3])
Df = f.jacobian(X)
# compute for actual values
Df0 = Df.subs([(x1, 1), (x2, 2), (x3, 3)])
pprint(Df)
pprint(Df0)
# a2()
# example linearization of a multivariable function via tangent plane
def a3():
x1, x2, x3 = sp.symbols('x1 x2 x3')
# define vectors
f1 = x1+x2**2-x3**2-13
f2 = sp.log(x2/4)+sp.exp(x3/2-1)
f3 = (x2-3)**2-x3**3+7
f = sp.Matrix([f1,f2,f3])
X = sp.Matrix([x1,x2,x3])
# Df(x) substitutes variables i.e 3*x1^2 inside partial derivative formuli with value
Df = f.jacobian(X)
# lambdify for use with numpy values
func = sp.lambdify([(x1,x2,x3)], f, "numpy")
jac = sp.lambdify([(x1,x2,x3)], Df, "numpy")
# evaluate f(v0)
v0 = np.array([1.5,3,2.5])
# evaluate jacobian at point v0 for slopes in all dimensions
DF = jac(v0).flatten()
# now choose v1 in vicinity of v0
offset = np.array([0.1,0.1,0.1])
v1 = v0+offset
# evaluating f(v0) and following the slope along the linearization
# should yeald a result close to actually evaluating f(v1) while still in v0's vicinity.
pprint(func(v1))
# tangent plane equation for linearization
# around neighbourhood of point v
def g(v_next,v): return (func(v).flatten()+jac(v)@(v_next-v)).flatten()
pprint(g(v1,v0))
a3()