読者です 読者をやめる 読者になる 読者になる

前進差分と中心差分で微分

perl 数値計算

次の関数を微分するプログラムをperlで書いてみました。

  • f(x)=2x^2+4
  • f(u)=\sqrt{\frac{1}{(1+e^{-2u})}
use strict;
use warnings;

difference();

sub f1 {
    my ($x) = @_;
    my $a = 2;
    my $b = 4;
    return $a*($x ** 2)+$b;
}

sub f2 {
    my ($u) = @_;
    return sqrt(1/(1+exp(-2*$u)));
}


sub difference {
    my $x = 2;
    my $hf = 2 ** -26;
    my $hc = 2 ** -52;

    my $fdiff1 = 0;
    my $cdiff1 = 0;
    print "f1(x) = " . f1($x) . "\n";
    print "f1(x+hf) = " . f1($x+$hf) . "\n";
    $fdiff1 = ( f1($x+$hf) - f1($x) ) / $hf;
    print "fdiff1 = " . $fdiff1 . "\n";
    $cdiff1 = ( f1($x+$hc) - f1($x-$hc) ) / (2*$hc);
    print "f1(x+hc) = " . f1($x+$hc) . "\n";
    print "f1(x-hc) = " . f1($x-$hc) . "\n";
    print "cdiff1 = " . $cdiff1 . "\n";
    print "TrueDiff = " . 4*$x . "\n";


    my $fdiff2 = 0;
    my $cdiff2 = 0;
    print "f2(x) = " . f2($x) . "\n";
    print "f2(x+hf) = " . f2($x+$hf) . "\n";
    $fdiff2 = ( f2($x+$hf) - f2($x) ) / $hf;
    print "fdiff2 = " . $fdiff2 . "\n";
    $cdiff2 = ( f2($x+$hc) - f2($x-$hc) ) / (2*$hc);
    print "f2(x+hc) = " . f2($x+$hc) . "\n";
    print "f2(x-hc) = " . f2($x-$hc) . "\n";
    print "cdiff2 = " . $cdiff2 . "\n";
    print "TrueDiff = " . exp(-2*$x) / (sqrt((1 + (exp(-2*$x))) ** 3)) . "\n";
}

結果はこちら。

f1(x) = 12
f1(x+hf) = 12.0000001192093
fdiff1 = 8
f1(x+hc) = 12
f1(x-hc) = 12
cdiff1 = 4
TrueDiff = 8
f2(x) = 0.99096608924721
f2(x+hf) = 0.990966089512804
fdiff2 = 0.017823725938797
f2(x+hc) = 0.99096608924721
f2(x-hc) = 0.99096608924721
cdiff2 = 0
TrueDiff = 0.0178237241465131


関数の性質と桁落ちによりうまく求められていない部分があります。