Path: utzoo!utgpu!news-server.csri.toronto.edu!cs.utexas.edu!wuarchive!udel!haven!umbc3!math15.math.umbc.edu!rouben From: rouben@math15.math.umbc.edu (Rouben Rostamian) Newsgroups: comp.sys.handhelds Subject: SPLINE for HP48 Message-ID: <4965@umbc3.UMBC.EDU> Date: 7 Feb 91 00:57:18 GMT Sender: newspost@umbc3.UMBC.EDU Reply-To: rouben@math15.math.umbc.edu Organization: Mathematics Department University of Maryland, Baltimore County Lines: 209 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 \<< 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 \>> \>> END