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

さまざまな言語で数値計算

Only Do What Only You Can Do

中点則

関数 $ f(x) $ の $ a, b $ 区間を $n$ 個に分割し, それらを長方形で近似し面積を合計する.

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

VBScript

Option Explicit

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

'中点則で積分
Dim n: n = 2
Dim i, j
For j = 1 To 10
    Dim h: h = (b - a) / n
    Dim s: s = 0
    Dim x: x = a + (h / 2)
    For i = 1 To n
        s = s + f(x)
        x = x + h
    Next
    s = h * s
    n = n * 2

    '結果を π と比較
    WScript.StdOut.Write Right(Space(2)  & j,                                   2) & " : "
    WScript.StdOut.Write Right(Space(13) & FormatNumber(s,      10, -1, 0, 0), 13) & ", "
    WScript.StdOut.Write Right(Space(13) & FormatNumber(s - PI, 10, -1, 0, 0), 13) & vbNewLine
Next

Private Function f(x)
    f = 4 / (1 + x * x)
End Function
Z:\>cscript //nologo Z:\0602.vbs
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

JScript

var a = 0
var b = 1

// 中点則で積分
var n = 2;
for (var j = 1; j <= 10; j++)
{
    var h = (b - a) / n
    var s = 0
    var x = a + (h / 2)
    for (var i = 1; i <= n; i++)
    {
        s += f(x)
        x += h
    }
    s *= h
    n *= 2

    // 結果を π と比較
    WScript.StdOut.Write(("  "            + j                        ).slice( -2) + " : ");
    WScript.StdOut.Write(("             " + s.toFixed(10)            ).slice(-13) + ", ");
    WScript.StdOut.Write(("             " + (s - Math.PI).toFixed(10)).slice(-13) + "\n" );
}

function f(x)
{
    return 4 / (1 + x * x)
}
Z:\>cscript //nologo Z:\0602.js
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

PowerShell

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

$a = 0
$b = 1

# 中点則で積分
$n = 2;
foreach ($j in 1..10)
{
    $h = ($b - $a) / $n
    $s = 0
    $x = $a + ($h / 2)
    foreach ($i in 1..$n)
    {
        $s += (f $x)
        $x += $h
    }
    $s *= $h
    $n *= 2

    # 結果を π と比較
    Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI))
}
Z:\>powershell -file Z:\0602.ps1
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Perl

use Math::Trig 'pi';

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

# 中点則で積分
my $n = 2;
for $j (1..10)
{
    my $h = ($b - $a) / $n;
    my $s = 0;
    my $x = $a + ($h / 2);
    for $i (1..$n)
    {
        $s += f($x);
        $x += $h;
    }
    $s *= $h;
    $n *= 2;

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - pi);
}

sub f
{
    my ($x) = @_;
    4 / (1 + $x * $x);
}
Z:\>perl Z:\0602.pl
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

PHP

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

# 中点則で積分
$n = 2;
foreach (range(1, 10) as $j)
{
    $h = ($b - $a) / $n;
    $s = 0;
    $x = $a + ($h / 2);
    foreach (range(1, $n) as $i)
    {
        $s += f($x);
        $x += $h;
    }
    $s *= $h;
    $n *= 2;

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - M_PI);
}

function f($x)
{
    return 4 / (1 + $x * $x);
}
?>
Z:\>php Z:\0602.php
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Python

# coding: Shift_JIS

import math

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

a = 0
b = 1

# 中点則で積分
n = 2
for j in range(1, 11):
    h = (b - a) / float(n)
    s = 0
    x = a + (h / 2)
    for i in range(1, (n + 1)):
        s += f(x)
        x += h
    s *= h
    n *= 2

    # 結果を π と比較
    print "%2d : %13.10f, %13.10f" % (j, s, s - math.pi)
Z:\>python Z:\0602.py
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Ruby

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

a = 0
b = 1

# 中点則で積分
n = 2
(1..10).each do |j|
    h = (b - a) / n.to_f
    s = 0
    x = a + (h / 2)
    (1..n).each do |i|
        s += f(x)
        x += h
    end
    s *= h
    n *= 2

    # 結果を π と比較
    printf("%2d : %13.10f, %13.10f\n", j, s, s - Math::PI)
end
Z:\>ruby Z:\0602.rb
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Groovy

Pascal

program Pas0602(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;
begin
    // 中点則で積分
    n := 2;
    for j := 1 to 10 do
    begin
        h := (b - a) / n;
        s := 0;
        x := a + (h / 2);
        for i := 1 to n do
        begin
            s := s + f(x);
            x := x + h;
        end;
        s := h * s;
        n := n * 2;

        // 結果を π と比較
        writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI]));
    end;
end.
Z:\>fpc -v0 -l- Pas0602.pp

Z:\>Pas0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

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 Ada0602 is
    a : Constant Long_Float := 0.0;
    b : Constant Long_Float := 1.0;

    h, s, x : Long_Float;
    n : Integer;

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

