Difference between revisions of "DPL94 derivatives"
Jump to navigation
Jump to search
(→sympy) |
|||
Line 47: | Line 47: | ||
------------------------------------------------------------------------------ | ------------------------------------------------------------------------------ | ||
""" % (d_f_d_R1, d_f_d_theta, d_f_d_R1rho_p, d_f_d_phi_ex, d_f_d_kex, d_f_d_we) ) | """ % (d_f_d_R1, d_f_d_theta, d_f_d_R1rho_p, d_f_d_phi_ex, d_f_d_kex, d_f_d_we) ) | ||
+ | |||
+ | #### Method 2 | ||
+ | # http://docs.sympy.org/0.7.2/modules/matrices/matrices.html | ||
+ | |||
+ | # The vectorial function. | ||
+ | X = Matrix([R1 * cos(theta)**2 + (R1rho_p + ( (phi_ex * kex) / (kex**2 + we**2) ) ) * sin(theta)**2]) | ||
+ | # What to derive for. | ||
+ | Y = Matrix([R1, theta, R1rho_p, phi_ex, kex, we]) | ||
+ | |||
+ | # Make the Jacobian | ||
+ | Jacobian = X.jacobian(Y) | ||
+ | |||
+ | jac_string = str(Jacobian) | ||
+ | jac_string_arr = jac_string.replace("Matrix", "array") | ||
+ | |||
+ | print("""Form the Jacobian matrix by: | ||
+ | ------------------------------------------------------------------------------ | ||
+ | from numpy import array, cos, sin, pi, transpose | ||
+ | |||
+ | R1 = 1.1 | ||
+ | theta = pi / 4 | ||
+ | R1rho_p = 10. | ||
+ | phi_ex = 1100. | ||
+ | kex = 2200. | ||
+ | we = 3300. | ||
+ | |||
+ | jacobian_matrix_2 = %s | ||
+ | |||
+ | print jacobian_matrix_2 | ||
+ | ------------------------------------------------------------------------------ | ||
+ | """ % (jac_string_arr) ) | ||
</source> | </source> | ||
+ | |||
=== output === | === output === | ||
output | output |
Revision as of 09:29, 1 September 2014
Jacobian
sympy
from sympy import *
# In contrast to other Computer Algebra Systems, in SymPy you have to declare symbolic variables explicitly:
R1 = Symbol('R1')
theta = Symbol('theta')
R1rho_p = Symbol('R1rho_p')
phi_ex = Symbol('phi_ex')
kex = Symbol('kex')
we = Symbol('we')
# Define function
f = R1 * cos(theta)**2 + (R1rho_p + ( (phi_ex * kex) / (kex**2 + we**2) ) ) * sin(theta)**2
print("Now calculate the Jacobian. The partial derivative matrix.\n")
print("Jacobian is m rows with function derivatives and n columns of parameters.")
d_f_d_R1 = diff(f, R1)
d_f_d_theta = diff(f, theta)
d_f_d_R1rho_p = diff(f, R1rho_p)
d_f_d_phi_ex = diff(f, phi_ex)
d_f_d_kex = diff(f, kex)
d_f_d_we = diff(f, we)
print("""Form the Jacobian matrix by:
------------------------------------------------------------------------------
from numpy import array, cos, sin, pi, transpose
R1 = 1.1
theta = pi / 4
R1rho_p = 10.
phi_ex = 1100.
kex = 2200.
we = 3300.
d_f_d_R1 = %s
d_f_d_theta = %s
d_f_d_R1rho_p = %s
d_f_d_phi_ex = %s
d_f_d_kex = %s
d_f_d_we = %s
jacobian_matrix = transpose(array( [d_f_d_R1 , d_f_d_theta, d_f_d_R1rho_p, d_f_d_phi_ex, d_f_d_kex, d_f_d_we] ) )
print jacobian_matrix
------------------------------------------------------------------------------
""" % (d_f_d_R1, d_f_d_theta, d_f_d_R1rho_p, d_f_d_phi_ex, d_f_d_kex, d_f_d_we) )
#### Method 2
# http://docs.sympy.org/0.7.2/modules/matrices/matrices.html
# The vectorial function.
X = Matrix([R1 * cos(theta)**2 + (R1rho_p + ( (phi_ex * kex) / (kex**2 + we**2) ) ) * sin(theta)**2])
# What to derive for.
Y = Matrix([R1, theta, R1rho_p, phi_ex, kex, we])
# Make the Jacobian
Jacobian = X.jacobian(Y)
jac_string = str(Jacobian)
jac_string_arr = jac_string.replace("Matrix", "array")
print("""Form the Jacobian matrix by:
------------------------------------------------------------------------------
from numpy import array, cos, sin, pi, transpose
R1 = 1.1
theta = pi / 4
R1rho_p = 10.
phi_ex = 1100.
kex = 2200.
we = 3300.
jacobian_matrix_2 = %s
print jacobian_matrix_2
------------------------------------------------------------------------------
""" % (jac_string_arr) )
output
output
from numpy import array, cos, sin, pi, transpose
R1 = 1.1
theta = pi / 4
R1rho_p = 10.
phi_ex = 1100.
kex = 2200.
we = 3300.
d_f_d_R1 = cos(theta)**2
d_f_d_theta = -2*R1*sin(theta)*cos(theta) + 2*(R1rho_p + kex*phi_ex/(kex**2 + we**2))*sin(theta)*cos(theta)
d_f_d_R1rho_p = sin(theta)**2
d_f_d_phi_ex = kex*sin(theta)**2/(kex**2 + we**2)
d_f_d_kex = (-2*kex**2*phi_ex/(kex**2 + we**2)**2 + phi_ex/(kex**2 + we**2))*sin(theta)**2
d_f_d_we = -2*kex*phi_ex*we*sin(theta)**2/(kex**2 + we**2)**2
jacobian_matrix = transpose(array( [d_f_d_R1 , d_f_d_theta, d_f_d_R1rho_p, d_f_d_phi_ex, d_f_d_kex, d_f_d_we] ) )
print jacobian_matrix