Hi everyone,

For my research, I needed a procedure to calculate an interpolant respecting the monotonicity of the given data. The curve fitting package of Maple 11 didn't help.

I'm pasting my code below.  I hope it helps some of you too.

Cheers,

Ozgur

PS: Thanks goes to Joe Riel for his help.

Cubic := proc(x, l, fl, u, fu, dl, du)
     local H1, H2, H3, H4, hi, pIi;
     option `Copyright (c) 2008 by Ozgur Inal`;     
     hi := u-l;
     #t := (x-l)/hi;
#Hi(x) are the Hermite Basis Functions as given in Fritsch and Carlson '80     
     H1(x) := 3*((u-x)/hi)^2-2*((u-x)/hi)^3;
     H2(x) := 3*((x-l)/hi)^2-2*((x-l)/hi)^3;
     H3(x) := -hi*(((u-x)/hi)^3-((u-x)/hi)^2);
     H4(x) := hi*(((x-l)/hi)^3-((x-l)/hi)^2);
     pIi(x) := fl*H1(x)+fu*H2(x)+dl*H3(x)+du*H4(x);
end proc:
 

MonotoneSpline := proc()
#enter the points to be interpolated:     
local i, k, j, n, a, b, D, d, t, p, S;
option `Copyright (c) 2008 by Ozgur Inal`;
n := nargs/2;
if irem(nargs,2)=1 then
  error "one point is missing!"
else
  #initialize d1:    
  if (args[4]-args[2])/(args[3]-args[1])=0 then
     d[1] := 0
  else
     d[1] := (args[4]-args[2])/(args[3]-args[1])
  end if;
  #initialize dn:
  if (args[nargs]-args[nargs-2])/(args[nargs-1]-args[nargs-3]) = 0 then
     d[n] := 0
  else
     d[n] := (args[nargs]-args[nargs-2])/(args[nargs-1]-args[nargs-3])
  end if;  
#initialize the remaining slope parameters:
  for k from 2 to (n-1) do
     d[k] := (args[2*k]-args[2*k-2])/(2*(args[2*k-1]-args[2*k-3]))+(args[2*k+2]-args[2*k])/(2*(args[2*k+1]-args[2*k-1]));
  end do;
  #check whether the slopes satisfy the necessary and sufficient contidionts:
  for i from 1 to (n-1) do
    D[i] := (args[2*i+2]-args[2*i])/(args[2*i+1]-args[2*i-1]);
      if D[i]=0 then     
        d[i] := 0;
        d[i+1] := 0;
      else
        a[i] := d[i]/D[i];
        b[i] := d[i+1]/D[i];
        if a[i]^2+b[i]^2>9 then
          t[i] := 3*((a[i]^2+b[i]^2)^(-0.5));
          d[i] := t[i]*d[i];
          d[i+1] := t[i]*d[i+1]; #Fritsch and Carlson '80 page 242, right before section 5
         end if;
      end if;
  end do;
  #calculate each cubic polynomial piece:
     for i from 1 to (nargs/2-1) do;
          p[i] := Cubic(x, args[2*i-1], args[2*i], args[2*i+1], args[2*i+2], d[i], d[i+1]);
     end do;
  #combine these polynomials as a piecewise function:
  S(x):=piecewise(seq([x<=args[2*i+1], p[i]][], i=1..n-1)); #thanks for this line goes to Joe Riel of mapleprimes.com
  #print the results:
  print(`slope of the spline at the knots are:`);
  print(seq([d[i]][],i=1..n));
  print(`and the function is:`);
  print(expand(S(x)));
  plot(S(x), x=args[1]..args[nargs-1]); 
end if;
end proc: