ベイリー法

home > さまざまな言語で数値計算 > 非線形方程式 >

ベイリー法

非線形方程式の解法（ベイリー法）を利用して2の平方根を求める .

関数 $ f(x) $ をテイラー展開すると,

と, 近似できる.
この式の二次の項までを左辺 $ = 0, x = x_{n+1} $ として変形すると,

この式の右辺の $ x_{n+1} $ をニュートン法で使った

で置き換えると

この式を漸化式として用いる.

VBScript

Option Explicit

Dim x: x = 2.0
WScript.StdOut.Write Right(Space(12) & FormatNumber(bailey(x), 10, -1, 0, 0), 12) & vbNewLine

Private Function bailey(ByVal x0)
    Dim x1
    Do While(True)
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        WScript.StdOut.Write Right(Space(12) & FormatNumber(x1,          10, -1, 0, 0), 12) & vbTab
        WScript.StdOut.Write Right(Space(12) & FormatNumber(x1 - Sqr(2), 10, -1, 0, 0), 12) & vbNewLine

        If Abs(x1 - x0) < 0.0000000001 Then Exit Do
        x0 = x1
    Loop

    bailey = x1
End Function

Private Function f0(ByVal x)
    f0 = x * x - 2
End Function

Private Function f1(ByVal x)
    f1 = 2 * x
End Function

Private Function f2(ByVal x)
    f2 = 2
End Function

Z:\>cscript //nologo Z:\0905.vbs
1.4285714286    0.0143578662
1.4142139268    0.0000003644
1.4142135624    0.0000000000
1.4142135624    0.0000000000
1.4142135624

JScript

var x = 2
WScript.StdOut.Write(("            " + bailey(x).toFixed(10) ).slice(-12) + "\n")

function bailey(x0) {
    var x1
    while (true) {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        WScript.StdOut.Write(("            " + x1.toFixed(10)                  ).slice(-12) + "\t")
        WScript.StdOut.Write(("            " + (x1 - Math.sqrt(2)).toFixed(10) ).slice(-12) + "\n")

        if (Math.abs(x1 - x0) < 0.0000000001) break
        x0 = x1
    }
    return x1
}

function f0(x) {
    return x * x - 2
}

function f1(x) {
    return 2 * x
}

function f2(x) {
    return 2
}

Z:\>cscript //nologo Z:\0905.js
1.4285714286    0.0143578662
1.4142139268    0.0000003644
1.4142135624    0.0000000000
1.4142135624    0.0000000000
1.4142135624

PowerShell

function bailey($x0)
{
    while ($true)
    {
        $x1 = $x0 - ((f0 $x0) / ((f1 $x0) - ((f0 $x0) * (f2 $x0) / (2 * (f1 $x0)))))
        Write-Host ([String]::Format("{0,12:F10}`t{1,13:F10}", $x1, $x1 - [Math]::Sqrt(2)))

        if ([Math]::Abs($x1 - $x0) -lt 0.0000000001)
        {
            break
        }
        $x0 = $x1
    }
    return $x1
}

function f0($x)
{
    return $x * $x - 2
}

function f1($x)
{
    return 2 * $x
}

function f2($x)
{
    return 2
}

$x = 2
Write-Host ([String]::Format("{0,12:F10}", (bailey $x)))

Z:\>powershell -file Z:\0905.ps1
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Perl

my $x = 2;
printf("%12.10f\n", bailey($x));

sub bailey
{
    my ($x0) = @_;
    my $x1;
    while (1)
    {
        $x1 = $x0 - (f0($x0) / (f1($x0) - (f0($x0) * f2($x0) / (2 * f1($x0)))));
        printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));

        if (abs($x1 - $x0) < 0.0000000001)
        {
            last;
        }
        $x0 = $x1;
    }
    $x1;
}

sub f0
{
    my ($x) = @_;
    $x * $x - 2;
}

sub f1
{
    my ($x) = @_;
    2 * $x;
}

sub f2
{
    my ($x) = @_;
    2;
}

Z:\>perl Z:\0905.pl
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

PHP

<?php
$x = 2;
printf("%12.10f\n", bailey($x));

