asy 1.79

git-svn-id: svn://tug.org/texlive/trunk@14078 c570f23f-e606-0410-a88d-b1316a301751

1 files changed, 188 insertions, 67 deletions

diff --git a/Build/source/utils/asymptote/base/ode.asy b/Build/source/utils/asymptote/base/ode.asy
index fe65f2f1c8e..77ab8d6e589 100644
--- a/Build/source/utils/asymptote/base/ode.asy
+++ b/Build/source/utils/asymptote/base/ode.asy
@@ -1,7 +1,8 @@
-real stepfactor=2.0; // Maximum dynamic step size adjustment factor.
+real stepfactor=2; // Maximum dynamic step size adjustment factor.
 
 struct coefficients
 {
+  real[] steps;
   real[] factors;
   real[][] weights;
   real[] highOrderWeights;
@@ -11,7 +12,6 @@ struct coefficients
 struct RKTableau
 {
   int order;
-  real[] steps;
   coefficients a;
   void stepDependence(real h, real c, coefficients a) {}
    
@@ -25,8 +25,8 @@ struct RKTableau
                          return sum(weights[i]);},weights.length),
                      void stepDependence(real, real, coefficients)=null) {
     this.order=order;
-    this.steps=steps;
-    a.factors=array(steps.length+1,1);
+    a.steps=steps;
+    a.factors=array(a.steps.length+1,1);
     a.weights=weights;
     a.highOrderWeights=highOrderWeights;
     a.lowOrderWeights=lowOrderWeights;
@@ -39,24 +39,113 @@ struct RKTableau
   }
 }
 
+real[] Coeff={1,1/2,1/6,1/24,1/120,1/720,1/5040,1/40320,1/362880,1/3628800,
+              1/39916800.0,1/479001600.0,1/6227020800.0,1/87178291200.0,
+              1/1307674368000.0,1/20922789888000.0,1/355687428096000.0,
+              1/6402373705728000.0,1/121645100408832000.0,
+              1/2432902008176640000.0,1/51090942171709440000.0,
+              1/1124000727777607680000.0};
+
+real phi1(real x) {return x != 0 ? expm1(x)/x : 1;}
+
+real phi2(real x)
+{
+  real x2=x*x;
+  if(fabs(x) > 1) return (exp(x)-x-1)/x2;
+  real x3=x2*x;
+  real x5=x2*x3;
+  if(fabs(x) < 0.1) 
+    return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5]
+      +x5*Coeff[6]+x3*x3*Coeff[7]+x5*x2*Coeff[8]+x5*x3*Coeff[9];
+    else {
+      real x7=x5*x2;
+      real x8=x7*x;
+      return Coeff[1]+x*Coeff[2]+x2*Coeff[3]+x3*Coeff[4]+x2*x2*Coeff[5]
+        +x5*Coeff[6]+x3*x3*Coeff[7]+x7*Coeff[8]+x8*Coeff[9]
+        +x8*x*Coeff[10]+x5*x5*Coeff[11]+x8*x3*Coeff[12]+x7*x5*Coeff[13]+
+        x8*x5*Coeff[14]+x7*x7*Coeff[15]+x8*x7*Coeff[16]+x8*x8*Coeff[17];
+    }
+}
+
+real phi3(real x)
+{
+  real x2=x*x;
+  real x3=x2*x;
+  if(fabs(x) > 1.6) return (exp(x)-0.5*x2-x-1)/x3;
+  real x5=x2*x3;
+  if(fabs(x) < 0.1) 
+    return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5]
+      +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9]
+      +x5*x3*Coeff[10];
+    else {
+      real x7=x5*x2;
+      real x8=x7*x;
+      real x16=x8*x8;
+      return Coeff[2]+x*Coeff[3]+x2*Coeff[4]+x3*Coeff[5]
+        +x2*x2*Coeff[6]+x5*Coeff[7]+x3*x3*Coeff[8]+x5*x2*Coeff[9]
+        +x5*x3*Coeff[10]+x8*x*Coeff[11]
+        +x5*x5*Coeff[12]+x8*x3*Coeff[13]+x7*x5*Coeff[14]
+        +x8*x5*Coeff[15]+x7*x7*Coeff[16]+x8*x7*Coeff[17]+x16*Coeff[18]
+        +x16*x*Coeff[19]+x16*x2*Coeff[20];
+    }
+}
+
+void expfactors(real x, coefficients a) 
+{
+  for(int i=0; i < a.steps.length; ++i)
+    a.factors[i]=exp(x*a.steps[i]);
+  a.factors[a.steps.length]=exp(x);
+}
+      
 // First-Order Euler
 RKTableau Euler=RKTableau(1,new real[][],
                           new real[] {1});
 
