エルミート補間

home > さまざまな言語で数値計算 > 関数の近似 >

エルミート補間

２個の関数値 $ f(x_0), f(x_1) $ と, それぞれの微分係数 $ f'(x_0), f'(x_1) $ とを与えられたとき, 与えられた２個の関数値を通る３次式を求めるには, まず次のような表を作る.

このとき,

あとは, ニュートン補間と同じ.

第$n$差分商を

で表すと, 与えられた $n$点を通る $2n-1$次式は次のように表すことができる.

この式を使って, 与えられた点以外の点の値を求める.

例題として,

を近似する.

VBScript

Option Explicit

'データ点の数 - 1
Private Const N   =  6
Private Const Nx2 = 13

Dim x():  ReDim x(N)
Dim y():  ReDim y(N)
Dim yd(): ReDim yd(N)

'1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
Dim i
For i = 0 To N
    Dim d1: d1 = i * 1.5 - 4.5
    x(i)  = d1
    y(i)  = f(d1)
    yd(i) = fd(d1)
Next

'差分商の表を作る
Dim z(): ReDim z(Nx2)
Dim d(): ReDim d(Nx2, Nx2)
For i = 0 To Nx2
    Dim j: j = i \ 2
    z(i)   = x(j)
    d(0,i) = y(j)
Next

For i = 1 To Nx2
    For j = 0 To (Nx2 - i)
        If i = 1 And j Mod 2 = 0 Then
            d(i,j) = yd(j \ 2)
        Else
            d(i,j) = (d(i-1,j+1) - d(i-1,j)) / (z(j+i) - z(j))
        End If
    Next
Next

'ｎ階差分商
Dim a(): ReDim a(Nx2)
For j = 0 To Nx2
    a(j) = d(j,0)
Next

'0.5刻みで 与えられていない値を補間
For i = 0 To 18
    d1 = i * 0.5 - 4.5
    Dim d2: d2 = f(d1)
    Dim d3: d3 = hermite(d1, z, a)

    '元の関数と比較
    WScript.StdOut.Write     Right(Space(5) & FormatNumber(d1,      2, -1, 0, 0), 5) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d2,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.Write     Right(Space(8) & FormatNumber(d3,      5, -1, 0, 0), 8) & vbTab
    WScript.StdOut.WriteLine Right(Space(8) & FormatNumber(d2 - d3, 5, -1, 0, 0), 8)
Next

'元の関数
Private Function f(ByVal x)
    f = x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2)
End Function
'導関数
Private Function fd(ByVal x)
    fd = 1 - (x ^ 2) / 2 + (x ^ 4) / (4 * 3 * 2)
End Function

'Hermite (エルミート) 補間
Private Function hermite(ByVal d, ByVal z(), ByVal a())
    Dim sum: sum = a(0)
    Dim i, j
    For i = 1 To Nx2
        Dim prod: prod = a(i)
        For j = 0 To (i - 1)
            If j <> i Then
                prod = prod * (d - z(j))
            End If
        Next
        sum = sum + prod
    Next
    hermite = sum
End Function

Z:\>cscript //nologo Z:\0704.vbs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

JScript

// データ点の数
var N   =  7
var Nx2 = 14

var x  = []
var y  = []
var yd = []

// 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
for (var i = 0; i < N; i++)
{
    var d1 = i * 1.5 - 4.5
    x[i]   = d1
    y[i]   = f(d1)
    yd[i]  = fd(d1)
}

// 差分商の表を作る
var z = []
var d = []
d[0]  = []
for (var i = 0; i < Nx2; i++)
{
    var j   = parseInt(i / 2)
    z[i]    = x[j]
    d[0][i] = y[j]
}
for (var i = 1; i < Nx2; i++)
{
    d[i] = []
    for (var j = 0; j < Nx2 - i; j++)
    {
        if (i == 1 && j % 2 == 0)
            d[i][j] = yd[parseInt(j / 2)]
        else
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j])
    }
}

// ｎ階差分商
var a = []
for (var j = 0; j < Nx2; j++)
    a[j] = d[j][0]
// 0.5刻みで 与えられていない値を補間
for (var i = 0; i <= 18; i++)
{
    var d1 = i * 0.5 - 4.5
    var d2 = f(d1)
    var d3 = hermite(d1, z, a)

    // 元の関数と比較
    WScript.StdOut.Write(("     "    + d1.toFixed(2)       ).slice(-5) + "\t")
    WScript.StdOut.Write(("        " + d2.toFixed(5)       ).slice(-8) + "\t")
    WScript.StdOut.Write(("        " + d3.toFixed(5)       ).slice(-8) + "\t")
    WScript.StdOut.Write(("        " + (d2 - d3).toFixed(5)).slice(-8) + "\n")
}

// 元の関数
function f(x)
{
    return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2)
}
// 導関数
function fd(x)
{
    return 1 - Math.pow(x,2) / 2 + Math.pow(x,4) / (4 * 3 * 2)
}

// Hermite (エルミート) 補間
function hermite(d, z, a) {
    var sum = a[0]
    for (var i = 1; i < Nx2; i++)
    {
        var prod = a[i]
        for (var j = 0; j < i; j++)
            prod *= (d - z[j])
        sum += prod
    }

    return sum
}

Z:\>cscript //nologo Z:\0704.js
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

PowerShell

# データ点の数
set-variable -option constant -name N   -value  7
set-variable -option constant -name Nx2 -value 14

# 元の関数
function f($x)
{
    $x - [Math]::Pow($x, 3) / (3 * 2) + [Math]::Pow($x, 5) / (5 * 4 * 3 * 2)
}
# 導関数
function fd($x)
{
    1 - [Math]::Pow($x, 2) / 2 + [Math]::Pow($x, 4) / (4 * 3 * 2)
}

# Hermite (エルミート) 補間
function hermite($d, $z, $a)
{
    $sum = $a[0]
    foreach ($i in 1..($Nx2 - 1))
    {
        $prod = $a[$i]
        foreach ($j in 0..($i - 1))
        {
            $prod *= ($d - $z[$j])
        }
        $sum += $prod
    }
    $sum
}

