Zenn
📌

PHPを使って√2を求める

2020/10/20に公開
PHP
tech

概要

\sqrt{2}を求める問題は数値計算ではFuzzBuzz並の基礎問題である、方法はf(x)=x^2-2=0を数値的に解くというものである。
この解法として有名なのは二分法とニュートン(ラフソン)法がある。
私はPHPはほぼ書いたことがないがこれを練習問題にしてみようと思う。

設計のようなもの

方程式を解くということは解法によらないはずである。そのため解き方というアルゴリズムの実装をインターフェイスとして分離する。
また、例はf(x)=x^2-2=0であるがPHPはラムダ式が使えるのでf(x)自体を引数にすることで汎用性を高める。

解法について

二分法もニュートン法も繰り返し法の一種である誤差に対してそれよりも予測される誤差が小さくなったときにそれを解とする方法である。(なので真の\sqrt{2}ではなく一定の誤差がある)
紹介した2つの方法は簡単には以下のようであり、二分法は解がある限り安定的にその解を探すことができる。一方でニュートン法は収束が速い一方、解から遠く関数系が複雑なときには解には向かわずに振動する場合がある。

二分法

いわゆる挟み撃ちである、解が間にある2つの点を取るとf(x_0)<0かつf(x_1)>0となる原理を用いて中点をx_{center}=0.5(x_0+x_1)としてそれがf(x_{center})の符号で解の範囲を狭めていく方法である。

ニュートン法

微分を用いた方法である初期値から一次式で近似して次の解を求めて、その差(+関数の値)から収束しているかを判断して解を求める方法である。

実装

Solver

方程式を解くモノのインターフェイスとしてSolverを用意する、引数は第一引数を解かせたい方程式を関数オブジェクトとして第二引数に初期値にする。
また、打ち切り誤差は定数としてSolverに持たせておく。
その実装として二分法(BySelection)とニュートン法(Newton)クラスを用意する。
内部実装は以下のとおりである。

Solver.php
<?php
interface Solver{
    const threshold=1.0e-8;
    public function solve($func, $x);
}

class BySelection implements Solver{
    public function solve($func, $x){
        $x1=0.5*$x;
        $x2=1.5*$x;
        $y =$func($x);
        $y1=$func($x1);
        $y2=$func($x2);

        $count=0;
        if( $y1*$y2>0 ) throw new Exception('二分法、初期値エラー');
        while( true ){
            $count++;
            if( $y*$y1>0 ){
                $x1=$x;
                $y1=$y;
            }
            else{
                $x2=$x;
                $y2=$y;
            }
            $x=0.5*($x1+$x2);
            $y=$func($x);

            if( abs($y1)<Solver::threshold && abs($y2)<Solver::threshold && abs($x1-$x2)<Solver::threshold ) return $x;
            if( $count==1000 ) throw new Exception('二分法で収束しませんでした');
        }
    }
}

class Newton implements Solver{
    public function solve($func, $x){
        $count=0;
        while( true ){
            $count++;
            $y=$func($x);
            $dy=($func($x+Solver::threshold)-$func($x-Solver::threshold))/(2.0*Solver::threshold);

            if( abs($y)<Solver::threshold && abs($y/$dy)<Solver::threshold ) return $x;

            $x=$x-$y/$dy;
            if( $count==1000 ) throw new Exception('ニュートン法で収束しませんでした');
        }
    }
}
?>

二分法はアルゴリズムとしては単純なので読めると思う。
ニュートン法では直線近似のために微分を使うがそれは数値微分を用いた。

main関数

実際に使う側では以下のようになる。
$solverの実体をBySelectionかNewtonかで内部のアルゴリズムは変更できるが結果は変わらない（はずである）。

main.php
<?php
require_once('Solver.php');

// $func=function($x){ return cos($x); };
$func=function($x){ return x*x-2; };

// $sovler=new BySelection(); 
$solver=new Newton();

var_dump(2.0*$solver->solve($func, 2.0));
?>

不備など

一応、これで想定したケースはいけるが一般の関数系となると初期値などにより解が求まらないことがある。（二分法の初期値の雑な処理等が最たる例）
数値計算では収束条件(threshold)やイテレーション回数等はアルゴリズムや解きたい問題に大きく依存するのでここまでのブラックボック化が適切ではないとおもう。

あとがきと感想

PHPをほとんど使ったことがない筆者が単体で動かせるサンプルとして数値計算をさせてみた、PHPがその用途に向いてるとは思わないがインターフェイスやラムダ式を使うことでかなりの対応ができると思った。
また、PHPの定数の使い方や例外処理の飛ばし方がわかった（このレベル）。

参考文献

Discussion