function bailey($x0)
{
    while (TRUE)
    {
        $x1 = $x0 - (f0($x0) / (f1($x0) - (f0($x0) * f2($x0) / (2 * f1($x0)))));
        printf("%12.10f\t%13.10f\n", $x1, $x1 - sqrt(2));

        if (abs($x1 - $x0) < 0.0000000001) break;

        $x0 = $x1;
    }
    return $x1;
}

function f0($x)
{
    return $x * $x - 2;
}

function f1($x)
{
    return 2 * $x;
}

function f2($x)
{
    return 2;
}
?>

Z:\>php Z:\0905.php
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

Python

# coding: Shift_JIS

import math

def f0(x):
    return x * x - 2.0

def f1(x):
    return 2.0 * x

def f2(x):
    return 2.0

def bailey(x0):
    while True:
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        print "%12.10f\t%13.10f" % (x1, x1 - math.sqrt(2))

        if (abs(x1 - x0) < 0.0000000001):
            break
        x0 = x1

    return x1

x = 2.0
print "%12.10f" % bailey(x)

Z:\>python Z:\0905.py
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

Ruby

include Math

def f0(x)
    x * x - 2.0
end

def f1(x)
    2.0 * x
end

def f2(x)
    2.0
end

def bailey(x0)
    while true
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        printf("%12.10f\t%13.10f\n", x1, x1 - sqrt(2))

        if (x1 - x0).abs < 0.0000000001
            break
        end

        x0 = x1
    end

    x1
end

x = 2.0
printf("%12.10f\n", bailey(x))

Z:\>ruby Z:\0905.rb
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

Groovy

x = 2
printf ("%12.10f\n" , bailey(x))

def bailey(x0) {
    while (true) {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        printf("%12.10f\t%13.10f\n", x1, x1 - Math.sqrt(2))

        if (Math.abs(x1 - x0) < 0.0000000001) break
        x0 = x1
    }
    x1
}

def f0(x) {
    x * x - 2
}

def f1(x) {
    2 * x
}

def f2(x) {
    2
}

Z:\>groovy Groovy0905.groovy
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Pascal

program Pas0905(arg);
{$MODE delphi}

uses
    SysUtils, Math;

function f0(x:Double):Double;
begin
    result := x * x - 2;
end;

function f1(x:Double):Double;
begin
    result := 2 * x;
end;

function f2(x:Double):Double;
begin
    result := 2;
end;

function bailey(x0:Double):Double;
var
    x1: Double;
begin
    while true do
    begin
        x1 := x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
        writeln(format('%12.10f'#9'%13.10f', [x1, x1 - Sqrt(2)]));

        if Abs(x1 - x0) < 0.0000000001 then break;
        x0 := x1;
    end;

    result := x1;
end;

var
    x:  Double = 2.0;
begin
    writeln(format('%12.10f', [bailey(x)]));
end.

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

Z:\>Pas0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Ada

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

procedure Ada0905 is

    function f0(x:Long_Float) return Long_Float is
    begin
        return x * x - 2.0;
    end f0;

    function f1(x:Long_Float) return Long_Float is
    begin
        return 2.0 * x;
    end f1;

    function f2(x:Long_Float) return Long_Float is
    begin
        return 2.0;
    end f2;

    function bailey(x:Long_Float) return Long_Float is
        x0: Long_Float;
        x1: Long_Float;
    begin
        x0 := x;
        while true loop
            x1 := x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2.0 * f1(x0)))));
            Put(x1,             Fore=>2, Aft=>10, Exp=>0);
            Put(Ascii.HT);
            Put(x1 - Sqrt(2.0), Fore=>2, Aft=>10, Exp=>0);
            New_Line;

            if Abs(x1 - x0) < 0.0000000001 then
                exit;
            end if;
            x0 := x1;
        end loop;

        return x1;
    end bailey;

    x: Long_Float := 2.0;
begin
    Put(bailey(x), Fore=>2, Aft=>10, Exp=>0);
    New_Line;
end Ada0905;

xxxxxx@yyyyyy /Z
$ gnatmake Ada0905.adb

xxxxxx@yyyyyy /Z
$ Ada0905
 1.4285714286    0.0143578662
 1.4142139268    0.0000003644
 1.4142135624    0.0000000000
 1.4142135624    0.0000000000
 1.4142135624