$x  = New-Object double[] $N
$y  = New-Object double[] $N
$yd = New-Object double[] $N

# 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
foreach ($i in 0..($N - 1))
{
    $d1      = $i * 1.5 - 4.5
    $x[$i]   = $d1
    $y[$i]   = f($d1)
    $yd[$i]  = fd($d1)
}

# 差分商の表を作る
$z = New-Object  double[]   $Nx2
$d = New-Object "double[,]" $Nx2,$Nx2
foreach ($i in 0..($Nx2 - 1))
{
    $j       = [Math]::Floor($i / 2)
    $z[$i]   = $x[$j]
    $d[0,$i] = $y[$j]
}
foreach ($i in 1..($Nx2 - 1))
{
    foreach ($j in 0..($Nx2 - $i - 1))
    {
        if ($i -eq 1 -and $j % 2 -eq 0)
        {
            $d[$i,$j] = $yd[[Math]::Floor($j / 2)]
        }
        else
        {
            $d[$i,$j] = ($d[($i-1),($j+1)] - $d[($i-1),$j]) / ($z[($j+$i)] - $z[$j])
        }
    }
}
# ｎ階差分商
$a = New-Object double[] $Nx2
foreach ($j in 0..($Nx2 - 1))
{
    $a[$j] = $d[$j,0]
}

# 0.5刻みで 与えられていない値を補間
foreach ($i in 0..18)
{
    $d1 = $i * 0.5 - 4.5
    $d2 = f($d1)
    $d3 = (hermite $d1 $z $a)

    # 元の関数と比較
    Write-Host ([String]::Format("{0,5:F2}`t{1,8:F5}`t{2,8:F5}`t{3,8:F5}", $d1, $d2, $d3, ($d2 - $d3)))
}

Z:\>powershell -file Z:\0704.ps1
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

Perl

# データ点の数 - 1
use constant N   =>  6;
use constant Nx2 => 13;

my @x  = ();
my @y  = ();
my @yd = ();

# 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
for $i (0..N)
{
    my $d1  = $i * 1.5 - 4.5;
    $x[$i]  = $d1;
    $y[$i]  = f($d1);
    $yd[$i] = fd($d1);
}

# 差分商の表を作る
my @z = ();
my @d = ();
for $i (0..Nx2)
{
    $j        = int($i / 2);
    $z[$i]    = $x[$j];
    $d[0][$i] = $y[$j];
}
for $i (1..Nx2)
{
    for $j (0..(Nx2 - $i))
    {
        if ($i == 1 && $j % 2 == 0)
        {
            $d[$i][$j] = $yd[int($j / 2)];
        }
        else
        {
            $d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($z[$j+$i] - $z[$j]);
        }
    }
}

# ｎ階差分商
my @a = ();
for $j (0..Nx2)
{
    $a[$j] = $d[$j][0];
}
# 0.5刻みで 与えられていない値を補間
for $i (0..18)
{
    my $d1 = $i * 0.5 - 4.5;
    my $d2 = f($d1);
    my $d3 = hermite($d1, \@z, \@a);

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d1, $d2, $d3, $d2 - $d3);
}

# 元の関数
sub f
{
    my ($x) = @_;
    $x - ($x ** 3) / (3 * 2) + ($x ** 5) / (5 * 4 * 3 * 2);
}
# 導関数
sub fd
{
    my ($x) = @_;
    1 - ($x ** 2) / 2 + ($x ** 4) / (4 * 3 * 2);
}

# Hermite (エルミート) 補間
sub hermite
{
    my ($d, $z, $a) = @_;
    my $sum = $$a[0];
    for $i (1..Nx2)
    {
        my $prod = $$a[$i];
        for $j (0..($i - 1))
        {
            $prod *= ($d - $$z[$j]);
        }
        $sum += $prod;
    }
    $sum;
}

Z:\>perl Z:\0704.pl
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

PHP

<?php
# データ点の数 - 1
define("N",    6);
define("Nx2", 13);

$x  = array();
$y  = array();
$yd = array();

# 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
foreach (range(0, N) as $i)
{
    $d1     = $i * 1.5 - 4.5;
    $x[$i]  = $d1;
    $y[$i]  = f($d1);
    $yd[$i] = fd($d1);
}

# 差分商の表を作る
$z = array();
$d = array();
foreach (range(0, Nx2) as $i)
{
    $j        = (int)($i / 2);
    $z[$i]    = $x[$j];
    $d[0][$i] = $y[$j];
}
foreach (range(1, Nx2) as $i)
{
    foreach (range(0, Nx2 - $i) as $j)
    {
        if ($i == 1 && $j % 2 == 0)
            $d[$i][$j] = $yd[(int)($j / 2)];
        else
            $d[$i][$j] = ($d[$i-1][$j+1] - $d[$i-1][$j]) / ($z[$j+$i] - $z[$j]);
    }
}
# ｎ階差分商
$a = array();
foreach (range(0, Nx2) as $j)
{
    $a[$j] = $d[$j][0];
}

# 0.5刻みで 与えられていない値を補間
foreach (range(0, 18) as $i)
{
    $d1 = $i * 0.5 - 4.5;
    $d2 = f($d1);
    $d3 = hermite($d1, $z, $a);

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", $d1, $d2, $d3, $d2 - $d3);
}

# 元の関数
function f($x)
{
    return $x - pow($x, 3) / (3 * 2) + pow($x, 5) / (5 * 4 * 3 * 2);
}
# 導関数
function fd($x)
{
    return 1 - pow($x, 2) / 2 + pow($x, 4) / (4 * 3 * 2);
}

# Hermite (エルミート) 補間
function hermite($d, $z, $a)
{
    $sum = $a[0];
    foreach (range(1, Nx2) as $i)
    {
        $prod = $a[$i];
        foreach (range(0, $i - 1) as $j)
        {
            $prod *= ($d - $z[$j]);
        }
        $sum += $prod;
    }
    return $sum;
}
?>

Z:\>php Z:\0704.php
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Python

# coding: Shift_JIS

# データ点の数
N   =  7
Nx2 = 14

# 元の関数
def f(x):
    return x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2)

