Announcement

**Dean Schrage** · 15 Nov 2021, 06:22 AM

Code:

#COMPILE EXE "Cubic Spline.exe"

#DIM ALL

   $S_MSGBOX_HDR = "Cubic Spline"
   %T = -1
   %F = 00

   '  Supports tridiagonal matrix solution

   TYPE UDT_LDU

      '  Lower, diagonal and upper vectors

      L                                                     AS DOUBLE
      D                                                     AS DOUBLE
      U                                                     AS DOUBLE

      '  Solution vector x<k> and x<k+1>
      xk                                                    AS DOUBLE
      xkp1                                                  AS DOUBLE

      '  Right hand side source vector
      b                                                     AS DOUBLE

   END TYPE

   '  Supports tridiagonal matrix solver parameters

   TYPE UDT_LDU_PARM

      N_ITR_MAX                                             AS LONG
      e_ITR_MIN                                             AS DOUBLE

   END TYPE


   '  Supports cubic spline interpolation

   TYPE UDT_CSP

      '  Spline coefficients
      h                                                     AS DOUBLE
      b                                                     AS DOUBLE
      v                                                     AS DOUBLE

      '  Source vector
      u                                                     AS DOUBLE

      '  Solution state
      z                                                     AS DOUBLE

   END TYPE

   '  Supports cubic spline input parameters

   TYPE UDT_CSP_PARM

      '  References the input data

      N_DATA                                                AS LONG
     pt_DATA                                                AS DOUBLE POINTER
     py_DATA                                                AS DOUBLE POINTER

      '  References the interpolated data which forms the spline
      N_INTP                                                AS LONG
     pt_INTP                                                AS DOUBLE POINTER
     py_INTP                                                AS DOUBLE POINTER

      '  Supports tridiagonal coefficient solution

      N_LDU                                                 AS LONG
     pU_LDU                                                 AS UDT_LDU POINTER
      U_LDU_PARM                                            AS UDT_LDU_PARM

   END TYPE


FUNCTION PBMAIN () AS LONG

   '  Call the test function
   '  Note this sets the interpolated data to the clipboard
   CSP_TEST

END FUNCTION


