How to implement a spline subroutine in my code

Help
1

I have tried to implement a spline type subroutine. It correctly takes data from my text file (17x2 matrix), I write the following subroutine (provided by chatgpt, I’m just starting to program).

C gL

  DO i=1,N
    r_val=r_1(i)
    idx=1
    min_diff=ABS(r_val - matrizB(1,1))

    DO iB=2,NB
      diff=ABS(r_val - matriz(iB,1))
      IF (diff< min_diff) THEN
        idx=iB
        min_diff = diff
      END IF
    END DO

    CALL CUBIC_INTERPOLATION (r_val, matrizB(idx-1,1), matrizB(idx,1), matrizB(idx-1,2), matrizB(idx,2),h)
    interpolated_values(i)=h
  END DO

  SUBROUTINE CUBIC_INTERPOLATION(x,x0,x1,y0,y1,h)
    REAL, INTENT(IN) :: x,x0, x1, y0, y1
    REAL, INTENT(OUT) :: h

    REAL :: t, t2, t3

    t = (x -x0)/(x1-x0)
    t2=t*t
    t3=t2*t

    h =y0*(2*t3 - 3*t2+1) + y1*(-2*t3+3*t2)+(y1-y0)*(t3-2*t2+t)

  END SUBROUTINE CUBIC_INTERPOLATION

This is only the problem part of my extensive fortran program.

And finally I make a call to the subroutine to create me a vector of 58 elements from another of the same dimension.

EDIT: I’ve trying some other interpolation, and this version is the one with less errors. I think a spline or a cubic interpolation could work. This is my data (Tabla2.txt):
0.00 1.010
0.10 1.010
0.15 1.018
0.25 1.030
0.50 1.030
0.75 1.020
1.00 1.000
1.50 0.937
2.00 0.857
3.00 0.689
4.00 0.538
5.00 0.409
6.00 0.313
7.00 0.232
8.00 0.176
9.00 0.134
10.00 0.0957

EDIT 2: Following the recommendation vmagnin, I’ll provide error message:

prueba2.f:358:72:

358 | CALL CUBIC_INTERPOLATION (r_val, matrizB(idx-1,1), matrizB(idx,1), matrizB(idx-1,2), matrizB(idx,2),h)
| 1
Error: Invalid form of array reference at (1)

What I expect. I expect to obtain a interpolation of Tabla2.txt, then create a array g_l what take every element of r_1(58 elements) and give me the array g_l(58 elements)

** I’m using Windows 11 GNU Fortran (MinGW-W64 x86_64-ucrt-posix-seh, built by Brecht Sanders) 12.2.0

2

Welcome @rasty0013 to the Fortran Discourse!

The title of your post sounds like a question but you don’t explain your problem. If you expect help, you should explain precisely what problem occurs (what is expected and what is really obtained), and give any compiler error message or execution error message that appears. Your OS and the compiler you use may be also interesting.

3

If you need a subroutine that does interpolation why you don’t check the bspline-fortran package of @jacobwilliams

3 Likes
4

Also the OP should look at Carl De Boor’s PPPACK,
https://people.sc.fsu.edu/~jburkardt//f_src/pppack/pppack.html

Plus a better option than a pure cubic spline is a monotonic piece-wise cubic routine like PCHIPS (I think @jacobwilliams has a refactored version) and the Akima family of splines from the ACM repository. Also, just about all computer graphics and CAGD books that cover the math will have a spline routine of some kind.

Also, NASA’s cfdtools github directory has several spline routines in the interplib directory.

5

Your code does not get far before failing !
For the first itteration of the first DO loop : Note that when i = 2, referencing “y_vals(i-2)” is addressing out of bounds and “h(i-1)” is not previously defined.
Also “alpha” is a double backward difference, which is probably not the best approach.

You need to check the end points of the DO loop for valid and defined references.

A possible way to address the range of valid values is you could have local arrays, dx(n) and dy(n), and possibly estimate end points to avoid addressing arrays out of bounds, or modify the i range for addressing vals (1:n) to cope with cubic interpolation.

The following code could get you past the first DO loop

REAL :: dx(n), dy(n)

DO i = 2, n
  dx(i) = x_vals(i) - x_vals(i-1)
  dy(i) = y_vals(i) - y_vals(i-1)
END DO
!   note dx and dy only defined for the range 2:n
DO i = 2, n-1
  alpha(i) = 3.0 * ( dy(i+1)/ dx(i+1) - dy(i)/ dx(i) )
END DO
!   now defined alpha (2:n-1)  but for different i:i+1 interval to your code.

Don’t give up with your approach, as you have to start somewhere.

1 Like
6

thanks for the answer, I have edited my question.

1 Like
7

Your interpolation is basically based on two points (x0,y0) and (x1,y1),
All you can generate with this info is a linear interpolation.
You can not generate a cubic interpolation, without further knowledge of the slopes (dy/dx) at these two points.

You are using 3 shape functions, based on “t”; the spacing of x between x0 and x1.

This formula is basically:
yt = y0 * Shape0(t) + y1 * Shape1(t) + dy * shape01(t).

Unfortunately, any cubic interpoplation is based on the suitability of the 3 shape functions to dy/dx at the two end points.
shape01(t) must be related to the end slopes, not assumed.

If you plot your Shape0(t), Shape1(t) and shape01(t) functions in Excel, you can identify how shape01(t) is a particular solution for dy/dx at the 2 end points, rather than a general solution.

You need to provide (x,y, dy/dx) at 2 points to achieve cubic interpolation, which supplies estimates of dy/dx at the 2 points you are using.

Try Google, rather than ChatGPT

Alternatively, there must be a 4-point Lagrange interpolation formula for non even spacing of xi

8

There are a lot of good interpolation libraries available. No need to roll your own.

  • bspline-fortran (1D-6D B-spline interpolation)
  • finterp (1D-6D linear or nearest neighbor interpolation)
  • regridpack (1D-4D linear and cubic interpolation)
  • splpak (ND least-squares cubic spline fitting)
  • PCHIP (1D piecewise cubic Hermite interpolation)
  • fitpack (package for curve and surface fitting)
3 Likes