# 導関数
def fd(x):
    return 1 - (x ** 2) / 2 + (x ** 4) / (4 * 3 * 2)

# Hermite (エルミート) 補間
def hermite(d, z, a):
    sum = a[0]
    for i in range(1, Nx2):
        prod = a[i]
        for j in range(0, i):
            prod *= (d - z[j])
        sum += prod
    return sum

x  = [0 for i in range(N)]
y  = [0 for i in range(N)]
yd = [0 for i in range(N)]

# 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
for i in range(0, N):
    d = i * 1.5 - 4.5
    x[i]  = d
    y[i]  = f(d)
    yd[i] = fd(d)

# 差分商の表を作る
z =  [0 for i in range(Nx2)]
d = [[0 for j in range(Nx2)] for i in range(Nx2)]
for i in range(0, Nx2):
    j       = int(i / 2)
    z[i]    = x[j]
    d[0][i] = y[j]

for i in range(1, Nx2):
    for j in range(0, Nx2 - i):
        if (i == 1 and j % 2 == 0):
            d[i][j] = yd[int(j / 2)]
        else:
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j])

# ｎ階差分商
a = [0 for i in range(Nx2)]
for j in range(0, Nx2):
    a[j] = d[j][0]

# 0.5刻みで 与えられていない値を補間
for i in range(0, 19):
    d1 = i * 0.5 - 4.5
    d2 = f(d1)
    d3 = hermite(d1, z, a)

    # 元の関数と比較
    print "%5.2f\t%8.5f\t%8.5f\t%8.5f" % (d1, d2, d3, d2 - d3)

Z:\>python Z:\0704.py
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Ruby

# データ点の数 - 1
N   =  6
Nx2 = 13

# 元の関数
def f(x)
    x - (x ** 3) / (3 * 2) + (x ** 5) / (5 * 4 * 3 * 2)
end

# 導関数
def fd(x)
    1 - (x ** 2) / 2 + (x ** 4) / (4 * 3 * 2)
end

# Hermite (エルミート) 補間
def hermite(d, z, a)
    sum = a[0]
    (1..Nx2).each do |i|
        prod = a[i]
        (0..(i - 1)).each do |j|
            if j != i
                prod *= (d - z[j])
            end
        end
        sum += prod
    end
    sum
end

x  = Array.new(N)
y  = Array.new(N)
yd = Array.new(N)

# 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
(0..N).each do |i|
    d = i * 1.5 - 4.5
    x[i]  = d
    y[i]  = f(d)
    yd[i] = fd(d)
end

# 差分商の表を作る
z = Array.new(Nx2)
d = Array.new(Nx2+1) { Array.new(Nx2) }
(0..Nx2).each do |i|
    j = i / 2
    z[i]    = x[j]
    d[0][i] = y[j]
end

(1..Nx2).each do |i|
    (0..(Nx2-i)).each do |j|
        if i == 1 && j % 2 == 0
            d[i][j] = yd[j / 2]
        else
            d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j])
        end
    end
end

# ｎ階差分商
a = Array.new(Nx2)
(0..Nx2).each do |j|
    a[j] = d[j][0]
end

# 0.5刻みで 与えられていない値を補間
(0..18).each do |i|
    d1 = i * 0.5 - 4.5
    d2 = f(d1)
    d3 = hermite(d1, z, a)

    # 元の関数と比較
    printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d1, d2, d3, d2 - d3)
end

Z:\>ruby Z:\0704.rb
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Groovy

TOP

Pascal

program Pas0704(arg);
{$MODE delphi}

uses
    SysUtils, Math;

const
    // データ点の数 - 1
    N   =  6;
    Nx2 = 13;

// 元の関数
function f(x:Double):Double;
begin
    result := x - power(x,3) / (3 * 2) + power(x,5) / (5 * 4 * 3 * 2);
end;

// 導関数
function fd(x:Double):Double;
begin
    result := 1 - power(x,2) / 2 + power(x,4) / (4 * 3 * 2);
end;

// Hermite (エルミート) 補間
function hermite(d:Double; z:array of Double; a:array of Double):Double;
var
    sum, prod :Double;
    i, j      :Integer;
begin
    sum := a[0];
    for i := 1 to High(a) do
    begin
        prod := a[i];
        for j := Low(z) to i-1 do
            prod := prod * (d - z[j]);
        sum := sum + prod;
    end;
    result := sum;
end;

var
    i, j :Integer;
    x  :array [0..N] of Double;
    y  :array [0..N] of Double;
    yd :array [0..N] of Double;
    z  :array [0..Nx2] of Double;
    a  :array [0..Nx2] of Double;
    d  :array [0..Nx2, 0..Nx2] of Double;
    d1, d2, d3 :Double;
