Multidimensional Newton-Raphson method: Difference between revisions
Multidimensional Newton-Raphson method (view source)
Revision as of 16:07, 14 March 2020
, 4 years agoRename Perl 6 -> Raku, alphabetize, minor clean-up
Langurmonkey (talk | contribs) m ("metod" to "method") |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 5:
using Newton-Raphson method.
<br><br>
=={{header|C sharp|C#}}==
For matrix inversion and matrix and vector definitions - see C# source from [[Gaussian elimination]]
<lang csharp>
Line 351 ⟶ 352:
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
</pre>
=={{header|Julia}}==
Line 387:
* Jacobian Calls (df/dx): 4
</pre>
=={{header|Kotlin}}==
Line 579 ⟶ 577:
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
</pre>
Line 761 ⟶ 683:
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.04008929
</pre>
=={{header|Raku}}==
(formerly Perl 6)
<lang perl6>#!/usr/bin/env perl6
# Reference:
# https://github.com/pierre-vigier/Perl6-Math-Matrix
# Mastering Algorithms with Perl
# By Jarkko Hietaniemi, John Macdonald, Jon Orwant
# Publisher: O'Reilly Media, ISBN-10: 1565923987
# https://resources.oreilly.com/examples/9781565923980/blob/master/ch16/solve
use v6;
sub solve_funcs ($funcs, @guesses, $iterations, $epsilon) {
my ($error_value, @values, @delta, @jacobian); my \ε = $epsilon;
for 1 .. $iterations {
for ^+$funcs { @values[$^i] = $funcs[$^i](|@guesses); }
$error_value = 0;
for ^+$funcs { $error_value += @values[$^i].abs }
return @guesses if $error_value ≤ ε;
for ^+$funcs { @delta[$^i] = -@values[$^i] }
@jacobian = jacobian $funcs, @guesses, ε;
@delta = solve_matrix @jacobian, @delta;
loop (my $j = 0, $error_value = 0; $j < +$funcs; $j++) {
$error_value += @delta[$j].abs ;
@guesses[$j] += @delta[$j];
}
return @guesses if $error_value ≤ ε;
}
return @guesses;
}
sub jacobian ($funcs is copy, @points is copy, $epsilon is copy) {
my ($Δ, @P, @M, @plusΔ, @minusΔ);
my Array @jacobian; my \ε = $epsilon;
for ^+@points -> $i {
@plusΔ = @minusΔ = @points;
$Δ = (ε * @points[$i].abs) || ε;
@plusΔ[$i] = @points[$i] + $Δ;
@minusΔ[$i] = @points[$i] - $Δ;
for ^+$funcs { @P[$^k] = $funcs[$^k](|@plusΔ); }
for ^+$funcs { @M[$^k] = $funcs[$^k](|@minusΔ); }
for ^+$funcs -> $j {
@jacobian[$j][$i] = (@P[$j] - @M[$j]) / (2 * $Δ);
}
}
return @jacobian;
}
sub solve_matrix (@matrix_array is copy, @delta is copy) {
# https://github.com/pierre-vigier/Perl6-Math-Matrix/issues/56
{ use Math::Matrix;
my $matrix = Math::Matrix.new(@matrix_array);
my $vector = Math::Matrix.new(@delta.map({.list}));
die "Matrix is not invertible" unless $matrix.is-invertible;
my @result = ( $matrix.inverted dot $vector ).transposed;
return @result.split(" ");
}
}
my $funcs = [
{ 9*$^x² + 36*$^y² + 4*$^z² - 36 }
{ $^x² - 2*$^y² - 20*$^z }
{ $^x² - $^y² + $^z² }
];
my @guesses = (1,1,0);
my @solution = solve_funcs $funcs, @guesses, 20, 1e-8;
say "Solution: ", @solution;
</lang>
{{out}}
<pre>
Solution: [0.8936282344764825 0.8945270103905782 -0.04008928615915281]
</pre>
| |||