check/series_expansion.cpp

   1 /** @file series_expansion.cpp
   2  *
   3  *  Series expansion test (Laurent and Taylor series). */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <ginac/ginac.h>
  24 using namespace GiNaC;
  25
  26 static symbol x("x");
  27
  28 static unsigned check_series(const ex &e, const ex &point, const ex &d)
  29 {
  30         ex es = e.series(x, point, 8);
  31         ex ep = static_cast<series *>(es.bp)->convert_to_poly();
  32         if ((ep - d).compare(exZERO()) != 0) {
  33                 clog << "series expansion of " << e << " at " << point
  34              << " erroneously returned " << ep << " (instead of " << d
  35              << ")" << endl;
  36                 (ep-d).printtree(clog);
  37                 return 1;
  38         }
  39         return 0;
  40 }
  41
  42 // Series expansion
  43 static unsigned series1(void)
  44 {
  45         unsigned result = 0;
  46         ex e, d;
  47
  48         e = sin(x);
  49         d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
  50         result += check_series(e, exZERO(), d);
  51
  52         e = cos(x);
  53         d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
  54         result += check_series(e, exZERO(), d);
  55
  56         e = exp(x);
  57         d = 1 + x + pow(x, 2) / 2 + pow(x, 3) / 6 + pow(x, 4) / 24 + pow(x, 5) / 120 + pow(x, 6) / 720 + pow(x, 7) / 5040 + Order(pow(x, 8));
  58         result += check_series(e, exZERO(), d);
  59
  60         e = pow(1 - x, -1);
  61         d = 1 + x + pow(x, 2) + pow(x, 3) + pow(x, 4) + pow(x, 5) + pow(x, 6) + pow(x, 7) + Order(pow(x, 8));
  62         result += check_series(e, exZERO(), d);
  63
  64         e = x + pow(x, -1);
  65         d = x + pow(x, -1);
  66         result += check_series(e, exZERO(), d);
  67
  68         e = x + pow(x, -1);
  69         d = 2 + pow(x-1, 2) - pow(x-1, 3) + pow(x-1, 4) - pow(x-1, 5) + pow(x-1, 6) - pow(x-1, 7) + Order(pow(x-1, 8));
  70         result += check_series(e, exONE(), d);
  71
  72         e = pow(x + pow(x, 3), -1);
  73         d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
  74         result += check_series(e, exZERO(), d);
  75
  76         e = pow(pow(x, 2) + pow(x, 4), -1);
  77         d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
  78         result += check_series(e, exZERO(), d);
  79
  80         e = pow(sin(x), -2);
  81         d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
  82         result += check_series(e, exZERO(), d);
  83
  84         e = sin(x) / cos(x);
  85         d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
  86         result += check_series(e, exZERO(), d);
  87
  88         e = cos(x) / sin(x);
  89         d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
  90         result += check_series(e, exZERO(), d);
  91
  92         e = pow(numeric(2), x);
  93         ex t = log(ex(2)) * x;
  94         d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
  95         result += check_series(e, exZERO(), d.expand());
  96
  97         e = pow(Pi, x);
  98         t = log(Pi) * x;
  99         d = 1 + t + pow(t, 2) / 2 + pow(t, 3) / 6 + pow(t, 4) / 24 + pow(t, 5) / 120 + pow(t, 6) / 720 + pow(t, 7) / 5040 + Order(pow(x, 8));
 100         result += check_series(e, exZERO(), d.expand());
 101
 102         return result;
 103 }
 104
 105 // Series addition
 106 static unsigned series2(void)
 107 {
 108         unsigned result = 0;
 109         ex e, d;
 110
 111         e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
 112         d = Order(pow(x, 6));
 113         result += check_series(e, exZERO(), d);
 114
 115         return result;
 116 }
 117
 118 // Series multiplication
 119 static unsigned series3(void)
 120 {
 121         unsigned result = 0;
 122         ex e, d;
 123
 124         e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
 125         d = 1 + Order(pow(x, 7));
 126         result += check_series(e, exZERO(), d);
 127
 128         return result;
 129 }
 130
 131 unsigned series_expansion(void)
 132 {
 133         unsigned result = 0;
 134
 135         cout << "checking series expansion..." << flush;
 136         clog << "---------series expansion:" << endl;
 137
 138         result += series1();
 139         result += series2();
 140         result += series3();
 141
 142         if (!result) {
 143                 cout << " passed ";
 144                 clog << "(no output)" << endl;
 145         } else {
 146                 cout << " failed ";
 147         }
 148         return result;
 149 }