rouben@math15.math.umbc.edu (Rouben Rostamian) (02/07/91)
SPLINE is a program for generating cubic spline interpolants
on the HP48sx calculator. Given an array [y_1,y_2,...,y_n] of
numbers and an interval [a,b], it returns a piecewise-cubic,
twice differentiable function defined on the interval which takes
on the values y_i at the "node points". The node points are n
equally spaced points on the interval [a,b]. The independent
variable defaults to "X" unless it is specified otherwise in level 1.
SPLINE does not use, generate, modify, or delete any global variables.
It does not disturb the stack, nor does it set or modify any flags,
although the calculator should be in the symbolic mode for things to work.
Example 1: Enter
2: [ 0 1 4 9 ] # The array [y_1,y_2,...,y_n]
1: [0 3] # The interval [a,b]
and press SPLINE to get:
1: '(-1/10*BSPL(1+X)) + 1/10*BSPL(-1+X)...
+ 3/5*BSPL(-2+X) + 3/2*BSPL(-3+X) + 12/5*BSPL(-4+X)'
The graph of this function passes through the coordinate points
(0,0), (1,1), (2,4), (3,9). Plot it and compare with the graph of
x^2 on the interval [0,3].
Example 2: Enter
3: [ 0 1 4 9 ]
2: [0 3]
1: 'T2'
and press SPLINE to get the same function but with the independent
variable "T2":
1: '(-1/10*BSPL(1+T2)) + 1/10*BSPL(-1+T2)...
+ 3/5*BSPL(-2+T2) + 3/2*BSPL(-3+T2) + 12/5*BSPL(-4+T2)'
The function BSPL (and its derivative derBSPL) are included in the
program below. Although the interpolating function is defined for
all values of of its argument, its values outside the specified
interval [a,b] have little meaning.
[ Note: The result shown above were obtained after applying "6 FIX"
and "->Q"; normally the coefficients are presented as decimal expansions.]
It is possible to differentiate the interpolating function. Therefore
the f'(x) and EXTR routines in the graphic menu work just like they do
with the built-in functions.
Example 3: Let's look at the derivative of the interpolant:
3: '(-1/10*BSPL(1+X)) + 1/10*BSPL(-1+X)...
+ 3/5*BSPL(-2+X) + 3/2*BSPL(-3+X) + 12/5*BSPL(-4+X)'
2: 'X'
1: (insert the derivative symbol here)
the result is (sorry about this :-) ):
1: '-(1/10*(((SIGN(ABS(1+X))-SIGN(ABS(1+X)-1))*(9*(1+X)^2-12*ABS(1+X))
+(SIGN(ABS(1+X)-1)-SIGN(ABS(1+X)-2))*-(3*(2-ABS(1+X))^2))*SIGN(1+X)/2))
+1/10*(((SIGN(ABS(-1+X))-SIGN(ABS(-1+X)-1))*(9*(-1+X)^2-12*ABS(-1+X))
+(SIGN(ABS(-1+X)-1)-SIGN(ABS(-1+X)-2))*-(3*(2-ABS(-1+X))^2))*SIGN(-1+X)/2)
+3/5*(((SIGN(ABS(-2+X))-SIGN(ABS(-2+X)-1))*(9*(-2+X)^2-12*ABS(-2+X))
+(SIGN(ABS(-2+X)-1)-SIGN(ABS(-2+X)-2))*-(3*(2-ABS(-2+X))^2))*SIGN(-2+X)/2)
+3/2*(((SIGN(ABS(-3+X))-SIGN(ABS(-3+X)-1))*(9*(-3+X)^2-12* ABS(-3+X))
+(SIGN(ABS(-3+X)-1)-SIGN(ABS(-3+X)-2))*-(3*(2-ABS(-3+X))^2))*SIGN(-3+X)/2)
+12/5*(((SIGN(ABS(-4+X))-SIGN(ABS(-4+X)-1))*(9*(-4+X)^2-12*ABS(-4+X))
+(SIGN(ABS(-4+X)-1)-SIGN(ABS(-4+X)-2))*-(3*(2-ABS(-4+X))^2))*SIGN(-4+X)/2)'
--
Rouben Rostamian Telephone: (301) 455-2458
Department of Mathematics and Statistics e-mail:
University of Maryland Baltimore County bitnet: rostamian@umbc.bitnet
Baltimore, MD 21228, U.S.A. internet: rouben@math9.math.umbc.edu
==================== CUT HERE CUT HERE CUT HERE =====================
%%HP: T(3)A(D)F(.);
DIR
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@ SPLINE is the main routine:
@ Checksums:
@ #24469d
@ 388
SPLINE
\<< DUP
IF TYPE 6 \=/ @ If no independent variable specified on
THEN 'X' @ level 1, then use 'X'.
END @
@ Massage the input stack into a digestable form:
3 ROLLD OBJ\-> DROP ROT DUP SIZE OBJ\-> DROP SWAP
OBJ\-> OBJ\-> DROP \->LIST { } { } 0
@ Save it local variables:
\-> X a b n y x c h
\<<
@ Define the interval size h:
b a - n 1 - / 'h' STO
@ Generate the x_i nodes, i=-1,...,(n+1):
-1 n FOR k
a h k * +
NEXT
n 2 + \->LIST 'x' STO n
@ Generate the coefficient matrix and invert:
CMAT INV
@ Append and prepend 0 to the list of y values.
@ Convert the list 'y' to array 'y':
{ 0 } y + { 0 } + OBJ\-> \->ARRY
@ Compute the coefficients:
* 'c' STO
@ Multiply the B-splines by coefficients and add:
0
1 n 2 + FOR k
'c' k GET 'x' k GET NEG X + 1 \->LIST 'BSPL' APPLY * +
NEXT
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@ Definition of the function "Basic Spline" (B-spline)
@
@ BSPL(x) = 4-6x^+3x^3 if 0<x<1
@ (2-x)^3 if 1<x<2
@ 0 if 2<x
@ and BSPL(-x) = BSPL(x) for all x.
@
@ Checksums:
@ #55554d
@ 154.5
BSPL
\<< \-> x
\<< x ABS 'x' STO
CASE x 2 \>=
THEN 0
END x 1 <
THEN '4-6*x^2+3*x^3' EVAL
END 2 x - 3 ^
END
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@ Definition of the derivative of BSPL:
@
@ Checksums:
@ #35435d
@ 201
derBSPL
\<< \-> x dx
\<< '(SIGN(ABS(x))-SIGN(ABS(x)-1))*(9*x^2-12*ABS(x))' EVAL
'(SIGN(ABS(x)- 1)-SIGN(ABS(x)-2))* (-3*(2-ABS(x))^2)' EVAL
+ x SIGN * 2 / dx *
\>>
\>>
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@@@@@@@@@@@@@@@@@@@@@@@@@@@@
@ Definition of the coefficient matrix. For any n, CMAT generates
@ an (n+2)x(n+2) matrix of the form:
@ [ 6 -12 6
@ 1 4 1
@ 1 4 1
@ 1 4 1
@ ...
@ 1 4 1
@ 6 -12 6 ]
@
@ Checksums:
@ #49706d
@ 217.5
CMAT
\<< 2 + \-> n
\<< n IDN 2 *
2 n FOR i
i i 1 - 2 \->LIST 1 PUT
NEXT
DUP TRN +
{ 1 1 } 6 PUT
{ 1 2 } -12 PUT { 1 3 } 6 PUT
n DUP 2 - 2 \->LIST 6 PUT
n DUP 1 - 2 \->LIST -12 PUT
n DUP 2 \->LIST 6 PUT
\>>
\>>
ENDrouben@math13.math.umbc.edu (Rouben Rostamian) (02/07/91)
Darn! The spline program that I posted to this newsgroup earlier today was an older version and it does not work as advertised. Goes to show that one should not do these things when in a hurry :-( I will post the correct version tomorrow. My apologies to the net. -- Rouben Rostamian Telephone: (301) 455-2458 Department of Mathematics and Statistics e-mail: University of Maryland Baltimore County bitnet: rostamian@umbc.bitnet Baltimore, MD 21228, U.S.A. internet: rouben@math9.math.umbc.edu
rouben@math15.math.umbc.edu (Rouben Rostamian) (02/08/91)
I mentioned in a previous message that the program SPLINE that I had
posted yesterday was an older version and did not perform as advertised.
Here is the latest version. Please discard the old version.
---
SPLINE is a program for generating cubic spline interpolants
on the HP48sx calculator. Given an array [y_1,y_2,...,y_n] of
numbers and an interval [a,b], it returns a piecewise-cubic,
twice differentiable function defined on the interval which takes
on the values y_i at the "node points". The node points are n
equally spaced points on the interval [a,b]. The independent
variable defaults to "X" unless it is specified otherwise in level 1.
SPLINE does not use, generate, modify, or delete any global variables.
It does not disturb the stack, nor does it set or modify any flags,
although the calculator should be in the symbolic mode for things to work.
Example 1: Enter
2: [ 0 1 4 9 ] # The array [y_1,y_2,...,y_n]
1: [0 3] # The interval [a,b]
and press SPLINE to get:
1: '(-1/10*BSPL(1+X)) +1/10*BSPL(-1+X)...
+ 3/5*BSPL(-2+X) +3/2*BSPL(-3+X) +12/5*BSPL(-4+X)'
The graph of this function passes through the coordinate points
(0,0), (1,1), (2,4), (3,9). Plot it and compare with the graph of
x^2 on the interval [0,3].
Example 2: Enter
3: [ 0 1 4 9 ]
2: [0 3]
1: 'T2'
and press SPLINE to get the same function but with the independent
variable "T2":
1: '(-1/10*BSPL(1+T2)) + 1/10*BSPL(-1+T2)...
+3/5*BSPL(-2+T2) +3/2*BSPL(-3+T2) +12/5*BSPL(-4+T2)'
Example 3: Enter
2: [ 0 1 4 9 ] # The array [y_1,y_2,...,y_n]
1: [0 1] # The interval [a,b]
and press SPLINE to get:
1: '-(1/10*BSPL((1/3+X)*3)) +1/10*BSPL((-(1/3)+X)*3) ...
+3/5*BSPL((-(2/3)+X)*3) +3/2*BSPL((-1+X)*3) ...
+12/5*BSPL((-(4/3)+X)*3)'
The graph now passes through (0,0), (1/3,1), (2/3,4), (1,9). Compare
with the graph of 9x^2.
The function BSPL (and its derivative derBSPL) are included in the
program below. Although the interpolating function is defined for
all values of of its argument, its values outside the specified
interval [a,b] have little meaning.
[ Note: The result shown above were obtained after applying "6 FIX"
and "->Q"; normally the coefficients are presented as decimal expansions.]
It is possible to differentiate the interpolating function. Therefore
the f'(x) and EXTR routines in the graphic menu work just like they do
with the built-in functions.
Example 3: Let's look at the derivative of the interpolant:
3: '(-1/10*BSPL(1+X)) + 1/10*BSPL(-1+X)...
+ 3/5*BSPL(-2+X) + 3/2*BSPL(-3+X) + 12/5*BSPL(-4+X)'
2: 'X'
1: (insert the derivative symbol here)
the result is (sorry about this :-) ):
1: '-(1/10*(((SIGN(ABS(1+X))-SIGN(ABS(1+X)-1))*(9*(1+X)^2-12*ABS(1+X))
+(SIGN(ABS(1+X)-1)-SIGN(ABS(1+X)-2))*-(3*(2-ABS(1+X))^2))*SIGN(1+X)/2))
+1/10*(((SIGN(ABS(-1+X))-SIGN(ABS(-1+X)-1))*(9*(-1+X)^2-12*ABS(-1+X))
+(SIGN(ABS(-1+X)-1)-SIGN(ABS(-1+X)-2))*-(3*(2-ABS(-1+X))^2))*SIGN(-1+X)/2)
+3/5*(((SIGN(ABS(-2+X))-SIGN(ABS(-2+X)-1))*(9*(-2+X)^2-12*ABS(-2+X))
+(SIGN(ABS(-2+X)-1)-SIGN(ABS(-2+X)-2))*-(3*(2-ABS(-2+X))^2))*SIGN(-2+X)/2)
+3/2*(((SIGN(ABS(-3+X))-SIGN(ABS(-3+X)-1))*(9*(-3+X)^2-12* ABS(-3+X))
+(SIGN(ABS(-3+X)-1)-SIGN(ABS(-3+X)-2))*-(3*(2-ABS(-3+X))^2))*SIGN(-3+X)/2)
+12/5*(((SIGN(ABS(-4+X))-SIGN(ABS(-4+X)-1))*(9*(-4+X)^2-12*ABS(-4+X))
+(SIGN(ABS(-4+X)-1)-SIGN(ABS(-4+X)-2))*-(3*(2-ABS(-4+X))^2))*SIGN(-4+X)/2)'
--
Rouben Rostamian Telephone: (301) 455-2458
Department of Mathematics and Statistics e-mail:
University of Maryland Baltimore County bitnet: rostamian@umbc.bitnet
Baltimore, MD 21228, U.S.A. internet: rouben@math9.math.umbc.edu
==================== CUT HERE CUT HERE CUT HERE =====================
%%HP: T(3)A(R)F(.);
DIR
SPLINE
@ Checksum: #28772d
@ 424.5
\<< DUP
IF TYPE 6 \=/
THEN 'X'
END
3 ROLLD OBJ\-> DROP ROT DUP SIZE OBJ\-> DROP SWAP
OBJ\-> OBJ\-> DROP \->LIST { } { } 0 0
\-> X a b n y x c h h1
\<< b a - n 1 - DUP2 / 'h' STO SWAP / 'h1' STO
-1 n FOR k a h k * +
NEXT
n 2 + \->LIST 'x' STO n
CMAT INV { 0 } y + { 0 } + OBJ\-> \->ARRY * 'c' STO
0
1 n 2 +
FOR k 'c' k GET 'x' k GET NEG X + h1 * 1 \->LIST 'BSPL' APPLY * +
NEXT
\>>
\>>
BSPL
@ Checksum: #55554d
@ 154.5
\<< \-> x
\<< x ABS 'x' STO
CASE x 2 \>=
THEN 0
END x 1 <
THEN '4-6 *x^2+3*x^3' EVAL
END 2 x - 3 ^
END
\>>
\>>
derBSPL
@ Checksum: #35435d
@ 201
\<< \-> x dx
\<< '(SIGN(ABS(x))-SIGN(ABS(x)-1))*(9*x^2-12*ABS(x))' EVAL
'(SIGN(ABS(x)-1)-SIGN(ABS(x)-2))*(-3*(2-ABS(x))^2)' EVAL
+ x SIGN * 2 / dx *
\>>
\>>
CMAT
@ Checksum: #49706d
@ 217.5
\<< 2 + \-> n
\<< n IDN 2 *
2 n FOR i
i i 1 - 2 \->LIST 1 PUT
NEXT
DUP TRN +
{ 1 1 } 6 PUT
{ 1 2 } -12 PUT
{ 1 3 } 6 PUT
n DUP 2 - 2 \->LIST 6 PUT
n DUP 1 - 2 \->LIST -12 PUT
n DUP 2 \->LIST 6 PUT
\>>
\>>
END