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