VB.NET

Option Explicit

Module VB0905
    Public Sub Main()
        Dim x As Double = 2.0
        Console.WriteLine(String.Format("{0,12:F10}", bailey(x)))
    End Sub

    Private Function bailey(ByVal x0 As Double)
        Dim x1 As Double
        Do While(True)
            x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
            Console.WriteLine(String.Format("{0,12:F10}{2}{1,13:F10}", x1, x1 - Math.Sqrt(2), vbTab))

            If Math.Abs(x1 - x0) < 0.0000000001 Then Exit Do
            x0 = x1
        Loop

        Return x1
    End Function

    Private Function f0(ByVal x As Double) As Double
        Return x * x - 2
    End Function

    Private Function f1(ByVal x As Double) As Double
        Return 2 * x
    End Function

    Private Function f2(ByVal x As Double) As Double
        Return 2
    End Function
End Module

Z:\>vbc -nologo VB0905.vb

Z:\>VB0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

C#

using System;

public class CS0905
{
    public static void Main()
    {
        double x = 2;
        Console.WriteLine(string.Format("{0,12:F10}", bailey(x)));
    }

    private static double bailey(double x0)
    {
        double x1;
        while (true)
        {
            x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
            Console.WriteLine(string.Format("{0,12:F10}\t{1,13:F10}", x1, x1 - Math.Sqrt(2)));

            if (Math.Abs(x1 - x0) < 0.0000000001) break;
            x0 = x1;
        }
        return x1;
    }

    private static double f0(double x)
    {
        return x * x - 2;
    }
    private static double f1(double x)
    {
        return 2 * x;
    }
    private static double f2(double x)
    {
        return 2;
    }
}

Z:\>csc -nologo CS0905.cs

Z:\>CS0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Java

import static java.lang.System.out;

public class Java0905 {

    public static void main(String []args) {
        double x = 2;
        out.println(String.format("%12.10f", bailey(x)));
    }

    private static double bailey(double x0) {
        double x1;
        while (true) {
            x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
            out.println(String.format("%12.10f\t%13.10f", x1, x1 - Math.sqrt(2)));

            if (Math.abs(x1 - x0) < 0.0000000001) break;
            x0 = x1;
        }
        return x1;
    }

    private static double f0(double x) {
        return x * x - 2;
    }
    private static double f1(double x) {
        return 2 * x;
    }
    private static double f2(double x) {
        return 2;
    }
}

Z:\>javac Java0905.java

Z:\>java Java0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

C++

#include <iostream>
#include <iomanip>
#include <math.h>

using namespace std;

double f0(double x)
{
    return x * x - 2;
}
double f1(double x)
{
    return 2 * x;
}
double f2(double x)
{
    return 2;
}

double bailey(double x0)
{
    double x1;
    while (true)
    {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
        cout << setw(12) << fixed << setprecision(10) << x1 << "\t";
        cout << setw(13) << fixed << setprecision(10) << x1 - sqrt(2) << endl;

        if (fabs(x1 - x0) < 0.0000000001) break;
        x0 = x1;
    }
    return x1;
}

int main()
{
    double x = 2.0;
    cout << setw(12) << fixed << setprecision(10) << bailey(x) << endl;
    return 0;
}

Z:\>bcc32 -q CP0905.cpp
cp0905.cpp:

Z:\>CP0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Objective-C

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

double f0(double x)
{
    return x * x - 2;
}
double f1(double x)
{
    return 2 * x;
}
double f2(double x)
{
    return 2;
}

double bailey(double x0)
{
    double x1;
    while (YES)
    {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
        printf("%12.10f\t%13.10f\n", x1, x1 - sqrt(2));

        if (fabs(x1 - x0) < 0.0000000001) break;
        x0 = x1;
    }
    return x1;
}

int main()
{
    double x = 2;
    printf("%12.10f\n", bailey(x));
    return 0;
}

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

xxxxxx@yyyyyy /Z
$ OC0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

D

import std.stdio;
import std.math;

void main(string[] args)
{
    double x = 2;
    writefln("%12.10f", bailey(x));
}