begin
    // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for i := Low(x) to High(x) do
    begin
        d1    := i * 1.5 - 4.5;
        x[i]  := d1;
        y[i]  := f(d1);
        yd[i] := fd(d1);
    end;

    // 差分商の表を作る
    for i := Low(z) to High(z) do
    begin
        j      := i div 2;
        z[i]   := x[j];
        d[0,i] := y[j];
    end;

    for i := 1 to High(z) do
    begin
        for j := Low(z) to High(z)-i do
        begin
            if (i = 1) and (j mod 2 = 0) then
                d[i,j] := yd[j div 2]
            else
                d[i,j] := (d[i-1,j+1] - d[i-1,j]) / (z[j+i] - z[j]);
        end;
    end;

    // ｎ階差分商
    for j := Low(a) to High(a) do
        a[j] := d[j,0];

    // 0.5刻みで 与えられていない値を補間
    for i := 0 to 18 do
    begin
        d1 := i * 0.5 - 4.5;
        d2 := f(d1);
        d3 := hermite(d1, z, a);

        // 元の関数と比較
        writeln(format('%5.2f'#9'%8.5f'#9'%8.5f'#9'%8.5f', [d1, d2, d3, d2 - d3]));
    end;
end.

Z:\>fpc -v0 -l- Pas0704.pp

Z:\>Pas0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

Ada

with TEXT_IO, Ada.Long_Float_Text_IO;
use  TEXT_IO, Ada.Long_Float_Text_IO;

procedure Ada0704 is
    -- データ点の数 - 1
    N :   Constant Integer :=  6;
    Nx2 : Constant Integer := 13;

    type Long_Float_ArrayN   is array (0..N)   of Long_Float;
    type Long_Float_ArrayNx2 is array (0..Nx2) of Long_Float;
    x  : Long_Float_ArrayN;
    y  : Long_Float_ArrayN;
    yd : Long_Float_ArrayN;
    z  : Long_Float_ArrayNx2;
    a  : Long_Float_ArrayNx2;
    d : array (0..Nx2, 0..Nx2) of Long_Float;
    d1, d2, d3 : Long_Float;
    k : Integer;

    -- 元の関数
    function f(x:Long_Float) return Long_Float is
    begin
        return x - Long_Float(x ** 3) / Long_Float(3 * 2) + Long_Float(x ** 5) / Long_Float(5 * 4 * 3 * 2);
    end f;

    -- 導関数
    function fd(x:Long_Float) return Long_Float is
    begin
        return 1.0 - Long_Float(x ** 2) / Long_Float(2) + Long_Float(x ** 4) / Long_Float(4 * 3 * 2);
    end fd;

    -- Hermite (エルミート) 補間
    function hermite(d:Long_Float; z:Long_Float_ArrayNx2; a:Long_Float_ArrayNx2) return Long_Float is
        sum, prod :Long_Float;
    begin
        sum := a(0);
        for i in 1 .. a'Last loop
            prod := a(i);
            for j in z'First .. i-1 loop
                prod := prod * (d - z(j));
            end loop;
            sum := sum + prod;
        end loop;
        return sum;
    end;
begin
    -- 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for i in x'Range loop
        d1    := Long_Float(i) * 1.5 - 4.5;
        x(i)  := d1;
        y(i)  := f(d1);
        yd(i) := fd(d1);
    end loop;

    -- 差分商の表を作る
    for i in z'Range loop
        k      := i / 2;
        z(i)   := x(k);
        d(0,i) := y(k);
    end loop;

    for i in 1 .. z'Last loop
        for j in z'First .. z'Last-i loop
            if (i = 1) and (j mod 2 = 0) then
                d(i,j) := yd(j / 2);
            else
                d(i,j) := (d(i-1,j+1) - d(i-1,j)) / (z(j+i) - z(j));
            end if;
        end loop;
    end loop;

    -- ｎ階差分商
    for j in a'Range loop
        a(j) := d(j,0);
    end loop;

    -- 0.5刻みで 与えられていない値を補間
    for i in 0..18 loop
        d1  := Long_Float(i) * 0.5 - 4.5;
        d2 := f(d1);
        d3 := hermite(d1, z, a);

        -- 元の関数と比較
        Put(d1,      Fore=>2, Aft=>2, Exp=>0);
        Put(Ascii.HT);
        Put(d2,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d3,      Fore=>3, Aft=>5, Exp=>0);
        Put(Ascii.HT);
        Put(d2 - d3, Fore=>3, Aft=>5, Exp=>0);
        New_Line;
    end loop;
end Ada0704;

xxxxxx@yyyyyy /Z
$ gnatmake Ada0704.adb

xxxxxx@yyyyyy /Z
$ Ada0704
-4.50    -4.68984    -4.68984     0.00000
-4.00    -1.86667    -1.86667     0.00000
-3.50    -0.73099    -0.73099     0.00000
-3.00    -0.52500    -0.52500     0.00000
-2.50    -0.70964    -0.70964     0.00000
-2.00    -0.93333    -0.93333    -0.00000
-1.50    -1.00078    -1.00078     0.00000
-1.00    -0.84167    -0.84167     0.00000
-0.50    -0.47943    -0.47943     0.00000
 0.00     0.00000     0.00000    -0.00000
 0.50     0.47943     0.47943    -0.00000
 1.00     0.84167     0.84167    -0.00000
 1.50     1.00078     1.00078    -0.00000
 2.00     0.93333     0.93333    -0.00000
 2.50     0.70964     0.70964    -0.00000
 3.00     0.52500     0.52500    -0.00000
 3.50     0.73099     0.73099    -0.00000
 4.00     1.86667     1.86667    -0.00000
 4.50     4.68984     4.68984    -0.00000

TOP

VB.NET

Module VB0704
    'データ点の数 - 1
    Private Const N   As Integer =  6
    Private Const Nx2 As Integer = 13

    Public Sub Main()
        Dim x(N)  As Double
        Dim y(N)  As Double
        Dim yd(N) As Double

        '1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
        For i As Integer = 0 To N
            Dim d1 As Double = i * 1.5 - 4.5
            x(i)  = d1
            y(i)  = f(d1)
            yd(i) = fd(d1)
        Next

        '差分商の表を作る
        Dim z(Nx2)      As Double
        Dim d(Nx2, Nx2) As Double
        For i As Integer = 0 To Nx2
            Dim j As Integer = i \ 2
            z(i)   = x(j)
            d(0,i) = y(j)
        Next

        For i As Integer = 1 To Nx2
            For j As Integer = 0 To (Nx2 - i)
                If i = 1 AndAlso j Mod 2 = 0 Then
                    d(i,j) = yd(j \ 2)
                Else
                    d(i,j) = (d(i-1,j+1) - d(i-1,j)) / (z(j+i) - z(j))
                End If
            Next
        Next

        'ｎ階差分商
        Dim a(Nx2) As Double
        For j As Integer = 0 To Nx2
            a(j) = d(j,0)
        Next

        '0.5刻みで 与えられていない値を補間
        For i As Integer = 0 To 18
            Dim d1 As Double = i * 0.5 - 4.5
            Dim d2 As Double = f(d1)
            Dim d3 As Double = hermite(d1, z, a)

            '元の関数と比較
            Console.WriteLine(String.Format("{0,5:F2}{4}{1,8:F5}{4}{2,8:F5}{4}{3,8:F5}", d1, d2, d3, d2 - d3, vbTab))
        Next
    End Sub

    '元の関数
    Private Function f(ByVal x As Double) As Double
        Return x - (x ^ 3) / (3 * 2) + (x ^ 5) / (5 * 4 * 3 * 2)
    End Function
    '導関数
    Private Function fd(ByVal x As Double) As Double
        Return 1 - (x ^ 2) / 2 + (x ^ 4) / (4 * 3 * 2)
    End Function

    'Hermite (エルミート) 補間
    Private Function hermite(ByVal d As Double, ByVal z() As Double, ByVal a() As Double) As Double
        Dim sum As Double = a(0)
        For i As Integer = 1 To Nx2
            Dim prod As Double = a(i)
            For j As Integer = 0 To (i - 1)
                prod *= (d - z(j))
            Next
            sum += prod
        Next
        Return sum
    End Function
End Module

Z:\>vbc -nologo VB0704.vb

Z:\>VB0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

C#

using System;

public class CS0704
{
    // データ点の数
    private const int N   =  7;
    private const int Nx2 = 14;

    public static void Main()
    {
        double[] x  = new double[N];
        double[] y  = new double[N];
        double[] yd = new double[N];

        // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
        for (int i = 0; i < N; i++)
        {
            double d1 = i * 1.5 - 4.5;
            x[i]  = d1;
            y[i]  = f(d1);
            yd[i] = fd(d1);
        }

        // 差分商の表を作る
        double[]  z = new double[Nx2];
        double[,] d = new double[Nx2,Nx2];
        for (int i = 0; i < Nx2; i++)
        {
            int j  = i / 2;
            z[i]   = x[j];
            d[0,i] = y[j];
        }
        for (int i = 1; i < Nx2; i++)
        {
            for (int j = 0; j < Nx2 - i; j++)
            {
                if (i == 1 && j % 2 == 0)
                    d[i,j] = yd[j / 2];
                else
                    d[i,j] = (d[i-1,j+1] - d[i-1,j]) / (z[j+i] - z[j]);
            }
        }

        // ｎ階差分商
        double[] a = new double[Nx2];
        for (int j = 0; j < Nx2; j++)
            a[j] = d[j,0];

        // 0.5刻みで 与えられていない値を補間
        for (int i = 0; i <= 18; i++)
        {
            double d1 = i * 0.5 - 4.5;
            double d2 = f(d1);
            double d3 = hermite(d1, z, a);

            // 元の関数と比較
            Console.WriteLine(string.Format("{0,5:F2}\t{1,8:F5}\t{2,8:F5}\t{3,8:F5}", d1, d2, d3, d2 - d3));
        }
    }

    // 元の関数
    private static double f(double x)
    {
        return x - Math.Pow(x,3) / (3 * 2) + Math.Pow(x,5) / (5 * 4 * 3 * 2);
    }
    // 導関数
    private static double fd(double x)
    {
        return 1 - Math.Pow(x,2) / 2 + Math.Pow(x,4) / (4 * 3 * 2);
    }

    // Hermite (エルミート) 補間
    private static double hermite(double d, double[] z, double[] a)
    {
        double sum = a[0];
        for (int i = 1; i < Nx2; i++)
        {
            double prod = a[i];
            for (int j = 0; j < i; j++)
                prod *= (d - z[j]);
            sum += prod;
        }

        return sum;
    }
}

Z:\>csc -nologo CS0704.cs

Z:\>CS0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

Java

public class Java0704 {

    // データ点の数
    private static final int N   =  7;
    private static final int Nx2 = 14;

    public static void main(String []args) {
        double[] x  = new double[N];
        double[] y  = new double[N];
        double[] yd = new double[N];

        // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
        for (int i = 0; i < N; i++) {
            double d1 = i * 1.5 - 4.5;
            x[i]  = d1;
            y[i]  = f(d1);
            yd[i] = fd(d1);
        }

        // 差分商の表を作る
        double[]  z = new double[Nx2];
        double[][] d = new double[Nx2][Nx2];
        for (int i = 0; i < Nx2; i++) {
            int j   = i / 2;
            z[i]    = x[j];
            d[0][i] = y[j];
        }
        for (int i = 1; i < Nx2; i++) {
            for (int j = 0; j < Nx2 - i; j++)
            {
                if (i == 1 && j % 2 == 0)
                    d[i][j] = yd[j / 2];
                else
                    d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]);
            }
        }

        // ｎ階差分商
        double[] a = new double[Nx2];
        for (int j = 0; j < Nx2; j++)
            a[j] = d[j][0];

        // 0.5刻みで 与えられていない値を補間
        for (int i = 0; i <= 18; i++) {
            double d1 = i * 0.5 - 4.5;
            double d2 = f(d1);
            double d3 = hermite(d1, z, a);

            // 元の関数と比較
            System.out.println(String.format("%5.2f\t%8.5f\t%8.5f\t%8.5f", d1, d2, d3, d2 - d3));
        }
    }

    // 元の関数
    private static double f(double x) {
        return x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2);
    }
    // 導関数
    private static double fd(double x) {
        return 1 - Math.pow(x,2) / 2 + Math.pow(x,4) / (4 * 3 * 2);
    }

    // Hermite (エルミート) 補間
    private static double hermite(double d, double[] z, double[] a) {
        double sum = a[0];
        for (int i = 1; i < Nx2; i++) {
            double prod = a[i];
            for (int j = 0; j < i; j++)
                prod *= (d - z[j]);
            sum += prod;
        }

        return sum;
    }
}