begin
    -- 中点則で積分
    n := 2;
    for j in 1..10 loop
        h := (b - a) / Long_Float(n);
        s := 0.0;
        x := a + (h / 2.0);
        for i in 1..n loop
            s := s + f(x);
            x := x + h;
        end loop;
        s := h * s;
        n := n * 2;

        -- 結果を π と比較
        Put(j,                   Width=> 2);
        Put(" : ");
        Put(s ,                  Fore=>3, Aft=>10, Exp=>0);
        Put(", ");
        Put(s - Ada.Numerics.Pi, Fore=>3, Aft=>10, Exp=>0);
        New_Line;
    end loop;
end Ada0602;
xxxxxx@yyyyyy /Z
$ gnatmake Ada0602.adb

xxxxxx@yyyyyy /Z
$ Ada0602
 1 :   3.1623529412,   0.0207602876
 2 :   3.1468005184,   0.0052078648
 3 :   3.1428947296,   0.0013020760
 4 :   3.1419181743,   0.0003255207
 5 :   3.1416740338,   0.0000813802
 6 :   3.1416129986,   0.0000203451
 7 :   3.1415977399,   0.0000050863
 8 :   3.1415939252,   0.0000012716
 9 :   3.1415929715,   0.0000003179
10 :   3.1415927331,   0.0000000795

VB.NET

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

        '中点則で積分
        Dim n As Integer = 2
        For j As Integer = 1 To 10
            Dim h As Double = (b - a) / n
            Dim s As Double = 0
            Dim x As Double = a + (h / 2)
            For i As Integer = 1 To n
                s += f(x)
                x += h
            Next
            s = h * s
            n *= 2

            '結果を π と比較
            Console.WriteLine(String.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI))
        Next
    End Sub

    Private Function f(ByVal x As Double) As Double
        Return 4 / (1 + x * x)
    End Function
End Module
Z:\>vbc -nologo VB0602.vb

Z:\>VB0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

C#

using System;

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

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

        // 中点則で積分
        int n = 2;
        for (int j = 1; j <= 10; j++)
        {
            double h = (b - a) / n;
            double s = 0;
            double x = a + (h / 2);
            for (int i = 1; i <= n; i++)
            {
                s += f(x);
                x += h;
            }
            s *= h;
            n *= 2;

            // 結果を π と比較
            Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI));
        }
    }
}
Z:\>csc -nologo CS0602.cs

Z:\>CS0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Java

public class Java0602 {

    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;

        // 中点則で積分
        int n = 2;
        for (int j = 1; j <= 10; j++) {
            double h = (b - a) / n;
            double s = 0;
            double x = a + (h / 2);
            for (int i = 1; i <= n; i++) {
                s += f(x);
                x += h;
            }
            s *= h;
            n *= 2;

            // 結果を π と比較
            System.out.println(String.format("%2d : %13.10f, %13.10f", j, s, s - Math.PI));
        }
    }
}
Z:\>javac Java0602.java

Z:\>java Java0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

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;

    // 中点則で積分
    int n = 2;
    for (int j = 1; j <= 10; j++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a + (h / 2);
        for (int i = 1; i <= n; i++)
        {
            s += f(x);
            x += h;
        }
        s *= h;
        n *= 2;

        // 結果を π と比較
        cout << setw(2)  << j << ":";
        cout << setw(13) << fixed << setprecision(10) << s        << ", ";
        cout << setw(13) << fixed << setprecision(10) << s - M_PI << endl;
    }

    return 0;
}
Z:\>bcc32 -q CP0602.cpp
cp0602.cpp:

Z:\>CP0602
 1: 3.1623529412,  0.0207602876
 2: 3.1468005184,  0.0052078648
 3: 3.1428947296,  0.0013020760
 4: 3.1419181743,  0.0003255207
 5: 3.1416740338,  0.0000813802
 6: 3.1416129986,  0.0000203451
 7: 3.1415977399,  0.0000050863
 8: 3.1415939252,  0.0000012716
 9: 3.1415929715,  0.0000003179
10: 3.1415927331,  0.0000000795

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;

    // 中点則で積分
    int n = 2;
    int i, j;
    for (j = 1; j <= 10; j++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a + (h / 2);
        for (i = 1; i <= n; i++)
        {
            s += f(x);
            x += h;
        }
        s *= h;
        n *= 2;

        // 結果を π と比較
        printf("%2d : %13.10f, %13.10f\n", j, s, s - M_PI);
    }

    return 0;
}
xxxxxx@yyyyyy /Z
$ gcc -o OC0602 OC0602.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS

xxxxxx@yyyyyy /Z
$ OC0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

D

import std.stdio;
import std.math;

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

    // 中点則で積分
    int n = 2;
    for (int j = 1; j <= 10; j++)
    {
        double h = (b - a) / n;
        double s = 0;
        double x = a + (h / 2);
        for (int i = 1; i <= n; i++)
        {
            s += f(x);
            x += h;
        }
        s *= h;
        n *= 2;

        // 結果を π と比較
        writefln("%2d : %13.10f, %13.10f", j, s, s - PI);
    }
}

double f(double x)
{
    return 4 / (1 + x * x);
}
Z:\>dmd D0602.d

