さまざまな言語で数値計算
Only Do What Only You Can Do
シンプソン則
関数 $ f(x) $ の $ a, b $ 間の $n$ 個の微小区間を, 二次曲線で近似し面積を合計する.
3点 $ (x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2)) $ を通る二次曲線の方程式 $ g(x) $ は,
で表される(ラグランジュ補間).
$ x_1 - x_0 $ を $h$ として, 右辺第1項を積分すると
同様に右辺第2項は
同じく第3項は
よって, 関数 $ g(x) $ の $ x_0, x_2 $ 区間の定積分は
と表すことができる.
まとめると, 関数 $ 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 'Simpson則で積分 Dim n: n = 2 Dim i, j For j = 1 To 5 Dim h: h = (b - a) / n Dim s2: s2 = 0 Dim s4: s4 = 0 Dim x: x = a + h For i = 1 To n \ 2 s4 = s4 + f(x) x = x + h s2 = s2 + f(x) x = x + h Next s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 = s4 * 4 Dim s: s = (s2 + s4) * h / 3 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:\0603.vbs 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
JScript
var a = 0 var b = 1 // Simpson則で積分 var n = 2 for (var j = 1; j <= 5; j++) { var h = (b - a) / n var s2 = 0 var s4 = 0 var x = a + h for (var i = 1; i <= n / 2; i++) { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s = (s2 + s4) * h / 3 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:\0603.js 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
PowerShell
function f($x) { 4 / (1 + $x * $x) } $a = 0 $b = 1 # Simpson則で積分 $n = 2; foreach ($j in 1..5) { $h = ($b - $a) / $n $s2 = 0 $s4 = 0 $x = $a + $h foreach ($i in 1..($n / 2)) { $s4 += (f $x) $x += $h $s2 += (f $x) $x += $h } $s2 = ($s2 - (f $b)) * 2 + (f $a) + (f $b) $s4 *= 4 $s = ($s2 + $s4) * $h / 3 $n *= 2 # 結果を π と比較 Write-Host ([String]::Format("{0,2:D} : {1,13:F10}, {2,13:F10}", $j, $s, $s - [Math]::PI)) }
Z:\>powershell -file Z:\0603.ps1 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, 0.0000000000
Perl
use Math::Trig 'pi'; my $a = 0; my $b = 1; # Simpson則で積分 my $n = 2; for $j (1..5) { my $h = ($b - $a) / $n; my $s2 = 0; my $s4 = 0; my $x = $a + $h; for $i (1..($n / 2)) { $s4 += f($x); $x += $h; $s2 += f($x); $x += $h; } $s2 = ($s2 - f($b)) * 2 + f($a) + f($b); $s4 *= 4; my $s = ($s2 + $s4) * $h / 3; $n *= 2; # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", $j, $s, $s - pi); } sub f { my ($x) = @_; 4 / (1 + $x * $x); }
Z:\>perl Z:\0603.pl 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
PHP
<?php $a = 0; $b = 1; # Simpson則で積分 $n = 2; foreach (range(1, 5) as $j) { $h = ($b - $a) / $n; $s2 = 0; $s4 = 0; $x = $a + $h; foreach (range(1, ($n / 2)) as $i) { $s4 += f($x); $x += $h; $s2 += f($x); $x += $h; } $s2 = ($s2 - f($b)) * 2 + f($a) + f($b); $s4 *= 4; $s = ($s2 + $s4) * $h / 3; $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:\0603.php 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
Python
# coding: Shift_JIS import math def f(x): return 4 / (1 + x * x) a = 0 b = 1 # Simpson則で積分 n = 2 for j in range(1, 6): h = (b - a) / float(n) s2 = 0 s4 = 0 x = a + h for i in range(1, (n / 2 + 1)): s4 += f(x) x += h s2 += f(x) x += h s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 s = (s2 + s4) * h / 3 n *= 2 # 結果を π と比較 print "%2d : %13.10f, %13.10f" % (j, s, s - math.pi)
Z:\>python Z:\0603.py 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
Ruby
def f(x) 4 / (1 + x * x) end a = 0 b = 1 # Simpson則で積分 n = 2 (1..5).each do |j| h = (b - a) / n.to_f s2 = 0 s4 = 0 x = a + h (1..(n / 2)).each do |i| s4 += f(x) x += h s2 += f(x) x += h end s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 s = (s2 + s4) * h / 3 n *= 2 # 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, s, s - Math::PI) end
Z:\>ruby Z:\0603.rb 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
Groovy
Pascal
program Pas0603(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; s2, s4 :Double; begin // Simpson則で積分 n := 2; for j := 1 to 5 do begin h := (b - a) / n; s2 := 0; s4 := 0; x := a + h; for i := 1 to (n div 2) do begin s4 := s4 + f(x); x := x + h; s2 := s2 + f(x); x := x + h; end; s2 := (s2 - f(b)) * 2 + f(a) + f(b); s4 := s4 * 4; s := (s2 + s4) * h / 3; n := n * 2; // 結果を π と比較 writeln(format('%2d : %13.10f, %13.10f', [j, s, s - PI])); end; end.
Z:\>fpc -v0 -l- Pas0603.pp Z:\>Pas0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 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 Ada0603 is a : Constant Long_Float := 0.0; b : Constant Long_Float := 1.0; h, s, x : Long_Float; s2, s4 : Long_Float; n : Integer; function f(x:Long_Float) return Long_Float is begin return 4.0 / (1.0 + x * x); end f; begin -- Simpson則で積分 n := 2; for j in 1..5 loop h := (b - a) / Long_Float(n); s2 := 0.0; s4 := 0.0; x := a + h; for i in 1..(n / 2) loop s4 := s4 + f(x); x := x + h; s2 := s2 + f(x); x := x + h; end loop; s2 := (s2 - f(b)) * 2.0 + f(a) + f(b); s4 := s4 * 4.0; s := (s2 + s4) * h / 3.0; 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 Ada0603;
xxxxxx@yyyyyy /Z $ gnatmake Ada0603.adb xxxxxx@yyyyyy /Z $ Ada0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
VB.NET
Module VB0603 Public Sub Main() Const a As Double = 0 Const b As Double = 1 'Simpson則で積分 Dim n As Integer = 2 For j As Integer = 1 To 5 Dim h As Double = (b - a) / n Dim s2 As Double = 0 Dim s4 As Double = 0 Dim x As Double = a + h For i As Integer = 1 To n \ 2 s4 += f(x) x += h s2 += f(x) x += h Next s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 Dim s As Double = (s2 + s4) * h / 3 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 VB0603.vb Z:\>VB0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, 0.0000000000
C#
using System; public class CS0603 { private static double f(double x) { return 4 / (1 + x * x); } public static void Main() { const double a = 0; const double b = 1; // Simpson則で積分 int n = 2; for (int j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (int i = 1; i <= n / 2; i++) { s4 += f(x); x += h; s2 += f(x); x += h; } s2 = (s2 - f(b)) * 2 + f(a) + f(b); s4 *= 4; double s = (s2 + s4) * h / 3; n *= 2; // 結果を π と比較 Console.WriteLine(string.Format("{0,2:D} : {1,13:F10}, {2,13:F10}", j, s, s - Math.PI)); } } }
Z:\>csc -nologo CS0603.cs Z:\>CS0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, 0.0000000000
Java
public class Java0603 { 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; // Simpson則で積分 int n = 2; for (int j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (int i = 1; i <= n / 2; i++) { s4 += f(x); x += h; s2 += f(x); x += h; } s2 = (s2 - f(b)) * 2 + f(a) + f(b); s4 *= 4; double s = (s2 + s4) * h / 3; n *= 2; // 結果を π と比較 System.out.println(String.format("%2d : %13.10f, %13.10f", j, s, s - Math.PI)); } } }
Z:\>javac Java0603.java Z:\>java Java0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 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; // Simpson則で積分 int n = 2; for (int j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (int i = 1; i <= n / 2; i++) { s4 += f(x); x += h; s2 += f(x); x += h; } s2 = (s2 - f(b)) * 2 + f(a) + f(b); s4 *= 4; double s = (s2 + s4) * h / 3; 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 CP0603.cpp cp0603.cpp: Z:\>CP0603 1: 3.1333333333, -0.0082593203 2: 3.1415686275, -0.0000240261 3: 3.1415925025, -0.0000001511 4: 3.1415926512, -0.0000000024 5: 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; // Simpson則で積分 int n = 2; int i, j; for (j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (i = 1; i <= n / 2; i++) { s4 += f(x); x += h; s2 += f(x); x += h; } s2 = (s2 - f(b)) * 2 + f(a) + f(b); s4 *= 4; double s = (s2 + s4) * h / 3; n *= 2; // 結果を π と比較 printf("%2d : %13.10f, %13.10f\n", j, s, s - M_PI); } return 0; }
xxxxxx@yyyyyy /Z $ gcc -o OC0603 OC0603.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
D
import std.stdio; import std.math; void main(string[] args) { const double a = 0; const double b = 1; // Simpson則で積分 int n = 2; for (int j = 1; j <= 5; j++) { double h = (b - a) / n; double s2 = 0; double s4 = 0; double x = a + h; for (int i = 1; i <= n / 2; i++) { s4 += f(x); x += h; s2 += f(x); x += h; } s2 = (s2 - f(b)) * 2 + f(a) + f(b); s4 *= 4; double s = (s2 + s4) * h / 3; n *= 2; // 結果を π と比較 writefln("%2d : %13.10f, %13.10f", j, s, s - PI); } } double f(double x) { return 4 / (1 + x * x); }
Z:\>dmd D0603.d Z:\>D0603 1 : 3.1333333333, -0.0082593202 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000023 5 : 3.1415926536, -0.0000000000
Go
package main import "fmt" import "math" func main() { const a = 0 const b = 1 // Simpson則で積分 var n int = 2 for j := 1; j <= 5; j++ { var h float64 = (b - a) / float64(n) var s2 float64 = 0 var s4 float64 = 0 var x float64 = a + h for i := 1; i <= n / 2; i++ { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s float64 = (s2 + s4) * h / 3 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 GO0603.go Z:\>8l -o GO0603.exe GO0603.8 Z:\>GO0603 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
Scala
object Scala0603 { def main(args: Array[String]) { val a:Double = 0 val b:Double = 1 // Simpson則で積分 var n:Int = 2 for (j <- 1 to 5) { var h:Double = (b - a) / n var s2:Double = 0 var s4:Double = 0 var x:Double = a + h for (i <- 1 to (n / 2)) { s4 += f(x) x += h s2 += f(x) x += h } s2 = (s2 - f(b)) * 2 + f(a) + f(b) s4 *= 4 var s:Double = (s2 + s4) * h / 3 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 Scala0603.scala 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
F#
module Fs0603 open System let f(x:double):double = 4.0 / (1.0 + x * x) let a:double = 0.0 let b:double = 1.0 // Simpson則で積分 let mutable n = 2 for j in [1..5] do let h:double = (b - a) / (double n) let mutable s2:double = 0.0 let mutable s4:double = 0.0 let mutable x:double = a + h for i in [1..(n / 2)] do s4 <- s4 + (f x) x <- x + h s2 <- s2 + (f x) x <- x + h s2 <- (s2 - (f b)) * 2.0 + (f a) + (f b) s4 <- s4 * 4.0 let s = (s2 + s4) * h / 3.0 n <- n * 2 // 結果を π と比較 printfn "%3d : %13.10f, %13.10f" j s (s - System.Math.PI) exit 0
Z:\>fsi --nologo --quiet Fs0603.fs 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, 0.0000000000
Clojure
(defn f[x] (/ 4 (+ 1 (* x x)))) (def a 0) (def b 1) ; Simpson則で積分 (doseq [j (range 1 6)] (def n (Math/pow 2 j)) (def h (/ (- b a) n)) (def x (+ a h)) (def s2 0) (def s4 0) (doseq [i (range 1 (+ (/ n 2) 1))] (def s4 (+ s4 (f x))) (def x (+ x h)) (def s2 (+ s2 (f x))) (def x (+ x h)) ) (def w2 (double (+ (* (- s2 (f b)) 2) (f a) (f b)))) (def w4 (* s4 4)) (def t1 (double (/ (* (+ w2 w4) h) 3))) (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 Clj0603.clj 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000
Haskell
import Text.Printf import Control.Monad f::Double->Double f x = 4 / (1 + (x * x)) main = do forM_ ([1..5::Integer]) $ \j -> do let n = 2 ^ j let a = 0.0 let b = 1.0 let h = (b - a) / (fromIntegral n) -- シンプソン則で積分 let w4 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2 - 1) $ [1..(div n 2)] let w2 = sum $ map(\x -> (f ((fromIntegral x) * h))) $ map(\i -> i * 2) $ [1..(div n 2)] let t1 = (w2 - (f b)) * 2 + (f a) + (f b) let t2 = w4 * 4 let t3 = (t1 + t2) * h / 3 -- 結果を π と比較 let t4 = t3 - pi printf "%3d : %13.10f, %13.10f\n" j t3 t4
Z:\>runghc Hs0603.hs 1 : 3.1333333333, -0.0082593203 2 : 3.1415686275, -0.0000240261 3 : 3.1415925025, -0.0000001511 4 : 3.1415926512, -0.0000000024 5 : 3.1415926536, -0.0000000000