00001
00002
00003
00004 #include "diff_tools.h"
00005 #include "symbol_factory.h"
00006
00007 #include <stdexcept>
00008
00009 namespace SyFi
00010 {
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
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
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
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 }