ロンバーグ積分

home > さまざまな言語で数値計算 > 数値積分 >

ロンバーグ積分

台形則を使って近似を行ったあと，Richardsonの補外法を用いて近似をより正確にする.

まず, $ T_{1,1} $ から $ T_{n,1} $ を, 台形則を使って求める.
次に, $ T_{2,2} $ から $ T_{n,n} $ まで, 以下の計算を行う.

例題として, $ 4/(1+x^2) $ を $0$ から $1$ まで積分し $\pi$ を求める.

VBScript

Option Explicit

Const PI = 3.14159265359
Const a = 0
Const b = 1

Dim t(6, 6)

'台形則で積分
Dim n: n = 2
Dim i, j
For i = 1 To 6
    Dim h: h = (b - a) / n
    Dim s: s = 0
    Dim x: x = a
    For j = 1 To n - 1
        x = x + h
        s = s + f(x)
    Next
    '結果を保存
    t(i,1) = h * ((f(a) + f(b)) / 2 + s)
    n = n * 2
Next

'Richardsonの補外法
n = 4
For j = 2 To 6
    For i = j To 6
        t(i,j) = t(i, j - 1) + (t(i, j - 1) - t(i - 1, j - 1)) / (n - 1)
        If i = j Then
            '結果を π と比較
            WScript.StdOut.Write Right(Space(2)  & j,                                        2) & " : "
            WScript.StdOut.Write Right(Space(13) & FormatNumber(t(i,j),      10, -1, 0, 0), 13) & ", "
            WScript.StdOut.Write Right(Space(13) & FormatNumber(t(i,j) - PI, 10, -1, 0, 0), 13) & vbNewLine
        End If
    Next
    n = n * 4
Next

Private Function f(x)
    f = 4 / (1 + x * x)
End Function

Z:\>cscript //nologo Z:\0604.vbs
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

JScript

var a = 0
var b = 1

var t = []

// 台形則で積分
var n = 2;
for (var i = 1; i <= 6; i++)
{
    var h = (b - a) / n
    var s = 0
    var x = a
    for (var j = 1; j <= n - 1; j++)
    {
        x += h
        s += f(x)
    }
    // 結果を保存
    t[i]    = []
    t[i][1] = h * ((f(a) + f(b)) / 2 + s)
    n *= 2
}

// Richardsonの補外法
n = 4
for (var j = 2; j <= 6; j++)
{
    for (var i = j; i <= 6; i++)
    {
        t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1)
        if (i == j)
        {
            // 結果を π と比較
            WScript.StdOut.Write(("  "            + j                              ).slice( -2) + " : ");
            WScript.StdOut.Write(("             " + t[i][j].toFixed(10)            ).slice(-13) + ", ");
            WScript.StdOut.Write(("             " + (t[i][j] - Math.PI).toFixed(10)).slice(-13) + "\n" );
        }
    }
    n *= 4
}

function f(x)
{
    return 4 / (1 + x * x)
}

Z:\>cscript //nologo Z:\0604.js
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000000
 6 :  3.1415926536, -0.0000000000

PowerShell

function f($x)
{
    4 / (1 + $x * $x)
}

$a = 0
$b = 1

$t = New-Object "double[,]" 7,7

# 台形則で積分
$n = 2;
foreach ($i in 1..6)
{
    $h = ($b - $a) / $n
    $s = 0
    $x = $a
    foreach ($j in 1..($n - 1))
    {
        $x += $h
        $s += (f $x)
    }
    # 結果を保存
    $t[$i,1] = $h * (((f $a) + (f $b)) / 2 + $s)
    $n *= 2
}

# Richardsonの補外法
$n = 4
foreach ($j in 2..6)
{
    foreach ($i in $j..6)
    {
        $t1 = $t[ $i     ,($j - 1)]
        $t2 = $t[($i - 1),($j - 1)]
        $t[$i,$j] = $t1 + ($t1 - $t2) / ($n - 1)
        if ($i -eq $j)
        {
            # 結果を π と比較
            $s = $t[$i,$j]
            Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, ($s - [Math]::PI)))
        }
    }
    $n *= 4;
}

Z:\>powershell -file Z:\0604.ps1
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536,  0.0000000000

Perl

use Math::Trig 'pi';

my $a = 0;
my $b = 1;

my @t = ();

# 台形則で積分
my $n = 2;
for $i (1..6)
{
    my $h = ($b - $a) / $n;
    my $s = 0;
    my $x = $a;
    for $j (1..($n - 1))
    {
        $x += $h;
        $s += f($x);
    }
    # 結果を保存
    $t[$i][1] = $h * ((f($a) + f($b)) / 2 + $s);
    $n *= 2;
}

# Richardsonの補外法
$n = 4;
for $j (2..6)
{
    for $i ($j..6)
    {
        $t[$i][$j] = $t[$i][$j - 1] + ($t[$i][$j - 1] - $t[$i - 1][$j - 1]) / ($n - 1);
        if ($i == $j)
        {
            # 結果を π と比較
            printf("%2d : %13.10f, %13.10f\n", $j, $t[$i][$j], $t[$i][$j] - pi);
        }
    }
    $n *= 4;
}

sub f
{
    my ($x) = @_;
    4 / (1 + $x * $x);
}

Z:\>perl Z:\0604.pl
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

PHP

<?php
$a = 0;
$b = 1;

$t = array();

# 台形則で積分
$n = 2;
foreach (range(1, 6) as $i)
{
    $h = ($b - $a) / $n;
    $s = 0;
    $x = $a;
    foreach (range(1, ($n - 1)) as $j)
    {
        $x += $h;
        $s += f($x);
    }
    # 結果を保存
    $t[$i][1] = $h * ((f($a) + f($b)) / 2 + $s);
    $n *= 2;
}

# Richardsonの補外法
$n = 4;
foreach (range(2, 6) as $j)
{
    foreach (range($j, 6) as $i)
    {
        $t[$i][$j] = $t[$i][$j - 1] + ($t[$i][$j - 1] - $t[$i - 1][$j - 1]) / ($n - 1);
        if ($i == $j)
        {
            # 結果を π と比較
            printf("%2d : %13.10f, %13.10f\n", $j, $t[$i][$j], $t[$i][$j] - M_PI);
        }
    }
    $n *= 4;
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>

Z:\>php Z:\0604.php
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

Python

# coding: Shift_JIS

import math

def f(x):
    return 4 / (1 + x * x)

a = 0
b = 1

t = [[0 for j in range(7)] for i in range(7)]

# 台形則で積分
n = 2
for i in range(1, 7):
    h = (b - a) / float(n)
    s = 0
    x = a
    for j in range(1, n):
        x += h
        s += f(x)
    # 結果を保存
    t[i][1] = h * ((f(a) + f(b)) / 2 + s)
    n *= 2

# Richardsonの補外法
n = 4
for j in range(2, 7):
    for i in range(j, 7):
        t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1)
        if (i == j):
            # 結果を π と比較
            print "%2d : %13.10f, %13.10f" % (j, t[i][j], t[i][j] - math.pi)
    n *= 4

Z:\>python Z:\0604.py
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

Ruby

def f(x)
    4 / (1 + x * x)
end

a = 0
b = 1

t = Array.new(7) { Array.new(7) }

# 台形則で積分
n = 2
(1..6).each do |i|
    h = (b - a) / n.to_f
    s = 0
    x = a
    (1..(n - 1)).each do |j|
        x += h
        s += f(x)
    end
    # 結果を保存
    t[i][1] = h * ((f(a) + f(b)) / 2 + s)
    n *= 2
end

# Richardsonの補外法
n = 4
(2..6).each do |j|
    (j..6).each do |i|
        t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1)
        if i == j
            # 結果を π と比較
            printf("%2d : %13.10f, %13.10f\n", j, t[i][j], t[i][j] - Math::PI)
        end
    end
    n *= 4
end

Z:\>ruby Z:\0604.rb
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

Groovy

Pascal