double bailey(double x0)
{
    double x1;
    while (true)
    {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))));
        writefln("%12.10f\t%13.10f", x1, x1 - sqrt(2.0));

        if (fabs(x1 - x0) < 0.0000000001) break;
        x0 = x1;
    }
    return x1;
}

double f0(double x)
{
    return x * x - 2;
}
double f1(double x)
{
    return 2 * x;
}
double f2(double x)
{
    return 2;
}

Z:\>dmd D0905.d

Z:\>D0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Go

package main

import "fmt"
import "math"

func main() {
    var x float64 = 2
    fmt.Printf("%12.10f\n", bailey(x))
}

func bailey(x0 float64) float64 {
    var x1 float64
    for {
        x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        fmt.Printf("%12.10f\t%13.10f\n", x1, x1 - math.Sqrt(2))

        if math.Fabs(x1 - x0) < 0.0000000001 {
            break
        }
        x0 = x1
    }
    return x1
}

func f0(x float64) float64 {
    return x * x - 2
}
func f1(x float64) float64 {
    return 2 * x
}
func f2(x float64) float64 {
    return 2
}

Z:\>8g GO0905.go

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

Z:\>GO0905
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

Scala

object Scala0905 {

    def main(args: Array[String]) {
        val x = 2.0
        println("%12.10f".format(bailey(x)))
    }

    def bailey(x0:Double):Double = {
        val x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2 * f1(x0)))))
        println("%12.10f\t%13.10f".format(x1, x1 - Math.sqrt(2)))

        if (Math.abs(x1 - x0) < 0.0000000001)
            x1
        else
            bailey(x1)
    }

    def f0(x:Double) = {
        x * x - 2
    }
    def f1(x:Double) = {
        2 * x
    }
    def f2(x:Double) = {
        2
    }
}

Z:\>scala Scala0905.scala
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

F#

module Fs0905

open System

let f0 (x:double):double =
    x * x - 2.0

let f1 (x:double):double =
    2.0 * x

let f2 (x:double):double =
    2.0

let rec bailey (x0:double):double =
    let x1 = x0 - (f0(x0) / (f1(x0) - (f0(x0) * f2(x0) / (2.0 * f1(x0)))))
    printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0))

    if abs(x1 - x0) < 0.0000000001 then
        x1
    else
        bailey x1

let x = 2.0
printfn "%12.10f" (bailey x)

exit 0

Z:\>fsi  --nologo --quiet Fs0905.fs
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624     0.0000000000
1.4142135624     0.0000000000
1.4142135624

Clojure

(defn f0[x]
    (- (* x x) 2.0))

(defn f1[x]
    (* 2.0 x))

(defn f2[x]
    2.0)

(defn bailey [x0]
    (def x1 (- x0
               (/ (f0 x0)
                  (- (f1 x0)
                     (/ (* (f0 x0) (f2 x0))
                        (* 2.0 (f1 x0)))))))
    (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2))))

    (if (< (. Math abs (- x1 x0)) 0.00000000001)
        x1
        (bailey x1)))

(def x 2.0)
(println (format "%12.10f" (bailey x)))

Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0905.clj
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624     0.0000000000
1.4142135624

Haskell

import Text.Printf
import Debug.Trace

f0::Double->Double
f0 x = x * x - 2.0

f1::Double->Double
f1 x = 2.0 * x

f2::Double->Double
f2 x = 2.0

bailey::Double->Double
bailey x0 =
    let
        x1 = x0 - ((f0 x0) / ((f1 x0) - ((f0 x0) * (f2 x0) / (2.0 * (f1 x0)))))
    in
        if abs(x1 - x0) < 0.0000000001 then
            x1
        else
            trace (printf "%12.10f\t%13.10f" x1 (x1 - ((sqrt 2.0)::Double)))
            bailey x1

main = do
    let x = 2.0
    printf "%12.10f\n" (bailey x)

Z:\>runghc Hs0905.hs
1.4285714286     0.0143578662
1.4142139268     0.0000003644
1.4142135624    -0.0000000000
1.4142135624

VBScript
JScript
PowerShell
Perl
PHP
Python
Ruby
Groovy

Pascal
Ada
VB.NET

C#
Java
C++
Objective-C
D
Go

Scala
F#
Clojure
Haskell