I am trying to make a Jacobian matrix using sympy.
I want to make a matrix of 2nd order derivatives, which I will call a 2nd order Jacobian.
I did this:
gk_letters = ["Δ", "Γ", "v"]
def greek_symbols():
return symbols(" ".join(gk_letters))
def partials_1st_order():
rows = []
gks = greek_symbols()
for idx1 in range(len(gk_letters)):
gk1 = gks[idx1]
col = []
for idx2 in range(len(gk_letters)):
gk2 = gks[idx2]
if gk1 == gk2:
col.append(1)
else:
col.append(Derivative(gk1, gk2))
rows.append(col)
return gks, Matrix(rows)
def partials_2nd_order():
gks, jacobian = partials_1st_order()
y = zeros(len(gks)**2, len(gks))
input_rows, input_cols = jacobian.shape
output_chunk = 0
for den_idx in range(len(gks)):
denominator = gks[den_idx]
for input_row_num in range(input_rows):
for input_col_num in range(input_cols):
numerator = jacobian[input_row_num, input_col_num]
new_diff = Derivative(numerator, denominator)
y[output_chunk*len(gks) + input_row_num, input_col_num] = new_diff
output_chunk += 1
return y
When I call partials_2nd_order() at the command prompt, I get this:
⎡ 2 2 ⎤
⎢ d d d ⎥
⎢ ──(1) ─────(Δ) ─────(Δ)⎥
⎢ dΔ dΔ dΓ dΔ dv ⎥
⎢ ⎥
⎢ 2 2 ⎥
⎢ d d d ⎥
⎢ ───(Γ) ──(1) ─────(Γ)⎥
⎢ 2 dΔ dΔ dv ⎥
⎢ dΔ ⎥
⎢ ⎥
...
⎢ ⎥
⎢ 2 2 ⎥
⎢ d d d ⎥
⎢─────(v) ─────(v) ──(1) ⎥
⎣dv dΔ dv dΓ dv ⎦
I want to collapse to 0 everything that will always be 0, such as derivatives of constants.
However, when I try this:
[jj.is_constant() for jj in x]
I get this:
[True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, True]
In other words, sympy thinks everything in this matrix is constant.
What am I missing?
All the derivatives in your matrix are zero. The partial derivative of x with respect to y is zero:
In [1]: x, y = symbols('x, y')
In [2]: d = Derivative(x, y)
In [3]: d
Out[3]:
d
──(x)
dy
In [4]: d.is_constant()
Out[4]: True
In [5]: d.doit()
Out[5]: 0
In [6]: diff(x, y)
Out[6]: 0
If you want a symbol that is a function of other symbols then you need a function:
In [10]: f = Function('f')
In [11]: d = Derivative(f(x), x)
In [12]: d
Out[12]:
d
──(f(x))
dx
In [13]: d.doit()
Out[13]:
d
──(f(x))
dx
In [14]: d.is_constant()
Out[14]: False
There is no way to make a set of 3 symbols that are all functions of each other though which seems to be what you want here.
You can keep the derivatives unevaluated though and just remove the ones that have a constant as the expression to be differentiated:
In [23]: d = [Derivative(x, y), Derivative(1, x)]
In [24]: d
Out[24]:
⎡d d ⎤
⎢──(x), ──(1)⎥
⎣dy dx ⎦
In [25]: [di for di in d if not di.expr.is_constant()]
Out[25]:
⎡d ⎤
⎢──(x)⎥
⎣dy ⎦