さまざまな言語で数値計算
Only Do What Only You Can Do
ニュートン法
非線形方程式の解法(ニュートン法)を利用して2の平方根を求める .
1. まず, 関数 $ f(x) $ 上の点 $ (x_0,f(x_0)) $ を考える.
2. 点 $ (x_0,f(x_0)) $ における $ f(x) $ の接線と $x$ 軸との交点 $x_1$ は $x_0$ より解に近づいている.
3. この作業を繰り返して行くことで解を求める.
点 $ (x_0,f(x_0)) $ での接線の傾きは
なので,
この式を漸化式として用いる.
VBScript
Option Explicit Dim x: x = 2.0 WScript.StdOut.Write Right(Space(12) & FormatNumber(newton(x), 10, -1, 0, 0), 12) & vbNewLine Private Function newton(ByVal x0) Dim x1 Do While(True) x1 = x0 - (f0(x0) / 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 newton = x1 End Function Private Function f0(ByVal x) f0 = x * x - 2 End Function Private Function f1(ByVal x) f1 = 2 * x End Function
Z:\>cscript //nologo Z:\0904.vbs 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
JScript
var x = 2 WScript.StdOut.Write((" " + newton(x).toFixed(10) ).slice(-12) + "\n") function newton(x0) { var x1 while (true) { x1 = x0 - (f0(x0) / 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 }
Z:\>cscript //nologo Z:\0904.js 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
PowerShell
function newton($x0) { while ($true) { $x1 = $x0 - ((f0 $x0) / (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 } $x = 2 Write-Host ([String]::Format("{0,12:F10}", (newton $x)))
Z:\>powershell -file Z:\0904.ps1 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
Perl
my $x = 2; printf("%12.10f\n", newton($x)); sub newton { my ($x0) = @_; my $x1; while (1) { $x1 = $x0 - (f0($x0) / 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; }
Z:\>perl Z:\0904.pl 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
PHP
<?php $x = 2; printf("%12.10f\n", newton($x)); function newton($x0) { while (TRUE) { $x1 = $x0 - (f0($x0) / 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; } ?>
Z:\>php Z:\0904.php 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton(x0): while True: x1 = x0 - (f0(x0) / 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" % newton(x)
Z:\>python Z:\0904.py 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton(x0) while true x1 = x0 - (f0(x0) / 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", newton(x))
Z:\>ruby Z:\0904.rb 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
Groovy
x = 2 printf ("%12.10f\n" , newton(x)) def newton(x0) { while (true) { x1 = x0 - (f0(x0) / 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 }
Z:\>groovy Groovy0904.groovy 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
Pascal
program Pas0904(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 newton(x0:Double):Double; var x1: Double; begin while true do begin x1 := x0 - (f0(x0) / 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', [newton(x)])); end.
Z:\>fpc -v0 -l- Pas0904.pp Z:\>Pas0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 Ada0904 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 newton(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)); 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 newton; x: Long_Float := 2.0; begin Put(newton(x), Fore=>2, Aft=>10, Exp=>0); New_Line; end Ada0904;
xxxxxx@yyyyyy /Z $ gnatmake Ada0904.adb xxxxxx@yyyyyy /Z $ Ada0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
VB.NET
Option Explicit Module VB0904 Public Sub Main() Dim x As Double = 2.0 Console.WriteLine(String.Format("{0,12:F10}", newton(x))) End Sub Private Function newton(ByVal x0 As Double) Dim x1 As Double Do While(True) x1 = x0 - (f0(x0) / 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 End Module
Z:\>vbc -nologo VB0904.vb Z:\>VB0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
C#
using System; public class CS0904 { public static void Main() { double x = 2; Console.WriteLine(string.Format("{0,12:F10}", newton(x))); } private static double newton(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / 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; } }
Z:\>csc -nologo CS0904.cs Z:\>CS0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
Java
import static java.lang.System.out; public class Java0904 { public static void main(String []args) { double x = 2; out.println(String.format("%12.10f", newton(x))); } private static double newton(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / 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; } }
Z:\>javac Java0904.java Z:\>java Java0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / 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) << newton(x) << endl; return 0; }
Z:\>bcc32 -q CP0904.cpp cp0904.cpp: Z:\>CP0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton(double x0) { double x1; while (YES) { x1 = x0 - (f0(x0) / 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", newton(x)); return 0; }
xxxxxx@yyyyyy /Z $ gcc -o OC0904 OC0904.m -lobjc -lgnustep-base -I $INCLUDE -L $LIB $CFLAGS xxxxxx@yyyyyy /Z $ OC0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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", newton(x)); } double newton(double x0) { double x1; while (true) { x1 = x0 - (f0(x0) / 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; }
Z:\>dmd D0904.d Z:\>D0904 1.5000000000 0.0857864376 1.4166666667 0.0024531042 1.4142156863 0.0000021239 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", newton(x)) } func newton(x0 float64) float64 { var x1 float64 for { x1 = x0 - (f0(x0) / 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 }
Z:\>8g GO0904.go Z:\>8l -o GO0904.exe GO0904.8 Z:\>GO0904 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
Scala
object Scala0904 { def main(args: Array[String]) { val x = 2.0 println("%12.10f".format(newton(x))) } def newton(x0:Double):Double = { val x1 = x0 - (f0(x0) / f1(x0)) println("%12.10f\t%13.10f".format(x1, x1 - Math.sqrt(2))) if (Math.abs(x1 - x0) < 0.0000000001) x1 else newton(x1) } def f0(x:Double) = { x * x - 2 } def f1(x:Double) = { 2 * x } }
Z:\>scala Scala0904.scala 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624 0.0000000000 1.4142135624
F#
module Fs0904 open System let f0 (x:double):double = x * x - 2.0 let f1 (x:double):double = 2.0 * x let rec newton (x0:double):double = let x1 = x0 - (f0(x0) / f1(x0)) printfn "%12.10f\t%13.10f" x1 (x1 - Math.Sqrt(2.0)) if abs(x1 - x0) < 0.0000000001 then x1 else newton x1 let x = 2.0 printfn "%12.10f" (newton x) exit 0
Z:\>fsi --nologo --quiet Fs0904.fs 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton [x0] (def x1 (- x0 (/ (f0 x0) (f1 x0)))) (println (format "%12.10f\t%13.10f" x1 (- x1 (Math/sqrt 2)))) (if (< (. Math abs (- x1 x0)) 0.00000000001) x1 (newton x1))) (def x 2.0) (println (format "%12.10f" (newton x)))
Z:\>java -cp C:\ProgramFiles\clojure-1.5.1\clojure-1.5.1.jar clojure.main Clj0904.clj 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 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 newton::Double->Double newton x0 = let x1 = x0 - ((f0 x0) / (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))) newton x1 main = do let x = 2.0 printf "%12.10f\n" (newton x)
Z:\>runghc Hs0904.hs 1.5000000000 0.0857864376 1.4166666667 0.0024531043 1.4142156863 0.0000021239 1.4142135624 0.0000000000 1.4142135624