Z:\>javac Java0704.java

Z:\>java Java0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

C++

#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

// データ点の数
const int N   = 7;
const int Nx2 = 14;

// 元の関数
double f(double x);

// 導関数
double fd(double x);

// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[]);

int main()
{
    double x[N], y[N], yd[N];

    // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i]  = d;
        y[i]  = f(d);
        yd[i] = fd(d);
    }

    // 差分商の表を作る
    double z[Nx2];
    double d[Nx2][Nx2];
    for (int i = 0; i < Nx2; i++)
    {
        int j   = i / 2;
        z[i]    = x[j];
        d[0][i] = y[j];
    }

    for (int i = 1; i < Nx2; i++)
    {
        for (int j = 0; j < Nx2 - i; j++)
        {
            if (i == 1 && j % 2 == 0)
                d[i][j] = yd[j / 2];
            else
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]);
        }
    }

    // ｎ階差分商
    double a[Nx2];
    for (int j = 0; j < Nx2; j++)
        a[j] = d[j][0];

    // 0.5刻みで 与えられていない値を補間
    for (int i = 0; i <= 18; i++)
    {
        double d  = i * 0.5 - 4.5;
        double d1 = f(d);
        double d2 = hermite(d, z, a);

        // 元の関数と比較
        cout << setw(5) << fixed << setprecision(2) << d       << '\t';
        cout << setw(8) << fixed << setprecision(5) << d1      << '\t';
        cout << setw(8) << fixed << setprecision(5) << d2      << '\t';
        cout << setw(8) << fixed << setprecision(5) << d1 - d2 << endl;
    }

   return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// 導関数
double fd(double x)
{
    return 1 - pow(x,2) / 2 + pow(x,4) / (4 * 3 * 2);
}

// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[])
{
    double sum = a[0];
    for (int i = 1; i < Nx2; i++)
    {
        double prod = a[i];
        for (int j = 0; j < i; j++)
            prod *= (d - z[j]);
        sum += prod;
    }

    return sum;
}

Z:\>bcc32 -q CP0704.cpp
cp0704.cpp:

Z:\>CP0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000    -0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078    -0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Objective-C

#import <Foundation/Foundation.h>
#import <math.h>

// データ点の数
const int N   = 7;
const int Nx2 = 14;

// 元の関数
double f(double x);

// 導関数
double fd(double x);

// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[]);

int main()
{
    int i, j;
    double x[N], y[N], yd[N];

    // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for (i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i]  = d;
        y[i]  = f(d);
        yd[i] = fd(d);
    }

    // 差分商の表を作る
    double z[Nx2];
    double d[Nx2][Nx2];
    for (i = 0; i < Nx2; i++)
    {
        int j   = i / 2;
        z[i]    = x[j];
        d[0][i] = y[j];
    }

    for (i = 1; i < Nx2; i++)
    {
        for (j = 0; j < Nx2 - i; j++)
        {
            if (i == 1 && j % 2 == 0)
                d[i][j] = yd[j / 2];
            else
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]);
        }
    }

    // ｎ階差分商
    double a[Nx2];
    for (j = 0; j < Nx2; j++)
        a[j] = d[j][0];

    // 0.5刻みで 与えられていない値を補間
    for (i = 0; i <= 18; i++)
    {
        double d  = i * 0.5 - 4.5;
        double d1 = f(d);
        double d2 = hermite(d, z, a);

        // 元の関数と比較
        printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2);
    }

   return 0;
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}
// 導関数
double fd(double x)
{
    return 1 - pow(x,2) / 2 + pow(x,4) / (4 * 3 * 2);
}

// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[])
{
    int i, j;
    double sum = a[0];
    for (i = 1; i < Nx2; i++)
    {
        double prod = a[i];
        for (j = 0; j < i; j++)
            prod *= (d - z[j]);
        sum += prod;
    }

    return sum;
}

xxxxxx@yyyyyy /Z
$ gcc -o OC0704 OC0704.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333    -0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000    -0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167    -0.00000
 1.50    1.00078     1.00078    -0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

D

import std.stdio;
import std.math;

// データ点の数
const int N = 7;
const int Nx2 = 14;

void main(string[] args)
{
    double x[N];
    double y[N];
    double yd[N];

    // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for (int i = 0; i < N; i++)
    {
        double d = i * 1.5 - 4.5;
        x[i]  = d;
        y[i]  = f(d);
        yd[i] = fd(d);
    }

    // 差分商の表を作る
    double z[Nx2];
    double d[Nx2][Nx2];
    for (int i = 0; i < Nx2; i++)
    {
        int j   = i / 2;
        z[i]    = x[j];
        d[0][i] = y[j];
    }

    for (int i = 1; i < Nx2; i++)
    {
        for (int j = 0; j < Nx2 - i; j++)
        {
            if (i == 1 && j % 2 == 0)
                d[i][j] = yd[j / 2];
            else
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j]);
        }
    }

    // ｎ階差分商
    double a[Nx2];
    for (int j = 0; j < Nx2; j++)
        a[j] = d[j][0];

    // 0.5刻みで 与えられていない値を補間
    for (int i = 0; i <= 18; i++)
    {
        double d1 = i * 0.5 - 4.5;
        double d2 = f(d1);
        double d3 = hermite(d1, z, a);

        // 元の関数と比較
        writefln("%5.2f\t%8.5f\t%8.5f\t%8.5f", d1, d2, d3, d2 - d3);
    }
}

// 元の関数
double f(double x)
{
    return x - pow(x,3) / (3 * 2) + pow(x,5) / (5 * 4 * 3 * 2);
}

// 導関数
double fd(double x)
{
    return 1 - pow(x,2) / 2 + pow(x,4) / (4 * 3 * 2);
}

// Hermite (エルミート) 補間
double hermite(double d, double z[], double a[])
{
    double sum = a[0];
    for (int i = 1; i < Nx2; i++)
    {
        double prod = a[i];
        for (int j = 0; j < i; j++)
            prod *= (d - z[j]);
        sum += prod;
    }

    return sum;
}

Z:\>dmd D0704.d

Z:\>D0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099    -0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Go

package main

import "fmt"
import "math"

// データ点の数
const N   =  7
const Nx2 = 14

func main() {
    var x  [N]float64
    var y  [N]float64
    var yd [N]float64

    // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    for i := 0; i < N; i++ {
        var d float64 = float64(i) * 1.5 - 4.5
        x[i]  = d
        y[i]  = f(d)
        yd[i] = fd(d)
    }

    // 差分商の表を作る
    var z[Nx2]      float64
    var d[Nx2][Nx2] float64
    for i := 0; i < Nx2; i++ {
        j      := i / 2
        z[i]    = x[j]
        d[0][i] = y[j]
    }

    for i := 1; i < Nx2; i++ {
        for j := 0; j < Nx2 - i; j++ {
            if (i == 1 && j % 2 == 0) {
                d[i][j] = yd[j / 2]
            } else {
                d[i][j] = (d[i-1][j+1] - d[i-1][j]) / (z[j+i] - z[j])
            }
        }
    }

    // ｎ階差分商
    var a[Nx2] float64
    for j := 0; j < Nx2; j++ {
        a[j] = d[j][0]
    }

    // 0.5刻みで 与えられていない値を補間
    for i := 0; i <= 18; i++ {
        var d  float64 = float64(i) * 0.5 - 4.5
        var d1 float64 = f(d)
        var d2 float64 = hermite(d, z[:], a[:])

        // 元の関数と比較
        fmt.Printf("%5.2f\t%8.5f\t%8.5f\t%8.5f\n", d, d1, d2, d1 - d2)
    }
}