FUNCTION GS_LDU_xkp1   (                                       _
                               N_LDU      AS LONG              , _
                        BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS LONG

   '  Routine computes a Gauss Seidel single step iteration

   DIM i                                                 AS LONG

   DIM Lxk                                               AS DOUBLE
   DIM Uxk                                               AS DOUBLE

   IF       N_LDU = 00                                   OR _
            pU_LDU = 00                                  THEN

      FUNCTION = %F

   ELSE

      FOR i = 00 TO N_LDU - 01

         '  Form the matrix (L+U)[x.k] noting boundaries for lower and upper
         Lxk = IIF(i = 00        , %F, @pU_LDU[i].L * @pU_LDU[i - 01].xkp1)   '  Uses most advance computed [x], Gauss Seidel method
         Uxk = IIF(i = N_LDU - 01, %F, @pU_LDU[i].U * @pU_LDU[i + 01].xk)

         '  Solve [x.kp1]
         IF NOT @pU_LDU[i].D = 0.0#                      THEN

            @pU_LDU[i].xkp1 = 1.0# / @pU_LDU[i].D * (@pU_LDU[i].b - (Lxk + Uxk))

         END IF

      NEXT i

      FUNCTION = %T

   END IF

END FUNCTION


FUNCTION GS_LDU_BFR    (                                       _
                               N_LDU      AS LONG              , _
                        BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS LONG

   '  Routine computes a Gauss Seidel solution vector buffering step

   DIM i                                                 AS LONG

   IF       N_LDU = 00                                   OR _
            pU_LDU = 00                                  THEN

      FUNCTION = %F

   ELSE

      FOR i = 00 TO N_LDU - 01

         @pU_LDU[i].xk = @pU_LDU[i].xkp1

      NEXT i

      FUNCTION = %T

   END IF

END FUNCTION


FUNCTION GS_LDU_e_ITR  (                                       _
                               N_LDU      AS LONG              , _
                        BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS DOUBLE

   '  Routine computes a Gauss Seidel solution vector error norm

   DIM i                                                 AS LONG

   DIM e_ITR                                             AS DOUBLE


   IF       N_LDU = 00                                   OR _
            pU_LDU = 00                                  THEN

      FUNCTION = %F

   ELSE

      FOR i = 00 TO N_LDU - 01

         e_ITR += ABS(@pU_LDU[i].xkp1 - @pU_LDU[i].xk)

      NEXT i

      FUNCTION = e_ITR

   END IF

END FUNCTION


FUNCTION GS_LDU_SOLVE  (                                             _
                               N_LDU      AS LONG                    , _
                        BYVAL pU_LDU      AS UDT_LDU        POINTER  , _
                        BYVAL pU_LDU_PARM AS UDT_LDU_PARM   POINTER     )  AS LONG

   '  Routine computes a Gauss Seidel solution to the GS LDU parameters

   DIM k_ITR                                             AS LONG
   DIM e_ITR                                             AS DOUBLE


   IF       N_LDU = 00                                   OR _
            pU_LDU = 00                                  OR _
            pU_LDU_PARM = 00                             THEN

      FUNCTION = %F

   ELSE

      DO

         INCR k_ITR
         GS_LDU_xkp1    N_LDU, pU_LDU              '  Comute single step
         GS_LDU_e_ITR   N_LDU, pU_LDU TO e_ITR     '  Compute error norm
         GS_LDU_BFR     N_LDU, pU_LDU              '  Back buffer [x.k] = [x.kp1]

         #IF %F
         MSGBOX "k_ITR, e_ITR " + STR$(k_itr) + " " + STR$(e_itr) + $CRLF + GS_LDU_GET_S(N_LDU, pU_LDU), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
         #ENDIF

      LOOP WHILE (k_ITR < @pU_LDU_PARM.N_ITR_MAX) AND (e_ITR > @pU_LDU_PARM.e_ITR_MIN)

      FUNCTION = %T

   END IF

END FUNCTION


FUNCTION GS_LDU_GET_S  (                                       _
                               N_LDU      AS LONG              , _
                        BYVAL pU_LDU      AS UDT_LDU POINTER      )  AS STRING

   DIM i                                                 AS LONG
   DIM S_LDU                                             AS STRING


   IF       N_LDU = 00                                   OR _
            pU_LDU = 00                                  THEN

      FUNCTION = ""

   ELSE

      FOR i = 00 TO N_LDU - 01

         S_LDU   += TRIM$("x.kp1[" + STR$(i+1) + "]") + " " + STR$(@pU_LDU[i].xkp1) + $CRLF

      NEXT i

      FUNCTION = S_LDU

   END IF

END FUNCTION


FUNCTION GS_LDU_TEST   (                                          )  AS LONG

   '  Routine computes an Gauss Seidel test function

   DIM N_LDU                                             AS LONG
   DIM U_LDU()                                           AS UDT_LDU
   DIM U_LDU_PARM                                        AS UDT_LDU_PARM


   N_LDU = 3
   REDIM U_LDU(00 TO N_LDU - 01)

   U_LDU(00).L = %F
   U_LDU(01).L = 1.0#
   U_LDU(02).L = 1.0#

   U_LDU(00).D = 8.0#
   U_LDU(01).D = -7.0#
   U_LDU(02).D = 9.0#

   U_LDU(00).U = 1.0#
   U_LDU(01).U = 2.0#
   U_LDU(02).U = %F

   U_LDU(00).b = 8.0#
   U_LDU(01).b = -4.0#
   U_LDU(02).b = 12.0#

   U_LDU(00).xk = 0.0#
   U_LDU(01).xk = 0.0#
   U_LDU(02).xk = 0.0#

   U_LDU_PARM.N_ITR_MAX = 10
   U_LDU_PARM.e_ITR_MIN = 1.0E-6#

  'GS_LDU_xkp1    N_LDU, VARPTR(U_LDU(00))
   GS_LDU_SOLVE   N_LDU, VARPTR(U_LDU(00)), VARPTR(U_LDU_PARM)

   MSGBOX         GS_LDU_GET_S(N_LDU, VARPTR(U_LDU(00))), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR

END FUNCTION


FUNCTION CSP_FIND_INDEX(                                             _
                               N_CSP      AS LONG                    , _
                        BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                        BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER  , _
                               t_INTP     AS DOUBLE                     )  AS LONG

   '  Routine finds the ith index that an interpolated data point would lie between
   '  Note this ranges 00 to n - 01 - 01, thus returns -1 on fail which is tested

   DIM i                                                 AS LONG


   IF       N_CSP = 00                                   OR _
            pU_CSP = 00                                  OR _
            pU_CSP_PARM = 00                             THEN

      FUNCTION = -1

   ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
           @pU_CSP_PARM.pt_DATA = 00                     OR _
           @pU_CSP_PARM.py_DATA = 00                     OR _
           @pU_CSP_PARM.N_INTP = 00                      OR _
           @pU_CSP_PARM.pt_INTP = 00                     OR _
           @pU_CSP_PARM.py_INTP = 00                     OR _
           @pU_CSP_PARM.N_LDU = 00                       OR _
           @pU_CSP_PARM.pU_LDU = 00                      THEN


   ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size

      FUNCTION = -1

   ELSE

      '  Trap the boundaries

      IF       t_INTP <= @[email protected]_DATA[00]                THEN

         FUNCTION = 00

      ELSEIF   t_INTP >= @[email protected]_DATA[N_CSP - 01 - 01]   THEN

         FUNCTION = N_CSP - 01 - 01

      ELSE

         FOR i = 00 TO N_CSP - 01 - 01

            IF t_INTP >= @[email protected]_DATA[i]                 AND _
               t_INTP <  @[email protected]_DATA[i + 01]            THEN

               FUNCTION = i
               EXIT FUNCTION

            END IF

         NEXT i

         '  Fail condition is encountered
         FUNCTION = -1

      END IF

   END IF

END FUNCTION


FUNCTION CSP_SOLVE     (                                             _
                               N_CSP      AS LONG                    , _
                        BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                        BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG

   '  Routine computes a cubic spline solution to a set of input data by pointer

   DIM i                                                 AS LONG
   DIM i_INTP                                            AS LONG


   IF       N_CSP = 00                                   OR _
            pU_CSP = 00                                  OR _
            pU_CSP_PARM = 00                             THEN

      FUNCTION = %F

   ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
           @pU_CSP_PARM.pt_DATA = 00                     OR _
           @pU_CSP_PARM.py_DATA = 00                     OR _
           @pU_CSP_PARM.N_INTP = 00                      OR _
           @pU_CSP_PARM.pt_INTP = 00                     OR _
           @pU_CSP_PARM.py_INTP = 00                     OR _
           @pU_CSP_PARM.N_LDU = 00                       OR _
           @pU_CSP_PARM.pU_LDU = 00                      THEN


   ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size

      FUNCTION = %F

   ELSE

      '  Compute b and h coefficients

      FOR i = 00 TO N_CSP - 01 - 01

         @pU_CSP[i].h = @[email protected]_DATA[i + 01] - @[email protected]_DATA[i]

         IF NOT @pU_CSP[i].h = 0.0#                      THEN

            @pU_CSP[i].b = 1.0# / @pU_CSP[i].h * (@[email protected]_DATA[i + 01] - @[email protected]_DATA[i])

         END IF

      NEXT i


      '  Compute v and u coefficients

      FOR i = 01 TO N_CSP - 01 - 01

         @pU_CSP[i].v = 2.0# * (@pU_CSP[i - 01].h + @pU_CSP[i].h)
         @pU_CSP[i].u = 6.0# * (@pU_CSP[i].b - @pU_CSP[i - 01].b)

      NEXT i


      '  Load the LDU coefficients to prepare for the tridiagonal solution

      FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01

         '  Load the LDU vectors
         '  Note the IIF operator for the boundary condition (i = 00 OR i = N_CSP - 01) which processes S''(x) = 0 by setting L, U = 0 and D = 1
         @[email protected]_LDU[i].L = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = 01             , %F, @pU_CSP[i - 01].h))
         @[email protected]_LDU[i].D = IIF(i = 00 OR i = N_CSP - 01, 1.0#, @pU_CSP[i].v                                   )
         @[email protected]_LDU[i].U = IIF(i = 00 OR i = N_CSP - 01, 0.0#, IIF(i = N_CSP - 01 - 01, %F, @pU_CSP[i].h)     )

         '  Load the b solution source vector which in this case is the [u] vector
         @[email protected]_LDU[i].b = @pU_CSP[i].u

      NEXT i


      '  Solve the cubic spline [z] vector for S.i''(x) using the tridiagonal Gauss Seidel LDU solver

      GS_LDU_SOLVE   @pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU, VARPTR(@pU_CSP_PARM.U_LDU_PARM)
      #IF %F
      MSGBOX GS_LDU_GET_S(@pU_CSP_PARM.N_LDU, @pU_CSP_PARM.pU_LDU), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
      #ENDIF

      '  Range over the span of the interpolated t points

      FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01

         '  Create the uniform interpolation points
         @[email protected]_INTP[i_INTP] = @[email protected]_DATA[00] + (@[email protected]_DATA[@pU_CSP_PARM.N_DATA - 01] - @[email protected]_DATA[00]) * i_INTP / (@pU_CSP_PARM.N_INTP - 01)

         '  Find the data index that each point lies in
         CSP_FIND_INDEX    N_CSP                            , _
                           pU_CSP                           , _
                           pU_CSP_PARM                      , _
                           @[email protected]_INTP[i_INTP]    TO i

         '  Error check the index for -1 from the find index function
         IF NOT i = -1                                   THEN

            '  Compute the cubic spline interpolation y = ao + a1 t + a2 t^2 + a3 t^3
            '  Note the LDU solved [x]kp1 term replaces the [z] solution vector by variable name
            @[email protected]_INTP[i_INTP] = @[email protected]_LDU[i + 01].xkp1 / (6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i]) ^ 3                                                  + _
                                            @[email protected]_LDU[i].xkp1      / (6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP]) ^ 3                                             + _
                                            (@[email protected]_DATA[i + 01] / @pU_CSP[i].h - @[email protected]_LDU[i + 01].xkp1 / 6.0# * @pU_CSP[i].h) * (@[email protected]_INTP[i_INTP] - @[email protected]_DATA[i])       + _
                                            (@[email protected]_DATA[i]      / @pU_CSP[i].h - @[email protected]_LDU[i].xkp1      / 6.0# * @pU_CSP[i].h) * (@[email protected]_DATA[i + 01] - @[email protected]_INTP[i_INTP])

         END IF

      NEXT i_INTP

      FUNCTION = %T

   END IF

END FUNCTION


FUNCTION CSP_SET_CLIP  (                                             _
                               N_CSP      AS LONG                    , _
                        BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                        BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS LONG

   '  Routine sets cubic spline solution to the clipboard for easy plotting

   DIM i                                                 AS LONG
   DIM i_INTP                                            AS LONG

   DIM S_CLIP                                            AS STRING


   IF       N_CSP = 00                                   OR _
            pU_CSP = 00                                  OR _
            pU_CSP_PARM = 00                             THEN

      FUNCTION = %F

   ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
           @pU_CSP_PARM.pt_DATA = 00                     OR _
           @pU_CSP_PARM.py_DATA = 00                     OR _
           @pU_CSP_PARM.N_INTP = 00                      OR _
           @pU_CSP_PARM.pt_INTP = 00                     OR _
           @pU_CSP_PARM.py_INTP = 00                     OR _
           @pU_CSP_PARM.N_LDU = 00                       OR _
           @pU_CSP_PARM.pU_LDU = 00                      THEN


   ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size

      FUNCTION = %F

   ELSE

      '  Range over the span of the interpolated t points

      FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01

         S_CLIP  += STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF

      NEXT i_INTP

      CLIPBOARD SET TEXT S_CLIP

      FUNCTION = %T

   END IF

END FUNCTION


FUNCTION CSP_GET_S     (                                             _
                               N_CSP      AS LONG                    , _
                        BYVAL pU_CSP      AS UDT_CSP        POINTER  , _
                        BYVAL pU_CSP_PARM AS UDT_CSP_PARM   POINTER     )  AS STRING

   DIM i                                                 AS LONG
   DIM i_INTP                                            AS LONG
   DIM S_CSP                                             AS STRING


   IF       N_CSP = 00                                   OR _
            pU_CSP = 00                                  OR _
            pU_CSP_PARM = 00                             THEN

      FUNCTION = ""

   ELSEIF  @pU_CSP_PARM.N_DATA = 00                      OR _     '  Test all allocations and pointers
           @pU_CSP_PARM.pt_DATA = 00                     OR _
           @pU_CSP_PARM.py_DATA = 00                     OR _
           @pU_CSP_PARM.N_INTP = 00                      OR _
           @pU_CSP_PARM.pt_INTP = 00                     OR _
           @pU_CSP_PARM.py_INTP = 00                     OR _
           @pU_CSP_PARM.N_LDU = 00                       OR _
           @pU_CSP_PARM.pU_LDU = 00                      THEN


   ELSEIF  NOT @pU_CSP_PARM.N_DATA = @pU_CSP_PARM.N_LDU  THEN     '  LDU allocation must match the input data size

      FUNCTION = ""

   ELSE

      #IF 00
      FOR i = 00 TO N_CSP - 01

         S_CSP   += TRIM$("b[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].b) + " " + _
                    TRIM$("v[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].v) + " " + _
                    TRIM$("u[" + STR$(i+1) + "]") + " " + STR$(@pU_CSP[i].u) + $CRLF

      NEXT i
      #ENDIF

      #IF 00
      FOR i = 00 TO @pU_CSP_PARM.N_LDU - 01

         S_CSP   += TRIM$("L[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].L) + " " + _
                    TRIM$("D[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].D) + " " + _
                    TRIM$("U[" + STR$(i+1) + "]") + " " + STR$(@[email protected]_LDU[i].U) + $CRLF

      NEXT i
      #ENDIF

      #IF -1
      FOR i_INTP = 00 TO @pU_CSP_PARM.N_INTP - 01

         S_CSP   += "[i INTP, t_INTP, y_INTP] " + STR$(i_INTP) + " " + STR$(@[email protected]_INTP[i_INTP]) + " " + STR$(@[email protected]_INTP[i_INTP]) + $CRLF

      NEXT i_INTP
      #ENDIF

      FUNCTION = S_CSP

   END IF

END FUNCTION


FUNCTION CSP_TEST      (                                          )  AS LONG

   '  Routine computes an cubic spline test function

   DIM N_CSP                                             AS LONG
   DIM U_CSP()                                           AS UDT_CSP
   DIM U_CSP_PARM                                        AS UDT_CSP_PARM

   DIM N_DATA                                            AS LONG
   DIM t_DATA()                                          AS DOUBLE
   DIM y_DATA()                                          AS DOUBLE

   DIM N_INTP                                            AS LONG
   DIM t_INTP()                                          AS DOUBLE
   DIM y_INTP()                                          AS DOUBLE

   DIM N_LDU                                             AS LONG
   DIM U_LDU()                                           AS UDT_LDU


   '  Create a test input data set
   '  Note the assumption is that the load data is sorted
   '     0.9   0
   '     1.3   1.5
   '     1.9   1.85
   '     2.1   2.1
   N_DATA = 14
   REDIM t_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN t_DATA() = 0.00, 1.00, 3.00, 4.00, 5.00, 6.00, 8.00, 11.00, 12.00, 13.00, 14.00, 16.00, 17.00, 19.00
   REDIM y_DATA(00 TO N_DATA - 01) : ARRAY ASSIGN y_DATA() = 0.00, 2.60, 23.16, 27.57, 24.26, 16.63, 30.41, 47.20, 50.03, 60.33, 59.89, 71.18, 84.27, 77.69
   U_CSP_PARM.N_DATA    = N_DATA
   U_CSP_PARM.pt_DATA   = VARPTR(t_DATA(00))
   U_CSP_PARM.py_DATA   = VARPTR(y_DATA(00))

   '  Allocate the LDU solution data structure
   N_LDU = N_DATA
   REDIM U_LDU(00 TO N_LDU - 01)
   U_CSP_PARM.N_LDU     = N_DATA
   U_CSP_PARM.pU_LDU    = VARPTR(U_LDU(00))

   '  Allocate the cubic spline coefficients to the number of data points
   N_CSP = N_DATA
   REDIM U_CSP(00 TO N_CSP - 01)

   '  Set the LDU solve parameters
   U_CSP_PARM.U_LDU_PARM.N_ITR_MAX = 5
   U_CSP_PARM.U_LDU_PARM.e_ITR_MIN = 1.0E-3#

   '  Create a test input interpolation set
   '  Note this is larger than the input data set to create smoothness
   N_INTP = 100
   REDIM t_INTP(00 TO N_INTP - 01)
   REDIM y_INTP(00 TO N_INTP - 01)
   U_CSP_PARM.N_INTP    = N_INTP
   U_CSP_PARM.pt_INTP   = VARPTR(t_INTP(00))
   U_CSP_PARM.py_INTP   = VARPTR(y_INTP(00))


   '  Solve for cubic spline
   CSP_SOLVE      N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)
   CSP_SET_CLIP   N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)

   #IF %T
   MSGBOX      CSP_GET_S(N_CSP, VARPTR(U_CSP(00)), VARPTR(U_CSP_PARM)), %MB_ICONINFORMATION OR %MB_TASKMODAL OR %MB_TOPMOST, $S_MSGBOX_HDR
   #ENDIF

END FUNCTION

**Mannish Bhandari** · 20 Nov 2021, 12:59 PM

Can you show how to plot out the interpolated data graphically?

Announcement

Forum Guidelines

Cubic Spline

Cubic Spline

Comment

Comment