SyFi 0.3
diff_tools.cpp
Go to the documentation of this file.
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines