check/series_expansion.cpp

   1 // check/series_expansion.cpp
   2
   3 /* Series expansion test (Laurent and Taylor series). */
   4
   5 #include <GiNaC/ginac.h>
   6
   7 static symbol x("x");
   8
   9 static unsigned check_series(const ex &e, const ex &point, const ex &d)
  10 {
  11         ex es = e.series(x, point, 8);
  12         ex ep = static_cast<series *>(es.bp)->convert_to_poly();
  13         if ((ep - d).compare(exZERO()) != 0) {
  14                 clog << "series expansion of " << e << " at " << point
  15              << " erroneously returned " << ep << " (instead of " << d
  16              << ")" << endl;
  17                 (ep-d).printtree(clog);
  18                 return 1;
  19         }
  20         return 0;
  21 }
  22
  23 // Series expansion
  24 static unsigned series1(void)
  25 {
  26         unsigned result = 0;
  27         ex e, d;
  28
  29         e = sin(x);
  30         d = x - pow(x, 3) / 6 + pow(x, 5) / 120 - pow(x, 7) / 5040 + Order(pow(x, 8));
  31         result += check_series(e, exZERO(), d);
  32
  33         e = cos(x);
  34         d = 1 - pow(x, 2) / 2 + pow(x, 4) / 24 - pow(x, 6) / 720 + Order(pow(x, 8));
  35         result += check_series(e, exZERO(), d);
  36
  37         e = exp(x);
  38         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));
  39         result += check_series(e, exZERO(), d);
  40
  41         e = pow(1 - x, -1);
  42         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));
  43         result += check_series(e, exZERO(), d);
  44
  45         e = x + pow(x, -1);
  46         d = x + pow(x, -1);
  47         result += check_series(e, exZERO(), d);
  48
  49         e = x + pow(x, -1);
  50         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));
  51         result += check_series(e, exONE(), d);
  52
  53         e = pow(x + pow(x, 3), -1);
  54         d = pow(x, -1) - x + pow(x, 3) - pow(x, 5) + Order(pow(x, 7));
  55         result += check_series(e, exZERO(), d);
  56
  57         e = pow(pow(x, 2) + pow(x, 4), -1);
  58         d = pow(x, -2) - 1 + pow(x, 2) - pow(x, 4) + Order(pow(x, 6));
  59         result += check_series(e, exZERO(), d);
  60
  61         e = pow(sin(x), -2);
  62         d = pow(x, -2) + numeric(1,3) + pow(x, 2) / 15 + pow(x, 4) * 2/189 + Order(pow(x, 5));
  63         result += check_series(e, exZERO(), d);
  64
  65         e = sin(x) / cos(x);
  66         d = x + pow(x, 3) / 3 + pow(x, 5) * 2/15 + pow(x, 7) * 17/315 + Order(pow(x, 8));
  67         result += check_series(e, exZERO(), d);
  68
  69         e = cos(x) / sin(x);
  70         d = pow(x, -1) - x / 3 - pow(x, 3) / 45 - pow(x, 5) * 2/945 + Order(pow(x, 6));
  71         result += check_series(e, exZERO(), d);
  72
  73         e = pow(numeric(2), x);
  74         ex t = log(ex(2)) * x;
  75         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));
  76         result += check_series(e, exZERO(), d.expand());
  77
  78         e = pow(Pi, x);
  79         t = log(Pi) * x;
  80         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));
  81         result += check_series(e, exZERO(), d.expand());
  82
  83         return result;
  84 }
  85
  86 // Series addition
  87 static unsigned series2(void)
  88 {
  89         unsigned result = 0;
  90         ex e, d;
  91
  92         e = pow(sin(x), -1).series(x, exZERO(), 8) + pow(sin(-x), -1).series(x, exZERO(), 12);
  93         d = Order(pow(x, 6));
  94         result += check_series(e, exZERO(), d);
  95
  96         return result;
  97 }
  98
  99 // Series multiplication
 100 static unsigned series3(void)
 101 {
 102         unsigned result = 0;
 103         ex e, d;
 104
 105         e = sin(x).series(x, exZERO(), 8) * pow(sin(x), -1).series(x, exZERO(), 12);
 106         d = 1 + Order(pow(x, 7));
 107         result += check_series(e, exZERO(), d);
 108
 109         return result;
 110 }
 111
 112 unsigned series_expansion(void)
 113 {
 114         unsigned result = 0;
 115
 116         cout << "checking series expansion..." << flush;
 117         clog << "---------series expansion:" << endl;
 118
 119         result += series1();
 120         result += series2();
 121         result += series3();
 122
 123         if (!result) {
 124                 cout << " passed ";
 125                 clog << "(no output)" << endl;
 126         } else {
 127                 cout << " failed ";
 128         }
 129         return result;
 130 }