// 元の関数
func f(x float64) float64 {
    return x - math.Pow(x,3) / (3 * 2) + math.Pow(x,5) / (5 * 4 * 3 * 2)
}
// 導関数
func fd(x float64) float64 {
    return 1 - math.Pow(x,2) / 2 + math.Pow(x,4) / (4 * 3 * 2)
}

// Hermite (エルミート) 補間
func hermite(d float64, z []float64, a []float64) float64 {
    var sum float64 = a[0]
    for i := 1; i < Nx2; i++ {
        var prod float64 = a[i]
        for j := 0; j < i; j++ {
            prod *= (d - z[j])
        }
        sum += prod
    }

    return sum
}

Z:\>8g GO0704.go

Z:\>8l -o GO0704.exe GO0704.8

Z:\>GO0704
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

Scala

object Scala0704 {
    // データ点の数 - 1
    val N   =  6
    val Nx2 = 13

    def main(args: Array[String]) {
        // 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
        val x  = (0 to N).map(_ * 1.5 - 4.5)
        val y  = x.map(f)
        val yd = x.map(fd)

        // 差分商の表を作る
        val z = (0 to Nx2).map(_ / 2).map(x(_))
        val d = Array.ofDim[Double](Nx2 + 1, Nx2 + 1)
        for (i <- 0 to Nx2)
            d(0)(i) = y(i / 2)

        for (i <- 1 to Nx2) {
            for (j <- 0 to Nx2 - i)
                if (i == 1 && j % 2 == 0)
                    d(i)(j) = yd(j / 2)
                else
                    d(i)(j) = (d(i-1)(j+1) - d(i-1)(j)) / (z(j+i) - z(j))
        }

        // ｎ階差分商
        val a = (0 to Nx2).map(d(_)(0))

        // 0.5刻みで 与えられていない値を補間
        val d1 = (0 to 18).map(_ * 0.5 - 4.5)
        val d2 = d1.map(f)
        val d3 = d1.map(hermite(_, z, a))

        (d1 zip d2 zip d3).foreach {
            case ((d1, d2), d3) =>
                println("%5.2f\t%8.5f\t%8.5f\t%8.5f".format(d1, d2, d3, d2 - d3))
        }
    }

    // 元の関数
    def f(x:Double) = {
        x - Math.pow(x,3) / (3 * 2) + Math.pow(x,5) / (5 * 4 * 3 * 2)
    }
    // 導関数
    def fd(x:Double) = {
        1 - Math.pow(x,2) / 2 + Math.pow(x,4) / (4 * 3 * 2)
    }

    // Hermite (エルミート) 補間
    def hermite(d:Double, z:IndexedSeq[Double], a:IndexedSeq[Double]) = {
        var sum_list = List(a(0))
        for (i <- 1 to Nx2) {
            var prod_list = List(a(i))
            for (j <- 0 to i - 1) {
                prod_list = (d - z(j))::prod_list
            }
            sum_list = (prod_list.product)::sum_list
        }
        sum_list.sum
    }
}

Z:\>scala Scala0704.scala
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943    -0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

F#

module Fs0704

open System

// データ点の数 - 1
let N   =  6
let Nx2 = 13

// 元の関数
let f (x:double):double =
    x - Math.Pow(x,3.0) / (float (3 * 2)) + Math.Pow(x,5.0) / (float (5 * 4 * 3 * 2))
// 導関数
let fd (x:double):double =
    1.0 - Math.Pow(x,2.0) / (float 2) + Math.Pow(x,4.0) / (float (4 * 3 * 2))

// Hermite (エルミート) 補間
let hermite(d:double) (z:double list) (a:double list) =
    let mutable sum_list = [a.[0]]
    for i in [1..Nx2] do
        let mutable prod_list = [a.[i]]
        for j in [0..i-1] do
            prod_list <- (d - z.[j])::prod_list
        sum_list <- (prod_list |> List.reduce(*))::sum_list
    sum_list |> List.sum

// 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
let x  = [0..N] |> List.map(fun i -> (float i) * 1.5 - 4.5)
let y  = x |> List.map(f)
let yd = x |> List.map(fd)

// 差分商の表を作る
let z =  [0..Nx2] |> List.map(fun i -> i / 2) |> List.map(fun j -> x.[j])
let d = Array2D.zeroCreate<double> (Nx2+1) (Nx2+1)
for i in [0..Nx2] do
    d.[0,i] <- y.[i / 2]

for i in [1..Nx2] do
    for j in [0..Nx2-i] do
        if (i = 1 && j % 2 = 0) then
            d.[i,j] <- yd.[j / 2]
        else
            d.[i,j] <- (d.[i-1,j+1] - d.[i-1,j]) / (z.[j+i] - z.[j])

// ｎ階差分商
let a = [0..Nx2] |> List.map(fun i -> d.[i,0])

// 0.5刻みで 与えられていない値を補間
let d1 = [0..18] |> List.map(fun i -> (float i) * 0.5 - 4.5)
let d2 = d1 |> List.map(f)
let d3 = d1 |> List.map(fun d -> (hermite d z a))

