Octave changeset

[Craig Hudson]

# HG changeset patch
# User Craig Hudson <address@hidden>
# Date 1371220620 -3600
# Node ID 7620295b392a8752cb368915f98d98a5e535d3f3
# Parent  169ea45f260fb0442e750779df1e66c603fd947b
Alter colloc to handle Laguerre polynomials (see Numerical, items 5,6 on Projects page)

diff -r 169ea45f260f -r 7620295b392a libinterp/corefcn/colloc.cc
--- a/libinterp/corefcn/colloc.cc    Mon Jun 10 16:44:19 2013 +0800
+++ b/libinterp/corefcn/colloc.cc    Fri Jun 14 15:37:00 2013 +0100
@@ -36,18 +36,19 @@
 
 DEFUN (colloc, args, ,
   "-*- texinfo -*-\n\
address@hidden {Built-in Function} address@hidden, @var{amat}, @var{bmat}, @var{q}] =} colloc (@var{n}, \"left\", \"right\")\n\
address@hidden {Built-in Function} address@hidden, @var{amat}, @var{bmat}, @var{q}] =} colloc (@var{n}, \"left\", \"right\", \"legendre\"/\"laguerre\")\n\
 Compute derivative and integral weight matrices for orthogonal\n\
 collocation using the subroutines given in J. Villadsen and\n\
 M. L. Michelsen, @cite{Solution of Differential Equation Models by\n\
 Polynomial Approximation}.\n\
+\"legendre\" (default) or \"laguerre\" will select the appropriate polynomial type.\
 @end deftypefn")
 {
   octave_value_list retval;
 
   int nargin = args.length ();
 
-  if (nargin < 1 || nargin > 3)
+  if (nargin < 1 || nargin > 4)
     {
       print_usage ();
       return retval;
@@ -78,6 +79,7 @@
     }
 
   octave_idx_type ntot = ncol;
+  enum polynomial_type polynomial_type = legendre;
   octave_idx_type left = 0;
   octave_idx_type right = 0;
 
@@ -103,6 +105,14 @@
             {
               left = 1;
             }
+          else if (s == "legendre")
+            {
+              polynomial_type = legendre;
+            }
+          else if (s == "laguerre")
+            {
+              polynomial_type = laguerre;
+            }
           else
             {
               error ("colloc: unrecognized argument");
@@ -123,7 +133,7 @@
       return retval;
     }
 
-  CollocWt wts (ncol, left, right);
+  CollocWt wts (ncol, left, right, polynomial_type);
 
   ColumnVector r = wts.roots ();
   Matrix A = wts.first ();
diff -r 169ea45f260f -r 7620295b392a liboctave/numeric/CollocWt.cc
--- a/liboctave/numeric/CollocWt.cc    Mon Jun 10 16:44:19 2013 +0800
+++ b/liboctave/numeric/CollocWt.cc    Fri Jun 14 15:37:00 2013 +0100
@@ -53,8 +53,8 @@
 // They may be computed using jcobi.
 
 static void
-dif (octave_idx_type nt, double *root, double *dif1, double *dif2,
-     double *dif3)
+dif (octave_idx_type nt, double *root,
+     double *dif1, double *dif2, double *dif3)
 {
   // Evaluate derivatives of node polynomial using recursion formulas.
 
@@ -80,11 +80,14 @@
     }
 }
 
-// Compute the zeros of the Jacobi polynomial.
+// Compute the zeros of the Jacobi or Laguerre polynomial.
 //
-//    (alpha,beta)
-//   p  (x)
-//    n
+//    (alpha,beta)          (alpha)
+//   p  (x)          or    L  (x)
+//    n                     n
+//
+// NB function name could do with changing, now that it doesn't
+// compute exclusively Jacobi polynomials.
 //
 // Use dif to compute the derivatives of the node
 // polynomial
@@ -93,7 +96,7 @@
 //   p  (x)  =  (x)   *  p (x)         *  (1 - x)
 //    nt                   n
 //
