r25517 - /branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py -- September 01, 2014 - 21:45

Author: tlinnet
Date: Mon Sep  1 21:45:15 2014
New Revision: 25517

URL: http://svn.gna.org/viewcvs/relax?rev=25517&view=rev
Log:
For reference, added small script with symbolic derivation of TSMFK01.

task #7824(https://gna.org/task/index.php?7824): Model parameter ERROR 
estimation from Jacobian and Co-variance matrix of dispersion models.

Added:
    
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py

Added: 
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py
URL: 
http://svn.gna.org/viewcvs/relax/branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py?rev=25517&view=auto
==============================================================================
--- 
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py
    (added)
+++ 
branches/est_par_error/test_suite/shared_data/dispersion/estimate_par_err/tsmfk01/tsmfk01_der.py
    Mon Sep  1 21:45:15 2014
@@ -0,0 +1,54 @@
+from sympy import *
+
+# In contrast to other Computer Algebra Systems, in SymPy you have to 
declare symbolic variables explicitly:
+r20a = Symbol('r20a')
+dw = Symbol('dw')
+k_AB = Symbol('k_AB')
+tcp = Symbol('tcp')
+
+denom = dw * tcp
+numer = sin(denom)
+back_calc = r20a + k_AB - k_AB * numer / denom
+
+d_f_d_r20a = diff(back_calc, r20a)
+d_f_d_dw = diff(back_calc, dw)
+d_f_d_k_AB = diff(back_calc, k_AB)
+
+print("""Form the Jacobian matrix by:
+------------------------------------------------------------------------------
+
+d_f_d_r20a = %s
+d_f_d_dw = %s
+d_f_d_k_AB = %s
+jacobian_matrix = transpose(array( [d_f_d_r20a , d_f_d_dw, d_f_d_k_AB] ) )
+
+d_f_d_r20a = %s
+d_f_d_dw = %s
+d_f_d_k_AB = %s
+jacobian_matrix = transpose(array( [d_f_d_r20a , d_f_d_dw, d_f_d_k_AB] ) )
+
+------------------------------------------------------------------------------
+""" % (d_f_d_r20a, d_f_d_dw, d_f_d_k_AB, simplify(d_f_d_r20a), 
simplify(d_f_d_dw), simplify(d_f_d_k_AB)) )
+
+# Try again.
+
+# The vectorial function.
+X = Matrix([r20a + k_AB - k_AB * sin(dw * tcp) / dw * tcp])
+# What to derive for.
+Y = Matrix([r20a, dw, k_AB, tcp])
+
+# 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
+
+jacobian_matrix_2 = %s
+
+print jacobian_matrix_2
+------------------------------------------------------------------------------
+""" % (jac_string_arr) )