+// First-Order Exponential Euler
 RKTableau E_Euler=RKTableau(1,new real[][], new real[] {1},
-                            new void (real h, real c, coefficients a) {
+                            new void(real h, real c, coefficients a) {
                               real x=-c*h;
-                              a.factors[0]=exp(x);
-                              a.highOrderWeights[0]=x != 0 ? expm1(x)/x : 1;
+                              expfactors(x,a);
+                              a.highOrderWeights[0]=phi1(x);
                             });
 
 // Second-Order Runge-Kutta
 RKTableau RK2=RKTableau(2,new real[][] {{1/2}},
-                        new real[] {0,1});
+                        new real[] {0,1}, // 2nd order
+                        new real[] {1,0}); // 1st order
+
+// Second-Order Exponential Runge-Kutta
+RKTableau E_RK2=RKTableau(2,new real[][] {{1/2}},
+                          new real[] {0,1}, // 2nd order
+                          new real[] {1,0}, // 1st order
+                          new void(real h, real c, coefficients a) {
+                            real x=-c*h;
+                            expfactors(x,a);
+                            a.weights[0][0]=1/2*phi1(x/2);
+                            real w=phi1(x);
+                            a.highOrderWeights[0]=0;
+                            a.highOrderWeights[1]=w;
+                            a.lowOrderWeights[0]=w;
+                          });
 
 // Second-Order Predictor-Corrector
 RKTableau PC=RKTableau(2,new real[][] {{1}},
-                       new real[] {1/2,1/2});
+                       new real[] {1/2,1/2}, // 2nd order
+                       new real[] {1,0}); // 1st order
+
+// Second-Order Exponential Predictor-Corrector
+RKTableau E_PC=RKTableau(2,new real[][] {{1}},
+                         new real[] {1/2,1/2}, // 2nd order
+                         new real[] {1,0}, // 1st order
+                         new void(real h, real c, coefficients a) {
+                           real x=-c*h;
+                           expfactors(x,a);
+                           real w=phi1(x);
+                           a.weights[0][0]=w;
+                           a.highOrderWeights[0]=w/2;
+                           a.highOrderWeights[1]=w/2;
+                           a.lowOrderWeights[0]=w;
+                         });
 
 // Third-Order Classical Runge-Kutta
 RKTableau RK3=RKTableau(3,new real[][] {{1/2},{-1,2}},
@@ -67,6 +156,30 @@ RKTableau RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}},
                           new real[] {2/9,1/3,4/9}, // 3rd order
                           new real[] {7/24,1/4,1/3,1/8}); // 2nd order
 
+// Third-Order Exponential Bogacki-Shampine Runge-Kutta
+RKTableau E_RK3BS=RKTableau(3,new real[][] {{1/2},{0,3/4}},
+                            new real[] {2/9,1/3,4/9}, // 3rd order
+                            new real[] {7/24,1/4,1/3,1/8}, // 2nd order
+                            new void(real h, real c, coefficients a) {
+                              real x=-c*h;
+                              expfactors(x,a);
+                              real w=phi1(x);
+                              real w2=phi2(x);
+                              a.weights[0][0]=1/2*phi1(x/2);
+                              real a11=9/8*phi2(3/4*x)+3/8*phi2(x/2);
+                              a.weights[1][0]=3/4*phi1(3/4*x)-a11;
+                              a.weights[1][1]=a11;
+                              real a21=1/3*w;
+                              real a22=4/3*w2-2/9*w;
+                              a.highOrderWeights[0]=w-a21-a22;
+                              a.highOrderWeights[1]=a21;
+                              a.highOrderWeights[2]=a22;
+                              a.lowOrderWeights[0]=w-17/12*w2;
+                              a.lowOrderWeights[1]=w2/2;
+                              a.lowOrderWeights[2]=2/3*w2;
+                              a.lowOrderWeights[3]=w2/4;
+                            });
+
 // Fourth-Order Classical Runge-Kutta
 RKTableau RK4=RKTableau(4,new real[][] {{1/2},{0,1/2},{0,0,1}},
                         new real[] {1/6,1/3,1/3,1/6});