program Pas0604(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f(x:Double):Double;
begin
    result := 4 / (1 + x * x);
end;

const
    a:Double = 0;
    b:Double = 1;
var
    n, i, j:Integer;
    h, s, x:Double;
    t:array [1..6, 1..6] of Double;
begin
    // 台形則で積分
    n := 2;
    for i := 1 to 6 do
    begin
        h := (b - a) / n;
        s := 0;
        x := a;
        for j := 1 to n - 1 do
        begin
            x := x + h;
            s := s + f(x);
        end;
        // 結果を保存
        t[i,1] := h * ((f(a) + f(b)) / 2 + s);
        n      := n * 2;
    end;

    // Richardsonの補外法
    n := 4;
    for j := 2 to 6 do
    begin
        for i := j to 6 do
        begin
            t[i,j] := t[i, j - 1] + (t[i, j - 1] - t[i - 1, j - 1]) / (n - 1);
            if i = j then
            begin
                // 結果を π と比較
                writeln(format('%2d : %13.10f, %13.10f', [j, t[i,j], t[i,j] - PI]));
            end;
        end;
        n := n * 4;
    end;
end.

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

Z:\>Pas0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536,  0.0000000000

Ada

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

procedure Ada0604 is
    a : Constant Long_Float := 0.0;
    b : Constant Long_Float := 1.0;

    h, s, x : Long_Float;
    n : Integer;
    t : array (1..6, 1..6) of Long_Float;

    function f(x:Long_Float) return Long_Float is
    begin
        return 4.0 / (1.0 + x * x);
    end f;

begin
    -- 台形則で積分
    n := 2;
    for i in 1..6 loop
        h := (b - a) / Long_Float(n);
        s := 0.0;
        x := a;
        for j in 1..(n - 1) loop
            x := x + h;
            s := s + f(x);
        end loop;
        -- 結果を保存
        t(i,1) := h * ((f(a) + f(b)) / 2.0 + s);
        n := n * 2;
    end loop;

    -- Richardsonの補外法
    n := 4;
    for j in 2..6 loop
        for i in j..6 loop
            t(i,j) := t(i, j - 1) + (t(i, j - 1) - t(i - 1, j - 1)) / Long_Float(n - 1);
            if i = j then
                -- 結果を π と比較
                Put(j,                        Width=> 2);
                Put(" : ");
                Put(t(i,j) ,                  Fore=>3, Aft=>10, Exp=>0);
                Put(", ");
                Put(t(i,j) - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0);
                New_Line;
            end if;
        end loop;
        n := n * 4;
    end loop;
end Ada0604;

xxxxxx@yyyyyy /Z
$ gnatmake Ada0604.adb

xxxxxx@yyyyyy /Z
$ Ada0604
 2 :   3.1415686275,  -0.0000240261
 3 :   3.1415940941,   0.0000014405
 4 :   3.1415926384,  -0.0000000152
 5 :   3.1415926536,   0.0000000001
 6 :   3.1415926536,  -0.0000000000

VB.NET

Module VB0604
    Public Sub Main()
        Const a As Double = 0
        Const b As Double = 1

        Dim t(6, 6) As Double

        '台形則で積分
        Dim n As Integer = 2
        For i As Integer = 1 To 6
            Dim h As Double = (b - a) / n
            Dim s As Double = 0
            Dim x As Double = a
            For j As Integer = 1 To n - 1
                x += h
                s += f(x)
            Next
            '結果を保存
            t(i,1) = h * ((f(a) + f(b)) / 2 + s)
            n *= 2
        Next

        'Richardsonの補外法
        n = 4
        For j As Integer = 2 To 6
            For i As Integer = j To 6
                t(i,j) = t(i, j - 1) + (t(i, j - 1) - t(i - 1, j - 1)) / (n - 1)
                If i = j Then
                    '結果を π と比較
                    Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, t(i,j), t(i,j) - Math.PI))
                End If
            Next
            n *= 4
        Next
    End Sub

    Private Function f(ByVal x As Double) As Double
        Return 4 / (1 + x * x)
    End Function
End Module

Z:\>vbc -nologo VB0604.vb

Z:\>VB0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536,  0.0000000000

C#

using System;

public class CS0604
{
    private static double f(double x)
    {
        return 4 / (1 + x * x);
    }

