00001 /*************************************************************************** 00002 jacobian.h - description 00003 ------------------- 00004 begin : Fri Apr 21 2000 00005 copyright : (C) 2000 by Michael Peeters 00006 email : Michael.Peeters@vub.ac.be 00007 ***************************************************************************/ 00008 00009 /*************************************************************************** 00010 * * 00011 * This program is free software; you can redistribute it and/or modify * 00012 * it under the terms of the GNU General Public License as published by * 00013 * the Free Software Foundation; either version 2 of the License, or * 00014 * (at your option) any later version. * 00015 * * 00016 ***************************************************************************/ 00017 00018 #ifndef JACOBIAN_H 00019 #define JACOBIAN_H 00020 00021 #include "numerictypes.h" 00022 #include "utility.h" 00023 #include <stdexcept> 00024 00025 namespace MODEL { 00026 00032 template <integer dims, class NT = NumericTraits<number,dims> > 00033 class Jacobian { 00034 public: 00035 typedef typename NT::number numT; 00036 typedef typename NT::vect vect; 00037 typedef typename NT::matrix matrix; 00038 typedef typename NT::vf vf; 00039 00042 Jacobian(vf& rvf, numT JEps=1E-4, bool own=false) : 00043 owned(own), 00044 f((own)?(rvf.clone()):(&rvf)), 00045 epsilon(JEps) {} 00046 ~Jacobian() {if (owned) delete f;} 00048 Jacobian(const Jacobian& cj) : owned(false), f(cj.f) {} 00049 00051 numT calculate_didj(integer i, integer j, const vect& u); 00052 00054 vect calculate_di(integer i, const vect& u); 00055 00058 matrix calculate(const vect& u, const vect& fu); 00059 00063 matrix calculate_accurate(const vect& u, const vect& fu); 00064 00067 matrix calculate(const vect& u) 00068 { vect fu( (*f)(u) ); return calculate(u,fu);} 00069 00072 matrix calculate_accurate(const vect& u) 00073 { vect fu( (*f)(u) ); return calculate_accurate(u,fu);} 00074 00076 vf& get_function(void) 00077 { 00078 return *f; 00079 } 00080 00081 private: 00083 PRIVATE_ASSIGN(Jacobian); 00084 00085 private: 00086 bool owned; 00087 vf* f; 00089 number epsilon; 00090 }; 00091 00092 // This returns: 00093 // # Jacobian 00094 // 0 1 0 = df[0]/dx[j] 00095 // 0 0 1 = df[1]/dx[j] 00096 // 1 0 0 = df[2]/dx[j] 00097 // so the matrix is organised along rows: matrix[i] gives all 00098 // the derivatives of one component. 00099 // j[i][j]=df[i]dx[j] 00100 00101 template <integer dims, class NT> 00102 typename NT::matrix 00103 Jacobian<dims,NT>::calculate(const vect& u, const vect& fu) 00104 { 00105 // There might be a faster way to do this, but I haven found it yet 00106 matrix jac(0.); 00107 00108 for(integer j=0;j<dims;j++) 00109 { 00110 vect udu(u); 00111 00112 numT du = epsilon*abs(u[j]); 00113 00114 if(du==0.0) du = epsilon; // avoid numerical error 00115 udu[j] += du; // " 00116 du = udu[j] - u[j]; // " 00117 00118 00119 vect fudu( (*f)(udu) ); 00120 00121 for(integer i=0;i<dims;i++) 00122 { 00123 jac[i][j]=(fudu[i]-fu[i])/du; 00124 } 00125 } 00126 return jac; 00127 } 00128 00129 template <integer dims, class NT> 00130 typename NT::matrix 00131 Jacobian<dims,NT>::calculate_accurate(const vect& u, const vect& fu) 00132 { 00133 // There might be a faster way to do this, but I haven found it yet 00134 matrix jac(0.); 00135 00136 for(integer j=0;j<dims;j++) 00137 { 00138 vect updu(u),umdu(u); 00139 00140 numT du = epsilon*abs(u[j])/2.,dup,dum; 00141 00142 if(du==0.0) du = epsilon/2.; // avoid numerical error 00143 00144 updu[j] += du; // " 00145 dup = updu[j] - u[j]; // " 00146 00147 umdu[j] -= du; // " 00148 dum = u[j] - umdu[j]; // " 00149 00150 vect fupdu( (*f)(updu) ); 00151 vect fumdu( (*f)(umdu) ); 00152 00153 for(integer i=0;i<dims;i++) 00154 { 00155 jac[i][j]=(fupdu[i]-fumdu[i])/(dum+dup); 00156 } 00157 } 00158 return jac; 00159 } 00160 00161 00162 } // end namespace 00163 #endif 00164 00165 /********************************************************************* 00166 $Id: jacobian.h,v 1.1 2001/05/22 10:54:55 mpeeters Exp $ 00167 ********************************************************************** 00168 00169 $Log: jacobian.h,v $ 00170 Revision 1.1 2001/05/22 10:54:55 mpeeters 00171 Moved sources and headers for libModel to model/ subdirectory, in an attempt to rationalize the source tree. This should make things "netter". 00172 00173 Revision 1.3 2001/05/21 11:53:16 mpeeters 00174 Removed Makefile.in, which is automatically generated anyway. 00175 00176 Revision 1.2 2000/09/15 10:26:31 mpeeters 00177 Cleaned out header and added CVS tails to files, removed superfuous 00178 @author comments, inserted dates 00179 00180 *********************************************************************/
More Info? Michael Peeters. Also, check our research website: www.alna.vub.ac.be
Last update: June 2002.
Looking for Open Source? Check out
!