#N88BASIC
/***********************************************************
spline.c -- スプライン補間
***********************************************************/
/* 非周期関数用 */
Const N = 5
Dim x[ELM(N)] = [ 0, 10, 20, 30, 40 ] As Double
Dim y[ELM(N)] = [ 610.66, 1227.4, 2338.1, 4244.9, 7381.2 ] As Double
Dim z[ELM(N)] As Double
Dim h[ELM(N)] As Double, d[ELM(N)] As Double'static
Sub maketable(x As *Double, y As *Double, z As *Double)
Dim i As Long
Dim t As Double
z[0] = 0: z[N - 1] = 0 /* 両端点での y''(x) / 6 */
For i=0 To N-2
h[i ] = x[i + 1] - x[i]
d[i + 1] = (y[i + 1] - y[i]) / h[i]
Next
z[1] = d[2] - d[1] - h[0] * z[0]
d[1] = 2 * (x[2] - x[0])
For i = 1 To N - 3
t = h[i] / d[i]
z[i + 1] = d[i + 2] - d[i + 1] - z[i] * t
d[i + 1] = 2 * (x[i + 2] - x[i]) - h[i] * t
Next
z[N - 2] = z[N - 2] - h[N - 2] * z[N - 1]
For i = N - 2 To 1 Step -1
z[i] = (z[i] - h[i] * z[i + 1]) / d[i]
Next
End Sub
Function spline(t As Double, x As *Double, y As *Double, z As *Double) As Double
Dim i As Long, j As Long, k As Long
Dim d As Double, h As Double
i = 0: j = N - 1
While i < j
k = (i + j) / 2
If x[k] < t Then i = k + 1 Else j = k
Wend
If i > 0 Then i=i-1
h = x[i + 1] - x[i]: d = t - x[i]
spline = (((z[i + 1] - z[i]) * d / h + z[i] * 3) * d _
+ ((y[i + 1] - y[i]) / h _
- (z[i] * 2 + z[i + 1]) * h)) * d + y[i]
End Function
'main
Dim i As Long
maketable(x, y, z)
For i = 10 To 30
Print i, spline(i,x,y,z)
Next