    public static void Main()
    {
        const double a = 0;
        const double b = 1;

        double[,] t = new double[7,7];

        // 台形則で積分
        int n = 2;
        for (int i = 1; i <= 6; i++)
        {
            double h = (b - a) / n;
            double s = 0;
            double x = a;
            for (int j = 1; j <= n - 1; j++)
            {
                x += h;
                s += f(x);
            }
            // 結果を保存
            t[i,1] = h * ((f(a) + f(b)) / 2 + s);
            n *= 2;
        }

        // Richardsonの補外法
        n = 4;
        for (int j = 2; j <= 6; j++)
        {
            for (int i = j; i <= 6; i++)
            {
                t[i,j] = t[i, j - 1] + (t[i, j - 1] - t[i - 1, j - 1]) / (n - 1);
                if (i == j)
                {
                    // 結果を π と比較
                    Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, t[i,j], t[i,j] - Math.PI));
                }
            }
            n *= 4;
        }
    }
}

Z:\>csc -nologo CS0604.cs

Z:\>CS0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536,  0.0000000000

Java

public class Java0604 {

    private static double f(double x) {
        return 4 / (1 + x * x);
    }

    public static void main(String []args) {
        final double a = 0;
        final double b = 1;

        double[][] t = new double[7][7];

        // 台形則で積分
        int n = 2;
        for (int i = 1; i <= 6; i++) {
            double h = (b - a) / n;
            double s = 0;
            double x = a;
            for (int j = 1; j <= n - 1; j++) {
                x += h;
                s += f(x);
            }
            // 結果を保存
            t[i][1] = h * ((f(a) + f(b)) / 2 + s);
            n *= 2;
        }

        // Richardsonの補外法
        n = 4;
        for (int j = 2; j <= 6; j++) {
            for (int i = j; i <= 6; i++) {
                t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1);
                if (i == j) {
                    // 結果を π と比較
                    System.out.println(String.format("%2d : %13.10f, %13.10f", j, t[i][j], t[i][j] - Math.PI));
                }
            }
            n *= 4;
        }
    }
}

Z:\>javac Java0604.java

Z:\>java Java0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

C++

#include <iostream>
#include <iomanip>
#include <math>

using namespace std;

double f(double x)
{
    return 4 / (1 + x * x);
}

int main()
{
    const double a = 0;
    const double b = 1;

    double t[7][7];

    // 台形則で積分
    int n = 2;
    for (int i = 1; i <= 6; i++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a;
        for (int j = 1; j <= n - 1; j++)
        {
            x += h;
            s += f(x);
        }
        // 結果を保存
        t[i][1] = h * ((f(a) + f(b)) / 2 + s);
        n *= 2;
    }

    // Richardsonの補外法
    n = 4;
    for (int j = 2; j <= 6; j++)
    {
        for (int i = j; i <= 6; i++)
        {
            t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1);
            if (i == j)
            {
                // 結果を π と比較
                cout << setw(2)  << j << ":";
                cout << setw(13) << fixed << setprecision(10) << t[i][j]        << ", ";
                cout << setw(13) << fixed << setprecision(10) << t[i][j] - M_PI << endl;
            }
        }
        n *= 4;
    }

    return 0;
}

Z:\>bcc32 -q CP0604.cpp
cp0604.cpp:

Z:\>CP0604
 2: 3.1415686275, -0.0000240261
 3: 3.1415940941,  0.0000014405
 4: 3.1415926384, -0.0000000152
 5: 3.1415926536,  0.0000000001
 6: 3.1415926536, -0.0000000000

Objective-C

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

double f(double x)
{
    return 4 / (1 + x * x);
}

int main()
{
    const double a = 0;
    const double b = 1;

    double t[7][7];

    // 台形則で積分
    int n = 2;
    int i, j;
    for (i = 1; i <= 6; i++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a;
        for (j = 1; j <= n - 1; j++)
        {
            x += h;
            s += f(x);
        }
        // 結果を保存
        t[i][1] = h * ((f(a) + f(b)) / 2 + s);
        n *= 2;
    }

    // Richardsonの補外法
    n = 4;
    for (j = 2; j <= 6; j++)
    {
        for (i = j; i <= 6; i++)
        {
            t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1);
            if (i == j)
            {
                // 結果を π と比較
                printf("%2d : %13.10f, %13.10f\n", j, t[i][j], t[i][j] - M_PI);
            }
        }
        n *= 4;
    }

    return 0;
}

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

