Trusted, Bureaucrats, Administrators
1,660
edits
Changes
no edit summary
# http://docs.sympy.org/dev/gotchas.html
# https://github.com/sympy/sympy/wiki/Faq
== Sumpy python installation ==Consider for example installing [[Epd_canopy | Enthought Canopy]] == Tutorial ==Created by: Troels Emtekær Linnet <br>PhD student <br>Copenhagen University <br>SBiNLab run by: '''python sympy_test.py''' <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**2print("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 chi2.")# Function to calculate the chi2.# 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_decchi2 = (1.0 / errors * (values - back_calc))**2print("The chi2 function is:", chi2) 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 chi2 function for exponential decay, with respect to r2eff")d_chi2_d_r2eff = diff(chi2, r2eff)print("d_chi2_d_r2eff=", d_chi2_d_r2eff)print("Should be: (2 * i0 * times * (values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * self.errors**2) ")val_d_chi2_d_r2eff = ( 2 * i0 * times * (values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2)print("Are they equal?: ", simplify(d_chi2_d_r2eff) == simplify(val_d_chi2_d_r2eff) ) print("") print("Derivative of chi2 function for exponential decay, with respect to i0")d_chi2_d_i0 = diff(chi2, i0)print("d_chi2_d_i0=", d_chi2_d_i0)print("Should be: - ( 2. * ( values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2) ")val_d_chi2_d_i0 = - ( 2. * ( values - i0 / exp(r2eff * times) ) ) / ( exp(r2eff * times) * errors**2) print("Are they equal?: ", simplify(d_chi2_d_i0) == simplify(val_d_chi2_d_i0) ) print("\n") print("""Form the Jacobian matrix by:------------------------------------------------------------------------------from numpy import array, transpose d_chi2_d_r2eff = %sd_chi2_d_i0 = %sjacobian_matrix = transpose(array( [df_d_r2eff , df_d_i0] ) )------------------------------------------------------------------------------""" % (d_chi2_d_r2eff, d_chi2_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_chi2_d_r2eff_d_r2eff = diff(chi2, r2eff, r2eff)d2_chi2_d_r2eff_d_i0 = diff(chi2, r2eff, i0)d2_chi2_d_i0_d_r2eff = diff(chi2, i0, r2eff)d2_chi2_d_i0_d_i0 = diff(chi2, i0, i0)print("d2_chi2_d_r2eff_d_r2eff=", d2_chi2_d_r2eff_d_r2eff)print("d2_chi2_d_r2eff_d_i0=", d2_chi2_d_r2eff_d_i0)print("d2_chi2_d_i0_d_r2eff=", d2_chi2_d_i0_d_r2eff)print("d2_chi2_d_i0_d_i0=", d2_chi2_d_i0_d_i0) print("\n") print("""Form the Hessian matrix by:------------------------------------------------------------------------------from numpy import array, transpose d2_chi2_d_r2eff_d_r2eff = %sd2_chi2_d_r2eff_d_i0 = %sd2_chi2_d_i0_d_r2eff = %sd2_chi2_d_i0_d_i0 = %shessian_matrix = transpose(array( [d2_chi2_d_r2eff_d_r2eff, d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0] ) )------------------------------------------------------------------------------""" % (d2_chi2_d_r2eff_d_r2eff,d2_chi2_d_r2eff_d_i0, d2_chi2_d_i0_d_r2eff, d2_chi2_d_i0_d_i0) )</source>