Trusted, Bureaucrats, Administrators
1,660
edits
Changes
<source lang="Python">
# Tutorial from: http://scipy-lectures.github.io/advanced/sympy.html
from sympy import *
# In contrast to other Computer Algebra Systems, in SymPy you have to declare symbolic variables explicitly:
x = Symbol('x')
y = Symbol('y')
# You can differentiate any SymPy expression using diff(func, var). Examples:
print("First some tests")
f = sin(x)
print("Function is", f)
print("df/dx is:", diff(f, x))
f = x**2
print("Function is", f)
print("df/dx is:", diff(f, x))
print("")
# Define exponential decay.
i0 = Symbol('i0')
times = Symbol('times')
r2eff = Symbol('r2eff')
print("Function to calculate intensity decay at time[i].")
exp_dec = i0 * exp( -times * r2eff)
print("Function f(r2eff, i0) is:", simplify(exp_dec))
print("df/d_i0 should be: exp(-r2eff*times)", simplify(diff(exp_dec, i0)) )
print("df/d_r2eff should be: -i0*times*exp(-r2eff*times)", simplify(diff(exp_dec, r2eff)) )
print("")
print("Function to calculate the weighted difference between function evaluation with fitted parameters and measured values.")
# Function to calculate the func.
# The standard deviation of the measured values.
errors = Symbol('errors')
# The measured values.
values = Symbol('values')
# back_calc, the function which return the values, according to parameters.
back_calc = exp_dec
func = 1.0 / errors * (values - back_calc)
print("The func function is:", func)
print("")
print("Now calculate the Jacobian. The partial derivative matrix.\n")
print("Jacobian is m rows with function derivatives and n columns of parameters.")
# If you want to test for symbolic equality, one way is to subtract one expression from the other and run it through functions like expand(), simplify(), and trigsimp() and see if the equation reduces to 0.
# https://github.com/sympy/sympy/wiki/Faq - Why does SymPy say that two equal expressions are unequal?
print("Derivative of func function for exponential decay, with respect to r2eff")
d_func_d_r2eff = diff(func, r2eff)
print("d_func_d_r2eff=", d_func_d_r2eff)
print("")
print("Derivative of func function for exponential decay, with respect to i0")
d_func_d_i0 = diff(func, i0)
print("d_func_d_i0=", d_func_d_i0)
print("\n")
print("""Form the Jacobian matrix by:
------------------------------------------------------------------------------
from numpy import array, transpose
d_func_d_r2eff = %s
d_func_d_i0 = %s
jacobian_matrix = transpose(array( [d_func_d_r2eff , d_func_d_i0] ) )
------------------------------------------------------------------------------
""" % (d_func_d_r2eff, d_func_d_i0) )
print("")
print("Now calculate the Hessian. The second-order partial derivatives.\n")
print("See for example: http://maxima-online.org/articles/hessian.html")
print("If all second partial derivatives of f exist, then the Hessian matrix of f is the matrix:")
print("Hf(x)i,j = d**2 f(x) / ( dx_i * dx_j )")
print("We need to do:")
print("diff(f(r2eff, i0), r2eff, r2eff)")
print("diff(f(r2eff, i0), r2eff, i0)")
print("diff(f(r2eff, i0), i0, r2eff)")
print("diff(f(r2eff, i0), i0, i0)")
d2_func_d_r2eff_d_r2eff = diff(func, r2eff, r2eff)
d2_func_d_r2eff_d_i0 = diff(func, r2eff, i0)
d2_func_d_i0_d_r2eff = diff(func, i0, r2eff)
d2_func_d_i0_d_i0 = diff(func, i0, i0)
print("d2_func_d_r2eff_d_r2eff=", d2_func_d_r2eff_d_r2eff)
print("d2_func_d_r2eff_d_i0=", d2_func_d_r2eff_d_i0)
print("d2_func_d_i0_d_r2eff=", d2_func_d_i0_d_r2eff)
print("d2_func_d_i0_d_i0=", d2_func_d_i0_d_i0)
print("\n")
print("""Form the Hessian matrix by:
------------------------------------------------------------------------------
from numpy import array, transpose
d2_func_d_r2eff_d_r2eff = %s
d2_func_d_r2eff_d_i0 = %s
d2_func_d_i0_d_r2eff = %s
d2_func_d_i0_d_i0 = %s
hessian_matrix = transpose(array( [d2_func_d_r2eff_d_r2eff, d2_func_d_r2eff_d_i0, d2_func_d_i0_d_r2eff, d2_func_d_i0_d_i0] ) )
------------------------------------------------------------------------------
""" % (d2_func_d_r2eff_d_r2eff,d2_func_d_r2eff_d_i0, d2_func_d_i0_d_r2eff, d2_func_d_i0_d_i0) )
from sympy import Function, hessian, pprint
f = Function('f')(r2eff, i0)
</source>