xxxxxx@yyyyyy /Z
$ OC0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

D

import std.stdio;
import std.math;

void main(string[] args)
{
    const double a = 0;
    const double b = 1;

    double t[7][7];

    // 台形則で積分
    int n = 2;
    for (int i = 1; i <= 6; i++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a;
        for (int j = 1; j <= n - 1; j++)
        {
            x += h;
            s += f(x);
        }
        // 結果を保存
        t[i][1] = h * ((f(a) + f(b)) / 2 + s);
        n *= 2;
    }

    // Richardsonの補外法
    n = 4;
    for (int j = 2; j <= 6; j++)
    {
        for (int i = j; i <= 6; i++)
        {
            t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / (n - 1);
            if (i == j)
            {
                // 結果を π と比較
                writefln("%2d : %13.10f, %13.10f", j, t[i][j], t[i][j] - PI);
            }
        }
        n *= 4;
    }
}

double f(double x)
{
    return 4 / (1 + x * x);
}

Z:\>dmd D0604.d

Z:\>D0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000151
 5 :  3.1415926536,  0.0000000000
 6 :  3.1415926536, -0.0000000000

Go

package main

import "fmt"
import "math"

func main() {
    const a = 0
    const b = 1

    var t[7][7] float64

    // 台形則で積分
    var n int = 2
    for i := 1; i <= 6; i++ {
        var h float64 = (b - a) / float64(n)
        var s float64 = 0
        var x float64 = a
        for j := 1; j <= n - 1; j++ {
            x += h
            s += f(x)
        }
        // 結果を保存
        t[i][1] = h * ((f(a) + f(b)) / 2 + s)
        n *= 2
    }

    // Richardsonの補外法
    n = 4
    for j := 2; j <= 6; j++ {
        for i := j; i <= 6; i++ {
            t[i][j] = t[i][j - 1] + (t[i][j - 1] - t[i - 1][j - 1]) / float64(n - 1)
            if i == j {
                // 結果を π と比較
                fmt.Printf("%2d : %13.10f, %13.10f\n", j, t[i][j], t[i][j] - math.Pi)
            }
        }
        n *= 4
    }

}

func f(x float64) float64 {
    return (4 / (1 + x * x))
}

Z:\>8g GO0604.go

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

Z:\>GO0604
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

Scala

object Scala0604 {
    def main(args: Array[String]) {
        val a:Double = 0
        val b:Double = 1

        val t = Array.ofDim[Double](7, 7)

        // 台形則で積分
        var n:Int = 2
        for (i <- 1 to 6) {
            var h:Double = (b - a) / n
            var s:Double = 0
            var x:Double = a
            for (j <- 1 to (n - 1)) {
                x += h
                s += f(x)
            }
            // 結果を保存
            t(i)(1) = h * ((f(a) + f(b)) / 2 + s)
            n *= 2
        }

        // Richardsonの補外法
        n = 4
        for (j <- 2 to 6) {
            for (i <- j to 6) {
                t(i)(j) = t(i)(j - 1) + (t(i)(j - 1) - t(i - 1)(j - 1)) / (n - 1)
                if (i == j) {
                    // 結果を π と比較
                    println("%3d : %13.10f, %13.10f".format(j, t(i)(j), t(i)(j) - Math.PI))
                }
            }
            n *= 4
        }
    }

    def f(x:Double):Double = {
        4 / (1 + x * x)
    }
}

Z:\>scala Scala0604.scala
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415940941,  0.0000014405
  4 :  3.1415926384, -0.0000000152
  5 :  3.1415926536,  0.0000000001
  6 :  3.1415926536, -0.0000000000

F#

module Fs0604

open System

let f(x:double):double =
    4.0 / (1.0 + x * x)

let a:double = 0.0
let b:double = 1.0

let t = Array2D.zeroCreate<double> 7 7