Z:\>D0602
 1 :  3.1623529412,  0.0207602875
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203450
 7 :  3.1415977399,  0.0000050862
 8 :  3.1415939252,  0.0000012715
 9 :  3.1415929715,  0.0000003178
10 :  3.1415927331,  0.0000000794

Go

package main

import "fmt"
import "math"

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

    // 中点則で積分
    var n int = 2
    for j := 1; j <= 10; j++ {
        var h float64 = (b - a) / float64(n)
        var s float64 = 0
        var x float64 = a + (h / 2)
        for i := 1; i <= n; i++ {
            s += f(x)
            x += h
        }
        s *= h
        n *= 2

        // 結果を π と比較
        fmt.Printf("%2d : %13.10f, %13.10f\n", j, s, s - math.Pi)
    }
}

func f(x float64) float64 {
    return (4 / (1 + x * x))
}
Z:\>8g GO0602.go

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

Z:\>GO0602
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Scala

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

        // 中点則で積分
        var n:Int = 2
        for (j <- 1 to 10) {
            var h:Double = (b - a) / n
            var s:Double = 0
            var x:Double = a + (h / 2)
            for (i <- 1 to n) {
                s += f(x)
                x += h
            }
            s *= h
            n *= 2

            // 結果を π と比較
            println("%3d : %13.10f, %13.10f".format(j, s, s - Math.PI))
        }
    }

    def f(x:Double):Double = {
        4 / (1 + x * x)
    }
}
Z:\>scala Scala0602.scala
  1 :  3.1623529412,  0.0207602876
  2 :  3.1468005184,  0.0052078648
  3 :  3.1428947296,  0.0013020760
  4 :  3.1419181743,  0.0003255207
  5 :  3.1416740338,  0.0000813802
  6 :  3.1416129986,  0.0000203451
  7 :  3.1415977399,  0.0000050863
  8 :  3.1415939252,  0.0000012716
  9 :  3.1415929715,  0.0000003179
 10 :  3.1415927331,  0.0000000795

F#

module Fs0602

open System

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

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

// 中点則で積分
let mutable n = 2
for j in [1..10] do
    let         h:double = (b - a) / (double n)
    let mutable s:double = 0.0
    let mutable x:double = a + (h / 2.0)

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

    // 結果を π と比較
    printfn "%3d : %13.10f, %13.10f" j s (s - Math.PI)

exit 0
Z:\>fsi  --nologo --quiet Fs0602.fs
  1 :  3.1623529412,  0.0207602876
  2 :  3.1468005184,  0.0052078648
  3 :  3.1428947296,  0.0013020760
  4 :  3.1419181743,  0.0003255207
  5 :  3.1416740338,  0.0000813802
  6 :  3.1416129986,  0.0000203451
  7 :  3.1415977399,  0.0000050863
  8 :  3.1415939252,  0.0000012716
  9 :  3.1415929715,  0.0000003179
 10 :  3.1415927331,  0.0000000795

Clojure

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

(def a 0)
(def b 1)

; 中点則で積分
(doseq [j (range 1 11)]
    (def n (Math/pow 2 j))
    (def h (/ (- b a) n))
    (def x (+ a (/ h 2)))
    (def s 0)
    (doseq [i (range 1 (+ n 1))]
        (def s (+ s (f x)))
        (def x (+ x h)))
    (def t1 (double (* h s)))
    (def t2 (- t1 (. Math PI)))
    ; 結果を π と比較
    (println (format "%2d : %13.10f, %13.10f" j t1 t2)))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0602.clj
 1 :  3.1623529412,  0.0207602876
 2 :  3.1468005184,  0.0052078648
 3 :  3.1428947296,  0.0013020760
 4 :  3.1419181743,  0.0003255207
 5 :  3.1416740338,  0.0000813802
 6 :  3.1416129986,  0.0000203451
 7 :  3.1415977399,  0.0000050863
 8 :  3.1415939252,  0.0000012716
 9 :  3.1415929715,  0.0000003179
10 :  3.1415927331,  0.0000000795

Haskell

import Text.Printf
import Control.Monad

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

main = do
    forM_ ([1..10::Integer]) $ \j -> do
        let n = 2 ^ j
        let a = 0.0
        let b = 1.0
        let h = (b - a) / (fromIntegral n)
        -- 中点則で積分
        let a2 = a + (h / 2)
        let w1 = sum $ map(\x -> f x) $ map(\i -> (fromIntegral i) * h + a2) $ [0..(n-1)]
        let t1 = h * w1
        -- 結果を π と比較
        let t2 = t1 - pi
        printf "%3d : %13.10f, %13.10f\n" j t1 t2
Z:\>runghc Hs0602.hs
  1 :  3.1623529412,  0.0207602876
  2 :  3.1468005184,  0.0052078648
  3 :  3.1428947296,  0.0013020760
  4 :  3.1419181743,  0.0003255207
  5 :  3.1416740338,  0.0000813802
  6 :  3.1416129986,  0.0000203451
  7 :  3.1415977399,  0.0000050863
  8 :  3.1415939252,  0.0000012716
  9 :  3.1415929715,  0.0000003179
 10 :  3.1415927331,  0.0000000795
inserted by FC2 system