Multidimensional Newton-Raphson method: Difference between revisions
Content added Content deleted
Langurmonkey (talk | contribs) m ("metod" to "method") |
Thundergnat (talk | contribs) (Rename Perl 6 -> Raku, alphabetize, minor clean-up) |
||
Line 5: | Line 5: | ||
using Newton-Raphson method. |
using Newton-Raphson method. |
||
<br><br> |
<br><br> |
||
=={{header|C#}}== |
|||
=={{header|C sharp|C#}}== |
|||
For matrix inversion and matrix and vector definitions - see C# source from [[Gaussian elimination]] |
For matrix inversion and matrix and vector definitions - see C# source from [[Gaussian elimination]] |
||
<lang csharp> |
<lang csharp> |
||
Line 351: | Line 352: | ||
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893 |
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893 |
||
</pre> |
</pre> |
||
=={{header|Julia}}== |
=={{header|Julia}}== |
||
Line 387: | Line 387: | ||
* Jacobian Calls (df/dx): 4 |
* Jacobian Calls (df/dx): 4 |
||
</pre> |
</pre> |
||
=={{header|Kotlin}}== |
=={{header|Kotlin}}== |
||
Line 579: | Line 577: | ||
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893 |
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893 |
||
</pre> |
|||
=={{header|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> |
</pre> |
||
Line 761: | Line 683: | ||
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.04008929 |
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> |
</pre> |
||