TSperl 発言

極悪です。

[TSabc:175] Re: 自作数学関数 の続きで、eval を使って任意の関
数の根を二分法または内分法で求める Perl スクリプトです。

f(x) = 0 となるような x を求めるには、f(xm) * f(yp) < 0 とな
るような xm と xp を探し出して x1 = (xm + xp)/2 としてやれば、
f(x1) は f(xm) や f(xp) よりもよりゼロに近い（x1 は xm や xp
よりも根に近い）はず、という話。

作ってる最中に、「そういえばすでに CPAN に Math::Brent とい
うそれらしいモジュールがあったなぁ・・・」と思い出したのですが、
これは関数の極値（最小値）を求めるモジュールでした。
実はこれを使ってもＯＫなはずで、f(x) が 0 になる点を求めるに
は、f(x) の自乗が最小値をとるような x を求めればいいはずです。
まぁ、いいや。

さて、２の平方根 1.41421356 を求めるには、このスクリプトに
x*x-2 = 0 を解かせます。

%perl root.pl "$x*$x-2"
x = +1.414214   (+1.414214 +1.600000)
y = -0.000000   (-0.000000 +0.560000)

%

eval を使ってるので、任意の式を指定できます。

%perl root.pl "sin($x)/cos($x) - $x - 0.1"
x = +0.631659   (+0.631659 +0.800000)
y = -0.000000   (-0.000000 +0.129639)

%

# perl の標準関数には tan がないのか、ってエラーになってから
# 気がついた・・・

（注）答えが見つからないと、途中で諦めて終了する場合があります。
（注）答えが見つからないと、無限ループする場合もあります。



--^ root.pl
#!usr/local/bin/perl
use strict;

my $f       = shift;
if(not $f =~ m/\$x\b/){
    die qq(bad func : "$f"\n);
}
my $yabserr = shift || 1e-6;
my $xrelerr = shift || 1e-6;

my $fcount  = 0;
sub f{
    my $x   = shift;
    my $y   = eval $f;
    ++$fcount;
#   printf "*** %2d (%+lf %+lf)\n",$fcount,$x,$y;
    return $y;
}

###
### あてずっぽうで y = 0 を挟む２点 (xm,ym)〜(xp,yp) を探す
###
my($xp,$yp);
$xp = 0;
$yp = f($xp);
if($@ or $yp !~ m/\d/){
    die qq(eval error : "$f"\n);
}
my($x,$y);
$x  = 0.1;
while(1){
    $y  = f($x);
    if($@){
        die qq(eval error : "$f"\n);
    }
    last if $y == 0;
    last if $y * $yp < 0;
    
    if($y and abs $y < abs $yp){
        ($xp,$yp)   = ($x,$y);
    }
    $x  *= (0 < $x) ? -1 : -2;
}
my($xm,$ym);
if($yp < 0){
    ($xm,$ym)   = ($xp,$yp);
    ($xp,$yp)   = ($x,$y);
}else{
    ($xm,$ym)   = ($x,$y);
}

###
### ２点 (xm,ym)〜(xp,yp) 間の y=0 の点を求める
###
while($yabserr < abs $y and $xrelerr < abs( ($xp-$xm)/$x) ){
#   $x  = ($xp+$xm)/2;                      # 二分法
    $x  = ($xp-$xm)/($yp-$ym)*(0-$ym)+$xm;  # 内分点
    $y  = f($x);
    if($@){
        die qq(eval error : "$f"\n);
    }
    if(0 < $y){
        ($xp,$yp)   = ($x,$y);
    }else{
        ($xm,$ym)   = ($x,$y);
    }
    if(0 < $yp * $ym){
        die qq(method failure\n);
    }
}

printf "x = %+lf\t(%+lf %+lf)\n",$x,$xm,$xp;
printf "y = %+lf\t(%+lf %+lf)\n",$y,$ym,$yp;

__END__
--$

-- 
FZH01112@..., http://www1.u-netsurf.ne.jp/~dune/