[B] Das Newton-Verfahren in Perl

[Bonkers]

Hey,
Ich musste im Zuge eines anderen Skriptes das Newton-Verfahren Anwenden, um die Nullstellen einer Funktion zu errechnen.
Hab das ganze jetzt mal in ein eigenständiges Skript gepackt.
Das Newton-Verfahren nimmt man (glaube Ich), in der 11. Klasse durch.
Der erste Sprung über die X-Achse wird automatisch berechnet, um die benötigten Schritte stark einzuschränken!

Auch wenn ihr das Skript an sich nicht braucht freue ich mich über jede Kritik.

Usage:

Das erste Argument gibt die Anzahl an Schritten an. 0 ist default.
Wenn ihr 0 angebt, dann stoppt das Skript automatisch wenn die Nullstelle gefunden wurde! Deshalb ist 0 absolut zu empfehlen.
Die Funktion und die Ableitung sind fest im Code verankert, können aber in den 2 Subs beliebig verändert werden. Achtet darauf dass " ** " einem " ^ " entspricht

Im Skript ist

return $x**3+2*$x-1;

also

f(x) = x^3+2*x-1

So, hier das Skript. Viel Spaß!

our $start = -50; #
our $count = 100; # Range to test for switch on x-axis
our $buffer;


sub f {  # Your f(x). ** = ^
my $x = shift;
return $x**3+2*$x-1;
}
sub abl { # Derivation of f(x). ** = ^
my $x = shift;
return 3*$x**2+2;
}


print "\n".("\t"x3)."### Newton's method by Bonkers ### \n\n";
print "\n\tSearching for switch on x-axis (From ".$start." to ".$count.")...\n\n";

sub isPos {  # Is the value positive?
my $cc = shift;
#if ($cc > 0) { return 1 } else { return 0 };
return $cc > 0;
}

sub output { # To format Output
my $work = shift;
if (length($work) > 9) { return substr($work, 0,10); } else { return $work.(" "x(10-length($work))) }
}

for (my $i = $start; $i < ($start+$count); $i++) {
$buffer = f($i);
my $sec = f($i+1);
if (isPos($buffer) != isPos($sec)) {$buffer = $i - 1; $i = ($start+$count); }
}
print "\tDone. Selected Value to start with: xn = ".$buffer."\n\n";
our $maximum = $ARGV[0];
if ($maximum != 0) { if ($maximum == "") { $maximum = 100; }}
print "\tMaximum steps, taken from first Argument (Default is 0) : ".$maximum."\n\tUse 0 for Auto-Stop.\n\n\nPress Enter to start. "; <STDIN>;
print "\n\n\t ".("_"x62)."\n\t| Xn           | f(xn)         |  f'(xn)      |  1-xn         |"."\n\t|".("-"x62)."\n";
print "\t| ".output($buffer)."   |   ".output(f($buffer))."  |  ".output(abl($buffer))."  |  ".output(1-(f($buffer)/abl($buffer)))."   |   \n";


for (my $counter = -1; $counter < $maximum;) {
$buffer = 1-(f($buffer)/abl($buffer));
print "\t| ".output($buffer)."   |   ".output(f($buffer))."  |  ".output(abl($buffer))."  |  ".output(1-(f($buffer)/abl($buffer)))."   |   \n";


if ($maximum != 0) { $counter++; } else { if ($buffer eq (1-(f($buffer)/abl($buffer)))) { $maximum = 1; $counter = 1;} }
}
print "\t-".("-"x62).("\n"x3)."\t\t Done :)\n\n";

Beispiel-Output:

[root@archBox Desktop]# perl math.pl 20

			### Newton's method by Bonkers ### 


	Searching for switch on x-axis (From -50 to 100)...

	Done. Selected Value to start with: xn = -1

	Maximum steps, taken from first Argument (Default is 100) : 20


Press Enter to start. 


	 ______________________________________________________________
	| Xn           | f(xn)         |  f'(xn)      |  1-xn         |
	|--------------------------------------------------------------
	| -1           |   -4          |  5           |  1.8          |   
	| 1.8          |   8.432       |  11.72       |  1.8          |   
	| 0.28054607   |   -0.4168271  |  2.23611830  |  0.28054607   |   
	| 1.18640657   |   3.04275227  |  6.22268171  |  1.18640657   |   
	| 0.51102235   |   0.15549504  |  2.78343152  |  0.51102235   |   
	| 0.94413548   |   1.72986563  |  4.67417546  |  0.94413548   |   
	| 0.62990999   |   0.50975982  |  3.19035979  |  0.62990999   |   
	| 0.84021870   |   1.27360447  |  4.11790240  |  0.84021870   |   
	| 0.69071523   |   0.71096210  |  3.43126261  |  0.69071523   |   
	| 0.79279869   |   1.08389495  |  3.88558929  |  0.79279869   |   
	| 0.72104747   |   0.81697434  |  3.55972836  |  0.72104747   |   
	| 0.77049531   |   0.99840520  |  3.78098908  |  0.77049531   |   
	| 0.73594073   |   0.87047341  |  3.62482628  |  0.73594073   |   
	| 0.75985789   |   0.95844559  |  3.73215206  |  0.75985789   |   
	| 0.74319224   |   0.89687536  |  3.65700413  |  0.74319224   |   
	| 0.75475134   |   0.93944647  |  3.70894875  |  0.75475134   |   
	| 0.74670815   |   0.90976067  |  3.67271921  |  0.74670815   |   
	| 0.75229234   |   0.93033984  |  3.69783130  |  0.75229234   |   
	| 0.74840933   |   0.91601509  |  3.68034958  |  0.74840933   |   
	| 0.75110649   |   0.92595796  |  3.69248290  |  0.75110649   |   
	| 0.74923161   |   0.91904292  |  3.68404405  |  0.74923161   |   
	---------------------------------------------------------------


		 Done :)

