In an engineering project, I have been assigned to a task constructing smooth curve for given sets of sample points. This is accomplished by implementing Cubic Splines and Thomas algorithm to solve linear algebraic equation. So, I write PLSQL stored procedure as follows:

IN x_arr numeric[], IN y_arr numeric[], IN resolution numeric,

OUT xs numeric, OUT ys numeric)

RETURNS SETOF record AS

$BODY$

DECLARE

x numeric[];

y numeric[];

e numeric[] DEFAULT ‘{}’::numeric[];

f numeric[] DEFAULT ‘{}’::numeric[];

g numeric[] DEFAULT ‘{}’::numeric[];

r numeric[] DEFAULT ‘{}’::numeric[];

f2 numeric[] DEFAULT ‘{}’::numeric[];

count integer;

i integer;

j integer;

idx integer;

tmp numeric;

BEGIN

x:=x_arr;

y:=y_arr;

count:=ARRAY_UPPER(x,1);

IF count<3 THEN

RAISE EXCEPTION ‘sample size must be >=3’;

END IF;

IF ARRAY_UPPER(y,1)<>count THEN

RAISE EXCEPTION ‘x and y array size must be identical’;

END IF;

——–sort x, ascendingly:

FOR i IN 1..count-1 LOOP

idx:=i+1;

FOR j IN i+2..count LOOP

IF x[j]<x[idx] THEN

idx:=j;

END IF;

END LOOP;

IF x[idx]=x[i] THEN

RAISE EXCEPTION ‘x values must be unique’;

ELSIF x[idx]<x[i] THEN –swap

tmp=x[i];

x[i]=x[idx];

x[idx]=tmp;

tmp=y[i];

y[i]=y[idx];

y[idx]=tmp;

END IF;

END LOOP;

——–prepare variables for Thomas algorithm

FOR i IN 1..count LOOP

IF i=1 OR i=count THEN

e:= e || 0.0;

f:= f || 0.0;

g:= g || 0.0;

r:= r || 0.0;

f2:=f2 || 0.0;

CONTINUE;

END IF;

e:= e || x[i] – x[i-1];

f:= f || 2*(x[i+1]-x[i-1]);

g:= g || x[i+1] – x[i];

r:= r || ( (6/(x[i+1]-x[i])) * (y[i+1]-y[i]) )

+ ( (6/(x[i]-x[i-1])) * (y[i-1]-y[i]) );

f2:=f2 || 0.0;

END LOOP;

——–decomposition

FOR i IN 2..count LOOP

IF f[i-1]<>0 THEN

e[i]:=e[i]/f[i-1];

END IF;

f[i]:=f[i]-(e[i]*g[i-1]);

END LOOP;

——–forward substitution

FOR i IN 2..count LOOP

r[i]:=r[i]-(e[i]*r[i-1]);

END LOOP;

——–backward substitution

FOR i IN REVERSE count-1..1 LOOP

IF f[i]<>0 THEN

f2[i]:=(r[i] – (g[i]*f2[i+1])) / f[i];

END IF;

END LOOP;

——–do interpolation

FOR i IN 2..count LOOP

xs:=x[i-1]-resolution;

LOOP

xs:=xs+resolution;

EXIT WHEN i=count AND xs>x[i];

EXIT WHEN i=x[i];

ys:= ( (f2[i-1] * (x[i]-xs)^3)

+ (f2[i] * (xs-x[i-1])^3)

) / (6*(x[i]-x[i-1]));

ys:= ys + (

(y[i-1]/(x[i]-x[i-1])) –

(f2[i-1]*(x[i]-x[i-1])/6)

) * (x[i]-xs);

ys:= ys + (

(y[i]/(x[i]-x[i-1])) –

(f2[i]*(x[i]-x[i-1])/6)

) * (xs-x[i-1]);

RETURN NEXT;

EXIT WHEN xs>=x[i];

END LOOP;

END LOOP;

RETURN;

END;

$BODY$

LANGUAGE ‘plpgsql’ IMMUTABLE STRICT;

The stored procedure takes 3 input arguments and 2 outs. **x_arr** and ** y_arr** are sample points’ **x** and **y**. With **resolution**, you can determine the resolution of smoothed curve’s **x**.

Output arguments: **xs** and **ys** contain the set of **x** and **y** of smoothed curve, respectively.

Let us use the function. At first, we have sample points:

# | x | y |
---|---|---|

1 | 3.3 | 1.2 |

2 | 4.2 | -1.1 |

3 | 6.8 | 1.5 |

4 | 8.2 | 1.2 |

5 | 10.2 | 2.5 |

Call the function, supply x and y array as displayed on above table, set resolution to **0.1**:

‘{3.3,4.2,6.8,8.2,10.2}’::numeric[],

‘{1.2,-1.1,1.5,1.2,2.5}’::numeric[],

0.1);

We will have set of records containing **xs** and **ys**.:

3.3 1.199999999999999970

3.4 0.8860175045157231811

3.5 0.5764170295261004851

3.6 0.2755805955257859750

3.7 -0.0121097769905662565

3.8 -0.28227206752830211620

… and so on

Let us put them in graph:

The blue-y1 and red-y2 are sample and interpolation points as returned by the **spline** function, respectively.