-// at the interpolation points.
+// (or similar for Laguerre) at the interpolation points.
 //
 // See Villadsen and Michelsen, pages 131-132 and 418.
 //
@@ -104,6 +107,13 @@
 //   n      : the degree of the jacobi polynomial, (i.e. the number
 //            of interior interpolation points)
 //
+//   polynomial_type : the type of polynomial used in the collocation.
+//                     Jacobi and generalised Laguerre polynomials can be
+//                     calculated, but the colloc function currently only
+//                     passes (shifted) Legendre and non-generalised
+//                     Laguerre polynomials, i.e. the special cases where
+//                     alpha == beta == 0.
+//
 //   n0     : determines whether x = 0 is included as an
 //            interpolation point
 //
@@ -141,7 +151,8 @@
 //            of the node polynomial at the zeros
 
 static bool
-jcobi (octave_idx_type n, octave_idx_type n0, octave_idx_type n1,
+jcobi (octave_idx_type n, enum polynomial_type polynomial_type,
+       octave_idx_type n0, octave_idx_type n1,
        double alpha, double beta, double *dif1, double *dif2,
        double *dif3, double *root)
 {
@@ -153,40 +164,55 @@
   assert (nt > 1);
 
 // -- first evaluation of coefficients in recursion formulas.
-// -- recursion coefficients are stored in dif1 and dif2.
+// -- recursion coefficients are stored in dif1 and dif2
+// -- (additionally dif3 for Laguerre polynomials).
 
-  double ab = alpha + beta;
-  double ad = beta - alpha;
-  double ap = beta * alpha;
+  if (polynomial_type == laguerre)
+    {
+      for (octave_idx_type i = 0; i < n; i++)
+        {
+          dif1[i] = 2.0 * i + alpha + 1.0;
+          dif2[i] = i + alpha;
+          dif3[i] = 1.0 / (i + 1.0);
+        }
+    }
+  else // polynomial_type == legendre
+    {
+      double ab = alpha + beta;
+      double ad = beta - alpha;
+      double ap = beta * alpha;
 
-  dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0;
-  dif2[0] = 0.0;
+      dif1[0] = (ad / (ab + 2.0) + 1.0) / 2.0;
+      dif2[0] = 0.0;
 
-  if (n >= 2)
-    {
-      for (octave_idx_type i = 1; i < n; i++)
+      if (n >= 2)
         {
-          double z1 = i;
-          double z = ab + 2 * z1;
+          for (octave_idx_type i = 1; i < n; i++)
+            {
+              double z1 = i;
+              double z = ab + 2 * z1;
 
-          dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0;
+              dif1[i] = (ab * ad / z / (z + 2.0) + 1.0) / 2.0;
 
-          if (i == 1)
-            dif2[i] = (ab + ap + z1) / z / z / (z + 1.0);
-          else
-            {
-              z = z * z;
-              double y = z1 * (ab + z1);
-              y = y * (ap + y);
-              dif2[i] = y / z / (z - 1.0);
+              if (i == 1)
+                dif2[i] = (ab + ap + z1) / z / z / (z + 1.0);
+              else
+                {
+                  z = z * z;
+                  double y = z1 * (ab + z1);
+                  y = y * (ap + y);
+                  dif2[i] = y / z / (z - 1.0);
+                }
             }
         }
     }
 
   // Root determination by Newton method with suppression of previously
   // determined roots.
+  // Recursion relation for Ln'(0) fails (division by zero), hence different
+  // start point for Laguerre polynomials
 
-  double x = 0.0;
+  double x = (polynomial_type == laguerre ? 0.01 : 0.0);
 
   for (octave_idx_type i = 0; i < n; i++)
     {
@@ -203,8 +229,18 @@
 
           for (octave_idx_type j = 0; j < n; j++)
             {
-              double xp  = (dif1[j] - x) * xn  - dif2[j] * xd;
-              double xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn;
+              double xp, xp1;
+
+              if (polynomial_type == laguerre)
+                {
+                  xp = dif3[j] * ((dif1[j] - x ) * xn - dif2[j] * xd);
+                  xp1 = ((j + 1) * xp - (dif2[j] + 1.0) * xn) / x;
+                }
+              else // polynomial_type == legendre
+                {
+                  xp  = (dif1[j] - x) * xn  - dif2[j] * xd;
+                  xp1 = (dif1[j] - x) * xn1 - dif2[j] * xd1 - xn;
+                }
 
               xd  = xn;
               xd1 = xn1;
@@ -267,6 +303,9 @@
 //   n      : the degree of the jacobi polynomial, (i.e. the number
 //            of interior interpolation points)
 //
+//   polynomial_type : type of polynomial used in the collocation,
+//                     (shifted) Legendre or Laguerre
+//
 //   n0     : determines whether x = 0 is included as an
 //            interpolation point
 //
@@ -303,7 +342,8 @@
 //   vect   : one dimensional vector of computed weights
 
 static void
-dfopr (octave_idx_type n, octave_idx_type n0, octave_idx_type n1,
+dfopr (octave_idx_type n, enum polynomial_type polynomial_type,
+       octave_idx_type n0, octave_idx_type n1,
        octave_idx_type i, octave_idx_type id, double *dif1,
        double *dif2, double *dif3, double *root, double *vect)
 {
@@ -352,12 +392,12 @@
         {
           double x  = root[j];
 
-          double ax = x * (1.0 - x);
+          double ax = x * (polynomial_type == laguerre ? 1.0 : (1.0 - x));
 
           if (n0 == 0)
             ax = ax / x / x;
 
-          if (n1 == 0)
+          if (polynomial_type != laguerre && n1 == 0)
             ax = ax / (1.0 - x) / (1.0 - x);
 
           vect[j] = ax / (dif1[j] * dif1[j]);
@@ -418,7 +458,15 @@
       return;
     }
 
-  octave_idx_type nt = n + inc_left + inc_right;
+  // Ignore "right" parameter for Laguerre polynomials, as boundary points
+  // are 0 and Inf, not 0 and 1.
+
+  if (polynomial_type == laguerre)
+    {
+      inc_right = 0;
+    }
+
+    octave_idx_type nt = n + inc_left + inc_right;
 
   if (nt < 0)
     {
@@ -427,6 +475,11 @@
     }
   else if (nt == 0)
     return;
+  else if (nt == 1)
+    {
+      error ("total number of collocation points equal to one");
+      return;
+    } // assert (nt > 1) in jcobi and dfopr results in Octave crash
 
   Array<double> dif1 (dim_vector (nt, 1));
   double *pdif1 = dif1.fortran_vec ();
@@ -449,7 +502,8 @@
 
   // Compute roots.
 
-  if (! jcobi (n, inc_left, inc_right, Alpha, Beta, pdif1, pdif2, pdif3, pr))
+  if (! jcobi (n, polynomial_type, inc_left, inc_right, Alpha, Beta, pdif1,
+               pdif2, pdif3, pr))
     {
       error ("jcobi: newton iteration failed");
       return;
@@ -462,7 +516,8 @@
   id = 1;
   for (octave_idx_type i = 0; i < nt; i++)
     {
-      dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
+      dfopr (n, polynomial_type, inc_left, inc_right, i, id,
+             pdif1, pdif2, pdif3, pr, pvect);
 
       for (octave_idx_type j = 0; j < nt; j++)
         A(i,j) = vect(j);
@@ -473,7 +528,8 @@
   id = 2;
   for (octave_idx_type i = 0; i < nt; i++)
     {
-      dfopr (n, inc_left, inc_right, i, id, pdif1, pdif2, pdif3, pr, pvect);
+      dfopr (n, polynomial_type, inc_left, inc_right, i, id,
+             pdif1, pdif2, pdif3, pr, pvect);
 
       for (octave_idx_type j = 0; j < nt; j++)
         B(i,j) = vect(j);
@@ -483,7 +539,8 @@
 
   id = 3;
   double *pq = q.fortran_vec ();
-  dfopr (n, inc_left, inc_right, id, id, pdif1, pdif2, pdif3, pr, pq);
+  dfopr (n, polynomial_type, inc_left, inc_right, id, id,
+         pdif1, pdif2, pdif3, pr, pq);
 
   initialized = 1;
 }
diff -r 169ea45f260f -r 7620295b392a liboctave/numeric/CollocWt.h
--- a/liboctave/numeric/CollocWt.h    Mon Jun 10 16:44:19 2013 +0800
+++ b/liboctave/numeric/CollocWt.h    Fri Jun 14 15:37:00 2013 +0100
@@ -28,6 +28,9 @@
 #include "dMatrix.h"
 #include "dColVector.h"
 
+enum polynomial_type {legendre, laguerre};
+
+
 class
 OCTAVE_API
 CollocWt
@@ -35,31 +38,38 @@
 public:
 
   CollocWt (void)
-    : n (0), inc_left (0), inc_right (0), lb (0.0), rb (1.0),
-      Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { }
+    : n (0), polynomial_type (legendre), inc_left (0), inc_right (0),
+      lb (0.0), rb (1.0), Alpha (0.0), Beta (0.0), r (), q (), A (),
+      B (), initialized (false) { }
 
-  CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir)
-    : n (nc), inc_left (il), inc_right (ir), lb (0.0), rb (1.0),
-      Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { }
+  CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir,
+            enum polynomial_type pt = legendre)
+    : n (nc), polynomial_type (pt), inc_left (il), inc_right (ir), lb (0.0),
+      rb (1.0), Alpha (0.0), Beta (0.0), r (), q (), A (), B (),
+      initialized (false) { }
 
   CollocWt (octave_idx_type nc, octave_idx_type il, octave_idx_type ir,
             double l, double rr)
-    : n (nc), inc_left (il), inc_right (ir), lb (l), rb (rr),
-      Alpha (0.0), Beta (0.0), r (), q (), A (), B (), initialized (false) { }
+    : n (nc), polynomial_type (legendre), inc_left (il), inc_right (ir),
+      lb (l), rb (rr), Alpha (0.0), Beta (0.0), r (), q (), A (), B (),
+      initialized (false) { }
 
   CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il,
             octave_idx_type ir)
-    : n (nc), inc_left (il), inc_right (ir), lb (0.0), rb (1.0),
-      Alpha (a), Beta (b), r (), q (), A (), B (), initialized (false) { }
+    : n (nc), polynomial_type (legendre), inc_left (il), inc_right (ir),
+      lb (0.0), rb (1.0), Alpha (a), Beta (b), r (), q (), A (), B (),
+      initialized (false) { }
 
   CollocWt (octave_idx_type nc, double a, double b, octave_idx_type il,
             octave_idx_type ir,
                       double ll, double rr)
-    : n (nc), inc_left (il), inc_right (ir), lb (ll), rb (rr),
-      Alpha (a), Beta (b), r (), q (), A (), B (), initialized (false) { }
+    : n (nc), polynomial_type (legendre), inc_left (il), inc_right (ir),
+      lb (ll), rb (rr), Alpha (a), Beta (b), r (), q (), A (), B (),
+      initialized (false) { }
 
   CollocWt (const CollocWt& a)
-    : n (a.n), inc_left (a.inc_left), inc_right (a.inc_right),
+    : n (a.n), polynomial_type (legendre),
+      inc_left (a.inc_left), inc_right (a.inc_right),
       lb (a.lb), rb (a.rb), Alpha (a.Alpha), Beta (a.Beta),
       r (a.r), q (a.q), A (a.A), B (a.B),
       initialized (a.initialized) { }
@@ -69,6 +79,7 @@
       if (this != &a)
         {
           n = a.n;
+          polynomial_type = a.polynomial_type;
           inc_left = a.inc_left;
           inc_right = a.inc_right;
           lb = a.lb;
@@ -165,6 +176,8 @@
 
   octave_idx_type n;
 
+  enum polynomial_type polynomial_type;
+
   octave_idx_type inc_left;
   octave_idx_type inc_right;