Ich muss wohl noch eine Funktion zum Runden einbauen, aber nur wenn jemand Verwendung für das Skript hat.

Changelog:

19.3.2010:

- Added 0 as Maximum to activate automatic stopping, when the point has been found. You can now use 'perl math.pl 0'.

[blackberry]

Ich habe mir noch etwas überlegt, um das eingeben von Polynomen etwas zu vereinfachen.

Ein Polynom N-ten Grades lässt sich ja wie folgt darstellen:
f(x) = (aN) * x^(N) + (aN-1) * x^(N-1) + [...] + a2 * x^2 + a1 * x + a0

Dabei muss man ja nur die Exponenten und Koeffizienten speichern.
Das macht sich relativ gut in einem Hash, weil die Exponenten jeweils nur einmal vorkommen (gesetzt man hat richtig zusammengefasst... und nicht sowas stehen wie x+x anstatt 2x).

Über Ableitungsregeln kann man dann auch automatisch die Ableitung berechnen lassen.

Beispielcode:

#!/usr/bin/perl

use strict;
use warnings;

sub print_polynom
{
    my %poly = %{$_[0]};
    my %sign = (-1 => "-", 1 => "+");
    my $i = 0;
    
    foreach my $key (reverse sort keys %poly)
    {
        my $val = $poly{$key};
        next if $val == 0;
        print " " . $sign{$val / abs $val} . " " if $i++ != 0 || $val < 0;
        print abs($val) if $key == 0;
        print abs($val) . "*x" if $key == 1;
        print abs($val) . "*x^$key" if $key != 0 && $key != 1;
    }
    print "\n";
}

sub calc_derivate
{
    my %poly = %{$_[0]};
    
    foreach my $key (keys %poly)
    {
        # Verwendete Ableitungsregeln:
        # ----------------------------
        # f(x) = a * x^n --> f'(x) = an * x^(n-1)
        # f(x) = g(x) + h(x) --> f'(x) = g'(x) + h'(x)
        $_{$key - 1} = $key * $poly{$key};
    }
    return %_;
}

sub eval_polynom
{
    my %poly = %{$_[0]};
    my $x = $_[1];
    
    $_ = 0;
    foreach my $key (keys %poly)
    {
        $_ += $poly{$key} * ($x ** $key);
    }
    return $_;
}

my %poly = (
    # f(x) = 1*x^3 + 2*x - 1
    3 => 1,
    1 => 2,
    0 => -1
    # --> f'(x) = 3*x^2 + 2
);

# Polynom ausgeben...
print "f(x) = ";
&print_polynom(\%poly);

# Ableitung berechnen und ausgeben...
my %deri = &calc_derivate(\%poly);
print "f'(x) = ";
&print_polynom(\%deri);

# f(0) berechnen und ausgeben...
print "f(0) = ";
print &eval_polynom(\%poly, 0);

print "\n";

P.S.: ich weiß, dass dein Weg mehr generisch ist, weil man da auch andere Funktionenklassen mit deren Ableitung eintragen kann, jedoch muss man dafür auch die Ableitung selber berechnen

P.S.2: ansonsten sicher eine ganz nette Idee (obwohl das Newton-Verfahren eher nutzlos ist wenn man eh am PC sitzt und einen Funktionplotter zur Hand hat...)


MfG. BlackBerry

[Bonkers]

Danke für deine Antwort!

obwohl das Newton-Verfahren eher nutzlos ist wenn man eh am PC sitzt und einen Funktionplotter zur Hand hat...

Stimmt schon, aber erstens war das Skript ja "Silent" gedacht, und zweitens musste ich damals die Tabellen per Hand im Heft rechnen, da hätte ich das Skript gebrauchen können

mit deren Ableitung eintragen kann, jedoch muss man dafür auch die Ableitung selber berechnen

Die Ableitung zu "berechnen" sollte wohl nicht das Problem sein, oder? Vielleicht sollte einer von uns beiden ja eine "Ableitungs-Sub" schreiben, sollte per RegEX kein Problem darstellen.

[blackberry]

Die Ableitung zu "berechnen" sollte wohl nicht das Problem sein, oder?

Vielleicht sollte einer von uns beiden ja eine "Ableitungs-Sub" schreiben, sollte per RegEX kein Problem darstellen.

Das wärst dann wohl eher du, da ich nur wenig Perl-Erfahrung habe (Hashes subs zu übergeben musste ich mir z.B. auch vorher nochmal anschauen...).

[Bonkers]

Das wärst dann wohl eher du, da ich nur wenig Perl-Erfahrung habe

Ich werde mich auf jeden Fall mal dran versuchen, obwohl ich irgendwie das Gefühl habe dass das trotz RegEX nicht so einfach werden wird.

Zur Perl-Erfahrung:

Egal um was es geht, ich Code immer einfach drauf los, wozu habe ich schließlich immer Google zur Hand