:

## monotone polynomial interpolation

Maple

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:

﻿