@@ -108,7 +221,7 @@ RKTableau RK5DP=RKTableau(5,new real[][] {{1/5},
 
 real error(real error, real initial, real lowOrder, real norm, real diff) 
 {
-  if(initial != 0.0 && lowOrder != initial) {
+  if(initial != 0 && lowOrder != initial) {
     static real epsilon=realMin/realEpsilon;
     real denom=max(abs(norm),abs(initial))+epsilon;
     return max(error,max(abs(diff)/denom));
@@ -116,46 +229,39 @@ real error(real error, real initial, real lowOrder, real norm, real diff)
   return error;
 }
 
-real adjust(real h, real error, real t, real c, real tolmin, real tolmax,
-            real dtmin, real dtmax, RKTableau tableau, bool verbose=true) 
+void report(real old, real h, real t)
 {
-  real dt=h;
-  void report(real t) {
-    if(h != dt) {
-      tableau.stepDependence(h,c,tableau.a);
-      if(verbose)
-        write("Time step changed from "+(string) dt+" to "+(string) h+" at t="+
-              (string) t+".");
-    }
-  }
-  if(error > tolmax) {
-    h=max(h*max((tolmin/error)^tableau.pshrink,1/stepfactor),dtmin);
-    report(t);
-    return h;
-  }
-  if(error > 0 && error < tolmin) {
-    h=min(h*min((tolmin/error)^tableau.pgrow,stepfactor),dtmax);
-    report(t+dt);
-  }
+  write("Time step changed from "+(string) old+" to "+(string) h+" at t="+
+        (string) t+".");
+}
+
+real adjust(real h, real error, real tolmin, real tolmax, RKTableau tableau)
+{
+  if(error > tolmax)
+    h *= max((tolmin/error)^tableau.pshrink,1/stepfactor);
+  else if(error > 0 && error < tolmin)
+    h *= min((tolmin/error)^tableau.pgrow,stepfactor);
   return h;
 }
 
 // Integrate dy/dt+cy=f(t,y) from a to b using initial conditions y,
 // specifying either the step size h or the number of steps n.
-real integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
-               real h=0, int n=0, bool dynamic=false, real tolmin=0,
-               real tolmax=0, real dtmin=0, real dtmax=realMax,
-               RKTableau tableau, bool verbose=false)
+real[] integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
+                 real h=0, int n=0, bool dynamic=false, real tolmin=0,
+                 real tolmax=0, real dtmin=0, real dtmax=realMax,
+                 RKTableau tableau, bool verbose=false)
 {
-  real f(real t, real y)=(c == 0 || tableau.exponential) ? g :
-    new real(real t, real y) {return g(t,y)-c*y;};
+  real[] Y={y};
 
   if(h == 0) {
-    if(b == a) return y;
+    if(b == a) return Y;
     if(n == 0) abort("Either n or h must be specified");
     else h=(b-a)/n;
   }
 
+  real f(real t, real y)=(c == 0 || tableau.exponential) ? g :
+    new real(real t, real y) {return g(t,y)-c*y;};
+
   tableau.stepDependence(h,c,tableau.a);
       
   real t=a;
@@ -165,58 +271,64 @@ real integrate(real y, real c=0, real g(real t, real y), real a, real b=a,
     (tableau.a.lowOrderWeights.length > tableau.a.highOrderWeights.length);
   if(fsal) f0=f(t,y);
 
+  real dt=h;
   while(t < b) {
+    h=min(h,b-t);
+    if(t+h == t) break;
+    if(h != dt) {
+      if(verbose) report(dt,h,t);
+      tableau.stepDependence(h,c,tableau.a);
+      dt=h;
+    }
+ 
     real[] predictions={fsal ? f0 : f(t,y)};
-    for(int i=0; i < tableau.steps.length; ++i)
-      predictions.push(f(t+h*tableau.steps[i],
+    for(int i=0; i < tableau.a.steps.length; ++i)
+      predictions.push(f(t+h*tableau.a.steps[i],
                          tableau.a.factors[i]*y+h*dot(tableau.a.weights[i],
                                                       predictions)));
 
     real highOrder=h*dot(tableau.a.highOrderWeights,predictions);
-    real Y=tableau.a.factors[tableau.steps.length]*y;
+    real y0=tableau.a.factors[tableau.a.steps.length]*y;
     if(dynamic) {
       real f1;
       if(fsal) {
-        f1=f(t+h,Y+highOrder);
+        f1=f(t+h,y0+highOrder);
         predictions.push(f1);
       }
       real lowOrder=h*dot(tableau.a.lowOrderWeights,predictions);
       real error;
-      error=error(error,y,Y+lowOrder,Y+highOrder,highOrder-lowOrder);
-      real dt=h;
-      h=adjust(h,error,t,c,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+      error=error(error,y,y0+lowOrder,y0+highOrder,highOrder-lowOrder);
+      h=adjust(h,error,tolmin,tolmax,tableau);
       if(h >= dt) {
         t += dt;
-        y=Y+highOrder;
+        y=y0+highOrder;
+        Y.push(y);
         f0=f1;
       }
+      h=min(max(h,dtmin),dtmax);
     } else {
       t += h;
-      y=Y+highOrder;
-    }
-    real remain=b-t;
-    if(h > remain) {
-      h=remain;
-      tableau.stepDependence(h,c,tableau.a);
+      y=y0+highOrder;
+      Y.push(y);
     }
-    if(t >= b || t+h == t) break;
   }
-  return y;
+  return Y;
 }
 
 // Integrate a set of equations, dy/dt=f(t,y), from a to b using initial
 // conditions y, specifying either the step size h or the number of steps n.
-real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
-                 real h=0, int n=0, bool dynamic=false,
-                 real tolmin=0, real tolmax=0, real dtmin=0, real dtmax=realMax,
-                 RKTableau tableau, bool verbose=false)
+real[][] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
+                   real h=0, int n=0, bool dynamic=false,
+                   real tolmin=0, real tolmax=0, real dtmin=0,
+                   real dtmax=realMax, RKTableau tableau, bool verbose=false)
 {
+  real[][] Y={copy(y)};
+      
   if(h == 0) {
-    if(b == a) return y;
+    if(b == a) return Y;
     if(n == 0) abort("Either n or h must be specified");
     else h=(b-a)/n;
   }
-  real[] y=copy(y);
   real t=a;
   real[] f0;
   if(tableau.a.lowOrderWeights.length == 0) dynamic=false;
@@ -224,10 +336,18 @@ real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
     (tableau.a.lowOrderWeights.length > tableau.a.highOrderWeights.length);
   if(fsal) f0=f(t,y);
 
+  real dt=h;
   while(t < b) {
+    h=min(h,b-t);
+    if(t+h == t) break;
+    if(h != dt) {
+      if(verbose) report(dt,h,t);
+      dt=h;
+    }
+
     real[][] predictions={fsal ? f0 : f(t,y)};
-    for(int i=0; i < tableau.steps.length; ++i)
-      predictions.push(f(t+h*tableau.steps[i],
+    for(int i=0; i < tableau.a.steps.length; ++i)
+      predictions.push(f(t+h*tableau.a.steps[i],
                          y+h*tableau.a.weights[i]*predictions));
 
     real[] highOrder=h*tableau.a.highOrderWeights*predictions;
@@ -242,21 +362,21 @@ real[] integrate(real[] y, real[] f(real t, real[] y), real a, real b=a,
       for(int i=0; i < y.length; ++i)
         error=error(error,y[i],y[i]+lowOrder[i],y[i]+highOrder[i],
                     highOrder[i]-lowOrder[i]);
-      real dt=h;
-      h=adjust(h,error,t,0,tolmin,tolmax,dtmin,min(dtmax,b-t-h),tableau,verbose);
+      h=adjust(h,error,tolmin,tolmax,tableau);
       if(h >= dt) {
         t += dt;
         y += highOrder;
+        Y.push(y);
         f0=f1;
       }
+      h=min(max(h,dtmin),dtmax);
     } else {
       t += h;
       y += highOrder;
+      Y.push(y);
     }
-    h=min(h,b-t);
-    if(t >= b || t+h == t) break;
   }
-  return y;
+  return Y;
 }
 
 real[][] finiteDifferenceJacobian(real[] f(real[]), real[] t,
@@ -292,7 +412,8 @@ real[] solveBVP(real[] f(real, real[]), real a, real b=a, real h=0, int n=0,
                 real[] guess, RKTableau tableau, int iterations=100)
 {
   real[] g(real[] t) {
-    return discrepancy(integrate(initial(t),f,a,b,h,n,tableau));
+    real[][] y=integrate(initial(t),f,a,b,h,n,tableau);
+    return discrepancy(y[y.length-1]);
   }
   real[][] jacobian(real[] t) {return finiteDifferenceJacobian(g,t);}
   return initial(newton(iterations,g,jacobian,guess));