SyFi 0.3
|
00001 // Copyright (C) 2006-2009 Kent-Andre Mardal and Simula Research Laboratory. 00002 // Licensed under the GNU GPL Version 2, or (at your option) any later version. 00003 00004 #include "diff_tools.h" 00005 #include "symbol_factory.h" 00006 00007 #include <stdexcept> 00008 00009 namespace SyFi 00010 { 00011 00012 /* Alternative implementation for any vector representation 00013 GiNaC::ex div(GiNaC::ex v){ 00014 using SyFi::nsd; 00015 using SyFi::p; 00016 00017 if(nsd != v.nops()) 00018 { 00019 throw std::invalid_argument("In div(v): The number of elements must equal nsd."); 00020 } 00021 00022 GiNaC::ex ret = 0; 00023 for(int i=0; i<nsd; i++) 00024 { 00025 ret = ret + v.op(i).diff(p[i]); 00026 } 00027 return ret; 00028 } 00029 */ 00030 00031 GiNaC::ex div(GiNaC::ex v) 00032 { 00033 using SyFi::nsd; 00034 using SyFi::x; 00035 using SyFi::y; 00036 using SyFi::z; 00037 00038 GiNaC::ex ret; 00039 if (GiNaC::is_a<GiNaC::matrix>(v)) 00040 { 00041 GiNaC::matrix m = GiNaC::ex_to<GiNaC::matrix>(v); 00042 if ( m.cols() == 1 && m.rows() == nsd ) 00043 { 00044 if (nsd == 1) 00045 { 00046 ret = diff(m,x); 00047 } 00048 else if (nsd == 2) 00049 { 00050 ret = diff(m.op(0),x) + diff(m.op(1),y) ; 00051 } 00052 else if (nsd == 3) 00053 { 00054 ret = diff(m.op(0),x) + diff(m.op(1),y) + diff(m.op(2),z) ; 00055 } 00056 else 00057 { 00058 throw std::runtime_error("Invalid nsd"); 00059 } 00060 00061 } 00062 else 00063 { 00064 GiNaC::matrix retm = GiNaC::matrix(m.cols(),1); 00065 if ( nsd != m.rows() ) 00066 { 00067 throw(std::invalid_argument("The number of rows must equal nsd.")); 00068 } 00069 GiNaC::symbol xr; 00070 GiNaC::ex tmp; 00071 for (unsigned int c=0; c<m.cols(); c++) 00072 { 00073 for (unsigned int r=0; r<m.rows(); r++) 00074 { 00075 if (r+1 == 1) xr = x; 00076 if (r+1 == 2) xr = y; 00077 if (r+1 == 3) xr = z; 00078 retm(c,0) += diff(m(c,r), xr); 00079 } 00080 } 00081 ret = retm; 00082 } 00083 return ret; 00084 00085 } 00086 else if (GiNaC::is_a<GiNaC::lst>(v)) 00087 { 00088 GiNaC::lst l = GiNaC::ex_to<GiNaC::lst>(v); 00089 return div(l); 00090 } 00091 throw std::invalid_argument("v must be a matrix or lst."); 00092 } 00093 00094 GiNaC::ex div(GiNaC::ex v, GiNaC::ex G) 00095 { 00096 using SyFi::nsd; 00097 using SyFi::x; 00098 using SyFi::y; 00099 using SyFi::z; 00100 00101 GiNaC::ex ret; 00102 if (GiNaC::is_a<GiNaC::matrix>(v) && GiNaC::is_a<GiNaC::matrix>(G)) 00103 { 00104 GiNaC::matrix m = GiNaC::ex_to<GiNaC::matrix>(v); 00105 GiNaC::matrix GG = GiNaC::ex_to<GiNaC::matrix>(G); 00106 if ( m.cols() == 1 && m.rows() == nsd && GG.rows() == nsd && GG.cols() == nsd ) 00107 { 00108 if ( nsd == 1 || nsd == 2 || nsd == 3) 00109 { 00110 ret = GiNaC::numeric(0); 00111 GiNaC::symbol xj; 00112 for (unsigned int i=0; i< nsd; i++) 00113 { 00114 for (unsigned int j=0; j< nsd; j++) 00115 { 00116 if (j == 0) xj = x; 00117 if (j == 1) xj = y; 00118 if (j == 2) xj = z; 00119 ret += m.op(i).diff(xj)*GG(i,j); 00120 } 00121 } 00122 } 00123 else 00124 { 00125 throw std::runtime_error("Invalid nsd"); 00126 } 00127 } 00128 else 00129 { 00130 throw std::invalid_argument("This functions needs v and G on the form: v.cols()=1, v.rows()=G.rows()=G.cols()=nsd."); 00131 } 00132 } 00133 else if (GiNaC::is_a<GiNaC::lst>(v)) 00134 { 00135 GiNaC::lst l = GiNaC::ex_to<GiNaC::lst>(v); 00136 return div(l,G); 00137 } 00138 else 00139 { 00140 throw std::invalid_argument("v must be a matrix or lst."); 00141 } 00142 return ret; 00143 } 00144 00145 GiNaC::ex div(GiNaC::lst& v) 00146 { 00147 using SyFi::x; 00148 using SyFi::y; 00149 using SyFi::z; 00150 00151 using SyFi::nsd; 00152 nsd = v.nops(); 00153 GiNaC::ex ret; 00154 if (nsd == 1) 00155 { 00156 ret = v.op(0).diff(x); 00157 } 00158 else if (nsd == 2) 00159 { 00160 ret = v.op(0).diff(x) + v.op(1).diff(y); 00161 } 00162 else if (nsd == 3) 00163 { 00164 ret = v.op(0).diff(x) + v.op(1).diff(y) + v.op(2).diff(z); 00165 } 00166 return ret; 00167 } 00168 00169 GiNaC::ex div(GiNaC::lst& v, GiNaC::ex G) 00170 { 00171 using SyFi::x; 00172 using SyFi::y; 00173 using SyFi::z; 00174 00175 using SyFi::nsd; 00176 nsd = v.nops(); 00177 GiNaC::ex ret; 00178 if (GiNaC::is_a<GiNaC::matrix>(G)) 00179 { 00180 GiNaC::matrix GG = GiNaC::ex_to<GiNaC::matrix>(G); 00181 if ( nsd != GG.cols() || nsd != GG.rows()) 00182 { 00183 throw(std::invalid_argument("The number of rows and cols in G must equal the size of v.")); 00184 } 00185 if (nsd == 1 || nsd == 2 || nsd == 3 ) 00186 { 00187 GiNaC::symbol xj; 00188 ret = GiNaC::numeric(0); 00189 for (unsigned int i=0; i< nsd; i++) 00190 { 00191 for (unsigned int j=0; j< nsd; j++) 00192 { 00193 if (i == 0) xj = x; 00194 if (i == 1) xj = y; 00195 if (i == 2) xj = z; 00196 ret += v.op(i).diff(xj)*GG(i,j); 00197 } 00198 } 00199 } 00200 else 00201 { 00202 throw std::runtime_error("Invalid nsd"); 00203 } 00204 } 00205 else 00206 { 00207 throw std::invalid_argument("v must be a matrix."); 00208 } 00209 return ret; 00210 } 00211 00212 GiNaC::ex div(GiNaC::exvector& v) 00213 { 00214 using SyFi::nsd; 00215 using SyFi::x; 00216 using SyFi::y; 00217 using SyFi::z; 00218 00219 GiNaC::ex ret; 00220 if (nsd == 2) 00221 { 00222 ret = v[0].diff(x) + v[1].diff(y); 00223 } 00224 else if (nsd == 3) 00225 { 00226 ret = v[0].diff(x) + v[1].diff(y) + v[2].diff(z); 00227 } 00228 return ret; 00229 } 00230 00231 GiNaC::ex grad(GiNaC::ex f) 00232 { 00233 using SyFi::nsd; 00234 using SyFi::x; 00235 using SyFi::y; 00236 using SyFi::z; 00237 00238 if (GiNaC::is_a<GiNaC::matrix>(f)) 00239 { 00240 GiNaC::matrix m = GiNaC::ex_to<GiNaC::matrix>(f); 00241 GiNaC::matrix ret_m(nsd,m.rows()); 00242 for (unsigned int r=0; r< m.rows(); r++) 00243 { 00244 if (nsd == 1) 00245 { 00246 // ret_m(0,r) = diff(m.op(r),x); 00247 return diff(f, x); 00248 } 00249 else if ( nsd == 2) 00250 { 00251 ret_m(0,r) = diff(m.op(r),x); 00252 ret_m(1,r) = diff(m.op(r),y); 00253 } 00254 else if ( nsd == 3) 00255 { 00256 ret_m(0,r) = diff(m.op(r),x); 00257 ret_m(1,r) = diff(m.op(r),y); 00258 ret_m(2,r) = diff(m.op(r),z); 00259 } 00260 } 00261 return ret_m; 00262 } 00263 else 00264 { 00265 00266 if (nsd == 1) 00267 { 00268 // return GiNaC::matrix(nsd,1,GiNaC::lst(diff(f,x))); 00269 return diff(f,x); 00270 } 00271 else if ( nsd == 2) 00272 { 00273 return GiNaC::matrix(nsd,1,GiNaC::lst(diff(f,x), diff(f,y))); 00274 } 00275 else if ( nsd == 3) 00276 { 00277 return GiNaC::matrix(nsd,1,GiNaC::lst(diff(f,x), diff(f,y), diff(f,z))); 00278 } 00279 else 00280 { 00281 throw(std::invalid_argument("nsd must be either 1, 2, or 3.")); 00282 return GiNaC::matrix(); 00283 } 00284 } 00285 } 00286 00287 GiNaC::ex grad(GiNaC::ex f, GiNaC::ex G) 00288 { 00289 using SyFi::nsd; 00290 using SyFi::x; 00291 using SyFi::y; 00292 using SyFi::z; 00293 00294 GiNaC::symbol xr; 00295 if ( GiNaC::is_a<GiNaC::matrix>(G)) 00296 { 00297 GiNaC::matrix GG = GiNaC::ex_to<GiNaC::matrix>(G); 00298 00299 if (! (GG.rows() == nsd && GG.cols() == nsd )) 00300 { 00301 throw(std::invalid_argument("The number of cols/rows in G must equal nsd.")); 00302 } 00303 00304 if (GiNaC::is_a<GiNaC::matrix>(f) ) 00305 { 00306 GiNaC::matrix m = GiNaC::ex_to<GiNaC::matrix>(f); 00307 GiNaC::matrix ret_m(nsd,m.rows()); 00308 for (unsigned int k=0; k< m.rows(); k++) 00309 { 00310 for (unsigned int c=0; c<nsd; c++) 00311 { 00312 for (unsigned int r=0; r<nsd; r++) 00313 { 00314 if (r == 0) xr = x; 00315 if (r == 1) xr = y; 00316 if (r == 2) xr = z; 00317 ret_m(c,k) += diff(f,xr)*GG(r,c); 00318 } 00319 } 00320 } 00321 00322 return ret_m; 00323 } 00324 else 00325 { 00326 GiNaC::matrix ret_m(nsd,1); 00327 for (unsigned int c=0; c<nsd; c++) 00328 { 00329 for (unsigned int r=0; r<nsd; r++) 00330 { 00331 if (r == 0) xr = x; 00332 if (r == 1) xr = y; 00333 if (r == 2) xr = z; 00334 ret_m(c,0) += diff(f,xr)*GG(r,c); 00335 } 00336 } 00337 return ret_m; 00338 } 00339 } 00340 else 00341 { 00342 throw(std::invalid_argument("G must be a matrix.")); 00343 } 00344 } 00345 00346 } // namespace SyFi