(List.zip (List.zip d1 d2) d3)
|> List.iter (fun ((d1, d2), d3) -> printfn "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (d2 - d3))

exit 0

Z:\>fsi  --nologo --quiet Fs0704.fs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000     0.00000     0.00000
 0.50    0.47943     0.47943     0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

Clojure

(def N    7)
(def Nx2 14)

; 元の関数
(defn f[x]
    (+ (- x (/ (Math/pow x 3.0) (* 3 2))) (/ (Math/pow x 5.0) (* 5 (* 4 (* 3 2))))))

; 導関数
(defn fd[x]
    (+ (- 1.0 (/ (Math/pow x 2.0) 2)) (/ (Math/pow x 4.0) (* 4 (* 3 2)))))

; Hermite (エルミート) 補間
(defn hermite [d z a]
    (def sum_list (list (nth a 0)))
    (doseq [i (range 1 Nx2)]
        (def prod_list (list (nth a i)))
        (doseq [j (range 0 i)]
            (def prod_list (cons (- d (nth z j)) prod_list)))
        (def w (reduce * prod_list))
        (def sum_list (cons w sum_list)))
    (reduce + sum_list))

; 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
(def x  (map #(- (* % 1.5) 4.5) (range 0 N)))
(def y  (map #(f  %) x))
(def yd (map #(fd %) x))

; 差分商の表を作る
(def z (map #(nth x %) (map #(quot % 2) (range 0 Nx2))))
(def w (map #(nth y %) (map #(quot % 2) (range 0 Nx2))))
(def d (cons w nil))
(doseq [i (range 1 Nx2)]
    (def w (nth d 0))
    (def t (list))
    (doseq [j (range 0 (- Nx2 i))]
        (def t
            (cons
                (if (and (= i 1) (= (rem j 2) 0))
                    (nth yd (quot j 2))
                    (/ (- (nth w (+ j 1)) (nth w j))
                       (- (nth z (+ j i)) (nth z j)))
                )
                t
            )
        )
    )
    (def d (cons (reverse t) d))
)
(def d (reverse d))

; ｎ階差分商
(def a (map #(nth (nth d %) 0) (range 0 Nx2)))

; 0.5刻みで 与えられていない値を補間
(def d1 (map #(- (* % 0.5) 4.5) (range 0 19)))
(def d2 (map #(f %) d1))
(def d3 (map #(hermite % z a) d1))

(doseq [d (map list d1 d2 d3)]
    (def d1 (nth d 0))
    (def d2 (nth d 1))
    (def d3 (nth d 2))
    (println (format "%5.2f\t%8.5f\t%8.5f\t%8.5f" d1 d2 d3 (- d2 d3))))

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0704.clj
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943    -0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333     0.00000
 2.50    0.70964     0.70964    -0.00000
 3.00    0.52500     0.52500    -0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667     0.00000
 4.50    4.68984     4.68984     0.00000

TOP

Haskell

import Text.Printf
import Control.Monad

-- データ点の数 - 1
n   =  6 :: Int
nx2 = 13 :: Int

-- 元の関数
f::Double->Double
f x = x - (x ^ 3) / (fromIntegral (3 * 2)) + (x ^ 5) / (fromIntegral (5 * 4 * 3 * 2))

-- 導関数
fd::Double->Double
fd x = 1.0 - (x ^ 2) / (fromIntegral 2) + (x ^ 4) / (fromIntegral (4 * 3 * 2))

-- Hermite (エルミート) 補間
hermite::Double->[Double]->[Double]->Double
hermite d z a =
    let sum_list = map(\i -> do
        let prod_list = map(\j -> do
            d - z!!j
            ) $ [0..(i-1)::Int]
        product $ a!!i : prod_list
        ) [1..nx2::Int]
    in
        sum $ a!!0 : sum_list

-- 差分商の表を作る
make_table::[Double]->[Double]->[Double]->[Double]->Int->[Double]
make_table yd z d a i =
    let
        w = map(\j -> do
            if (i == 1 && (rem j 2) == 0) then
                yd!!(div j 2)
            else
                ((d!!(j+1) - d!!j) / (z!!(j+i) - z!!j))
            ) $ [0..(nx2-i)::Int]
        t = w!!0:a
    in
        -- ｎ階差分商
        if i == nx2 then t
                    else (make_table yd z w t (i + 1))

main = do
    -- 1.5刻みで -4.5～4.5 まで, ７点だけ値をセット
    let x  = map(\i -> (fromIntegral i) * 1.5 - 4.5) [0..n]
    let y  = map(\i -> f(i)) x
    let yd = map(\i -> fd(i)) x

    -- 差分商の表を作る
    let z = map(\i -> x!!(div i 2)) [0..nx2]
    let w = map(\i -> y!!(div i 2)) [0..nx2]
    let a = reverse (make_table yd z w [w!!0] 1)

    -- 0.5刻みで 与えられていない値を補間
    let d1 = map(\i -> (fromIntegral i) * 0.5 - 4.5) [0..18]
    let d2 = map(\i -> (f i)) d1
    let d3 = map(\i -> (hermite i z a)) d1

    forM_ (zip (zip d1 d2) d3) $ \((d1, d2), d3) -> do
        printf "%5.2f\t%8.5f\t%8.5f\t%8.5f\n" d1 d2 d3 (d2 - d3)

Z:\>runghc Hs0704.hs
-4.50   -4.68984    -4.68984     0.00000
-4.00   -1.86667    -1.86667     0.00000
-3.50   -0.73099    -0.73099     0.00000
-3.00   -0.52500    -0.52500     0.00000
-2.50   -0.70964    -0.70964     0.00000
-2.00   -0.93333    -0.93333     0.00000
-1.50   -1.00078    -1.00078     0.00000
-1.00   -0.84167    -0.84167     0.00000
-0.50   -0.47943    -0.47943     0.00000
 0.00    0.00000    -0.00000     0.00000
 0.50    0.47943     0.47943    -0.00000
 1.00    0.84167     0.84167     0.00000
 1.50    1.00078     1.00078     0.00000
 2.00    0.93333     0.93333    -0.00000
 2.50    0.70964     0.70964     0.00000
 3.00    0.52500     0.52500     0.00000
 3.50    0.73099     0.73099     0.00000
 4.00    1.86667     1.86667    -0.00000
 4.50    4.68984     4.68984    -0.00000

TOP

VBScript
JScript
PowerShell
Perl
PHP
Python
Ruby
Groovy

Pascal
Ada
VB.NET

C#
Java
C++
Objective-C
D
Go

Scala
F#
Clojure
Haskell