Mathcad 8 Pro. code conversion

2c4ad377-1966-4f15-87f5-4ee181e3216b 57 4 5 decimal

&[DATE] &[TIME] - &[FILENAME]

The 3 types of "cubic spline"

Essentially, the coefficients of the cubic spline consists in solving the basis tridiagonal matrix wrt the 'br' vector of derivatives . Independently of the type of cubic spline [type: 0, 1, 2 = l, p, c] two steps of preparation are needed: 1. Install the constant vector 'b' 2. Install the blank tridiagonal matrix.

XY 3.52 - 2.358 - 3 - 2.213 - 2 - 1.256 - 1 - 1.739 - 0 1.482 - 1 1.277 - 62 mat :

X XY 1 col :

Y XY 2 col :

N Y rows 1 - :

N 5

r 2 N range :

Prep.0 Install the constant vector 'br'

i r br i el Y i 1 + el Y i el - X i 1 + el X i el - / (Y i el Y i 1 - el - X i el X i 1 - el - / (- : br N 1 + el 0 : 21 line for

br 0 0.6782 1.44 - 0.74 0.052 - 061 mat

Prep.1 Install the tridiagonal matrix

M 11 el 1 :

M N 1 + N 1 + el 1 :

M 10000000000000000000000000000000000166 mat

Observe the population of the rows upper left row ---------------------- ----------------------lower right row respectively to linear/parabolic/cubic. The basis tridiagonal is the "linear", further modified to suit user l_p_c .

0. Construct the tridiagonal linear spline matrix [lspline]

i r M i i 1 - el X i el X i 1 - el - 6 / : for

i r M i i el X i 1 + el X i 1 - el - 3 / : for

i r M i i 1 + el X i 1 + el X i el - 6 / : for

s M 1 -^br * :

M 100000 0.0867 0.5067 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 00000166 mat

s 0 2.419 3.285 - 2.0808 0.5982 - 061 mat

1. Construct the tridiagonal parabolic spline matrix [pspline]

M 12 el 1 - :

M N 1 + N el 1 - :

s M 1 -^br * :

M 11 - 0000 0.0867 0.5067 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 00001 - 166 mat

s 2.0344 2.0344 3.1736 - 2.02 0.4664 - 0.4664 - 61 mat

2. Construct the tridiagonal cubic spline matrix [cspline]

M 12 el X 1 el X 3 el - (X 3 el X 2 el - / (:

M N 1 + N el X N el X N 2 - el - (X N 2 - el X N 1 - el - / :

M 13 el X 2 el X 1 el - (X 3 el X 2 el - (/ :

M N 1 + N 1 - el X N 1 - el X N el - (X N 1 - el X N el - (/ :

s M 1 -^br * :

M 1 1.52 - 0.52 000 0.0867 0.5067 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 0000 0.1667 0.6667 0.1667 00012 - 166 mat

s 7.1661 1.0623 2.8867 - 1.8447 0.052 - 1.9487 - 61 mat

Construct the general spline cubic polynomial

s 0 2.419 3.285 - 2.0808 0.5982 - 061 mat :

In there, plug the 's' of the selected type of spline. In this demo, we select "linear" type in order to view all segments on the XY data set

Linear 0 2.419 3.285 - 2.0808 0.5982 - 061 mat

Parabolic 2.0344 2.0344 3.1736 - 2.02 0.4664 - 0.4664 - 61 mat

Cubic 4.0766 1.6415 3.0412 - 1.8833 0.052 - 1.9873 - 61 mat

Collect the cubic spline components for C(j,u) abcd i 1 X rows 1 - range hs i el s i 1 + el s i el - : hy i el Y i 1 + el Y i el - : hx i el X i 1 + el X i el - : 31 line for i 1 X rows 1 - range a i el hs i el 6 hx i el * / : b i el s i el 2 / : c i el hy i el hx i el / hx i el 6 / s i 1 + el 2 s i el * + (* - (: d i el Y i el : 41 line for a b c d augment 31 line :

abcd 0.7753 0 0.0692 2.358 - 0.9507 - 1.2095 0.6982 2.213 - 0.8943 1.6425 - 0.2652 1.256 - 0.4465 - 1.0404 0.3369 - 1.739 - 0.0997 0.2991 - 0.4044 1.482 - 54 mat

Recover segments coeff's [a,b,c,d]

a abcd 1 col :

c abcd 3 col :

b abcd 2 col :

d abcd 4 col :

j u C i j : a i el u X i el - (3^* b i el u X i el - (2^* + c i el u X i el - (* + d i el + 21 line :

Global pts plot

Sy . :

size 14 :

clr black :

ptsSpec k 1 : i k X rows range data i el Sy size clr 13 mat : k k 1 + : 21 line for 21 line Unest P v data 1 el : j 2 data rows range v v data j el stack : for v 31 line 31 line :

Data X Y ptsSpec augment :

plot 1 x C 2 x C 3 x C 4 x C 5 x C Data 61 sys :

2 2.5 - C 1.680375 -

Cubic spline segments through data points

plot

support 3.52 - 2.358 - . 14 black 3 - 2.213 - . 14 black 2 - 1.256 - . 14 black 1 - 1.739 - . 14 black 0 1.482 - . 14 black 1 1.277 - . 14 black 65 mat :

Demo: collect the 'j' segments [if there are few]. u Spline j local assignment n X rows : u X 1 el ≤ n 1 : u X n el ≥ n n 1 - : i 1 n range u X i el ≤ break return 'j' < (i-1) if for if if i 1 - (41 line : j u C 21 line :

x segments x X 1 el ≥ (x X X rows el ≤ (& x Spline Not interpolated if :

"Proofing" the cubic spline project

x segments 4 x C 3.52 - 2.358 - . 14 black 3 - 2.213 - . 14 black 2 - 1.256 - . 14 black 1 - 1.739 - . 14 black 0 1.482 - . 14 black 1 1.277 - . 14 black 6 5 mat 3 1 sys

We want a more general tool for lot more segments

Mathcad 8 Pro. code conversion

s X Y t InterpSpline Mathcad 8 Pro. code conversion t X 1 el ≤ idx 1221 mat : t X N el ≥ idx N N 1 + 21 mat : i 1 : X i el t < i i 1 + : while idx i 1 - i 21 mat : 31 line if if hx X idx 2 el el X idx 1 el el - (: hy Y idx 2 el el Y idx 1 el el - (: hs s idx 2 el el s idx 1 el el - (: 31 line a hs 6 hx * / : b s idx 1 el el 2 / : c hy hx / hx 6 / s idx 2 el el 2 s idx 1 el el * + (* - (: d Y idx 1 el el : 41 line a t X idx 1 el el - (3^* b t X idx 1 el el - (2^* + c t X idx 1 el el - (* + d + 51 line :

x MySpline s X Y x InterpSpline :

x Smath X Y x ainterp :

collect/simplify plot x cubic X 1 el x ≤ (x X X rows el ≤ (& (x MySpline Don't interpolate if :

"Proofing" the cubic spline project

x cubic 3 x C X 1 el Y 1 el . 14 black X 2 el Y 2 el . 14 black X 3 el Y 3 el + 14 red X 4 el Y 4 el . 14 black X 5 el Y 5 el . 14 black X 6 el Y 6 el . 14 black 6 5 mat 3 1 sys

Discretise for export L H N xd U 0 : dx H L - N / : i 1 N 1 + range U i el L dx - (dx i * + : for U 41 line :

L 3.52 - : H 1 : Npts 100 : 31 line

U L H Npts xd :

Ω i 1 U rows range vy i el U i el MySpline : for U vy augment 21 line :

Discrete, user segment, Data

Ω 3 x C 5 x C Data 4 1 sys

ω i 1 U rows range vy i el U i el Smath : for U vy augment 21 line :

Ω Ω :

ω ω :

This plot compares Smath cinterp(X,Y,x) with the Western conventional linear/parabolic/cubic splines. The closest Smath compares is the "linear", but it is not especially in the end points regions. Akima [from ALGLIB [ainterp] concerns the ALGLIB casual case.

Ω ω Data 3 1 sys

5 - MySpline 5 - Smath 21 mat 4.973843 - 1.450077 - 21 mat

linear... 4.973843 - 4.613213 - 21 mat ≡ parabolic 0.24022 4.613213 - 21 mat ≡ cubic.... 5.48049 4.613213 - 21 mat ≡ 31 mat

Visit for Akima: http://www.alglib.net/interpolation/spline3.php

data Deriv data Derivative BackSubstitution X data 1 col : Y data 2 col : n X rows : deriv 1 el 0 : deriv n el 0 : i 2 n 1 - range deriv i el 0.5 Y i 1 - el Y i el - X i 1 - el X i el - / Y i el Y i 1 - el - X i el X i 1 - el - / + (* eval : for deriv 1 el deriv 2 el : deriv n el deriv deriv rows el : X deriv augment 101 line :

Dfinite Ω Deriv :

Finite differences derivative

Ω Deriv

The complimentary scalar dspline s X Y t dspline t X 2 el ≤ k 1 : t X X rows 1 - el ≥ k X rows 1 - : k 1 : t X k 1 + el > k k 1 + : while 21 line if if c Y k 1 + el Y k el - X k 1 + el X k el - / s k 1 + el 2 s k el * + (X k 1 + el X k el - 6 / (* - : s k 1 + el s k el - X k 1 + el X k el - / (t X k el - (2^2 / * s k el t X k el - (* + c + 31 line :

x D1 X 1 el x ≤ (x X 6 el ≤ (& s X Y x dspline Not interpolated if :

The scalar dspline

x D1 0 2 1 sys The while loop "s" interpolation

A shorther form of the spline interpolation function s X Y t Interpspline Find the knots that bracket 't' ct 1 : j X rows : j ct - 1 > i ct j + 2 / Floor : X i el t > j i : ct i : if 21 line while Evaluate the spline based on those knots h X j el X ct el - : a X j el t - h / : b t X ct el - h / : a Y ct el * b Y j el * + h 2^6 / a a 2^1 - (* s ct el * b b 2^1 - (* s j el * + (* + 91 line :

x f1 X 1 el x ≤ (x X 6 el ≤ (& s X Y x Interpspline Not interpolated if :

Compare the two interpolations

x f1 x MySpline 2 1 sys

X Ω 1 col : Y Ω 2 col : 21 line

1. ω Akima => crap 2. Ω lspline => excellent

N Y rows 1 - :

N 100

r 2 N range :

i r br i el Y i 1 + el Y i el - X i 1 + el X i el - / (Y i el Y i 1 - el - X i el X i 1 - el - / (- : br N 1 + el 0 : 21 line for

M M Clear :

M 1

M 11 el 1 :

M N 1 + N 1 + el 1 :

i r M i i 1 - el X i el X i 1 - el - 6 / : for

i r M i i el X i 1 + el X i 1 - el - 3 / : for

i r M i i 1 + el X i 1 + el X i el - 6 / : for

<= Linear

M M :

s M 1 -^br * eval :

2nd Derivative from 's' coefficients

CWD C:\SmathFile CurrentDirectory :

Spool to file the spline coefficcients s CWD splineD2 wfile 1

s CWD splineD2 rfile :

It is true that augment(X,s) represents the 2nd derivative of a splined data set. It was not discussed/revealed in the Mathsoft Collaboratory. It is true under the "linear type" of the 's' coefficients and that Ω is "true", i.e: from a true and clean Ω. In this proof, Ω was obtained from the linear cubic of the 6 original X,Y pairs. 'ω' from Smath cinterp comes close. 'ω' from Akima ainterp is stirctly a piece of crap.

X s augment

Read the suite: Spline l_p_c [dspline, fintegrate].sm

NO reconciliation: Akima=> best or crap !

plotG

vx vy a b c plotG a char ≡ b size ≡ c color ≡ 13 mat n vx length : plot vx 1 el vy 1 el a b c augment : i 2 n range plot plot vx i el vy i el a b c augment stack : for plot 51 line :

XY 0 0.00794 0.1 0.012 0.2 0.022 0.3 0.048 0.4 0.167 0.5 1 0.6 0.167 0.7 0.048 0.8 0.022 0.9 0.012 1 0.00794 112 mat :

X XY 1 col : Y XY 2 col : data X Y + 10 black plotG : 31 line

Resume the "No reconciliation": This symmetric Gaussian is best splined by Akima. Unlike this best case, if our 6 points demo from the top of the work sheet has been "Akima splined" and if from this splined/populated we further want the 2nd "s derivative", the result is just scrap !

L 0 : H 1 : Npts 100 : 31 line

U L H Npts xd :

Linear [s]

N Y rows 1 - :

N 10

r 2 N range :

i r br i el Y i 1 + el Y i el - X i 1 + el X i el - / (Y i el Y i 1 - el - X i el X i 1 - el - / (- : br N 1 + el 0 : 21 line for

M 11 el 1 :

M N 1 + N 1 + el 1 :

s s Clear :

0. Construct the tridiagonal linear spline matrix [lspline]

i r M i i 1 - el X i el X i 1 - el - 6 / : for

i r M i i el X i 1 + el X i 1 - el - 3 / : for

i r M i i 1 + el X i 1 + el X i el - 6 / : for

s M 1 -^br * eval :

s s :

x Akima X Y x ainterp :

x Cinterp X Y x cinterp :

x Linterp s X Y x Interpspline :

Akim i 1 U rows range vy i el U i el Akima eval : for U vy augment 21 line :

Cin i 1 U rows range vy i el U i el Cinterp eval : for U vy augment 21 line :

Lin i 1 U rows range vy i el U i el Linterp eval : for U vy augment 21 line :

From the above comparison, Smath 'cinterp' was the closest match to "conventional cubic linear". This hard case proves.

Akim Cin Lin data 4 1 sys Lin Cin data 3 1 sys