// 台形則で積分
let mutable n = 2
for i in [1..6] do
    let         h:double = (b - a) / (double n)
    let mutable s:double = 0.0
    let mutable x:double = a

    for j in [1..(n - 1)] do
        x <- x + h
        s <- s + (f x)

    // 結果を保存
    t.[i,1] <- h * (((f a) + (f b)) / 2.0 + s)
    n <- n * 2

// Richardsonの補外法
n <- 4
for j in [2..6] do
    for i in [j..6] do
        t.[i,j] <- t.[i, j - 1] + (t.[i, j - 1] - t.[i - 1, j - 1]) / (double (n - 1))
        if i = j then
            // 結果を π と比較
            printfn "%3d : %13.10f, %13.10f" j t.[i,j] (t.[i,j] - Math.PI)
    n <- n * 4

exit 0

Z:\>fsi  --nologo --quiet Fs0604.fs
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415940941,  0.0000014405
  4 :  3.1415926384, -0.0000000152
  5 :  3.1415926536,  0.0000000001
  6 :  3.1415926536,  0.0000000000

Clojure

(defn f[x]
    (/ 4 (+ 1 (* x x))))

(def a 0)
(def b 1)

(def t (list))
; 台形則で積分
(doseq [j (range 1 11)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x a)
    (def s 0)
    (doseq [i (range 1 n)]
        (def x (+ x h))
        (def s (+ s (f x))))
    (def t1 (double (* h (+ (/ (+ (f a) (f b)) 2) s))))
    (def t2 (- t1 (. Math PI)))
    ; 結果を保存
    (def t (cons t1 t)))

; Richardsonの補外法
(def t (reverse t))
(doseq [j (range 2 7)]
    (def n (Math/pow 4 (- j 1)))
    (def s (list))
    (doseq [i (range j 7)]
        (def t1 (+ (second t) (/ (- (second t) (first t)) (- n 1))))
        (def t2 (- t1 (. Math PI)))
        ; 結果を π と比較
        (if (= i j)
            (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
        (def t (rest t))
        (def s (cons t1 s)))
    (def t (reverse s)))

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0604.clj
 2 :  3.1415686275, -0.0000240261
 3 :  3.1415940941,  0.0000014405
 4 :  3.1415926384, -0.0000000152
 5 :  3.1415926536,  0.0000000001
 6 :  3.1415926536, -0.0000000000

Haskell

import Text.Printf
import Control.Monad

f::Double->Double
f x = 4 / (1 + (x * x))

-- Richardsonの補外法
rich_sub n []     ss = ss
rich_sub n (x:[]) ss = ss
rich_sub n (x:xs) ss =
    let
        x2 = (head xs)
        s  = x2 + (x2 - x) / (fromIntegral (n - 1))
    in
        rich_sub n xs (s:ss)

richardson n ([])   = 0
richardson n (x:[]) = x
richardson n (x:xs) =
    let
        s = rich_sub n (x:xs) []
    in
        richardson (n * 4) (reverse s)

main = do
    let a = 0.0
    let b = 1.0
    t <- forM ([1..6::Integer]) $ \j -> do
        let n = 2 ^ j
        let h = (b - a) / (fromIntegral n)
        -- 台形則で積分
        let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a) $ [1..(n - 1)]
        let w2 = ((f a) + (f b)) / 2
        let t1 = h * (w1 + w2)
        return t1

    -- Richardsonの補外法で積分して, 結果を π と比較
    forM_ ([2..6::Int]) $ \i -> do
        let s = richardson 4 (take i t)
        printf "%3d : %13.10f, %13.10f\n" i s (s - pi)

Z:\>runghc Hs0604.hs
  2 :  3.1415686275, -0.0000240261
  3 :  3.1415940941,  0.0000014405
  4 :  3.1415926384, -0.0000000152
  5 :  3.1415926536,  0.0000000001
  6 :  3.1415926536, -0.0000000000

VBScript
JScript
PowerShell
Perl
PHP
Python
Ruby
Groovy

Pascal
Ada
VB.NET

C#
Java
C++
Objective-C
D
Go

Scala
F#
Clojure
Haskell