「スプライン補間」の編集履歴（バックアップ）一覧はこちら

「スプライン補間」(2010/02/17 (水) 12:15:43) の最新版変更点

追加された行は緑色になります。

削除された行は赤色になります。

#asciiart(blockquote){
#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
}