Multidimensional Newton-Raphson method: Difference between revisions
Multidimensional Newton-Raphson method (view source)
Revision as of 09:48, 29 July 2019
, 4 years ago→{{header|Phix}}
(julia example) |
|||
Line 655:
<pre>
Solution: [0.8936282344764825 0.8945270103905782 -0.04008928615915281]
</pre>
=={{header|Phix}}==
{{trans|Go}}
Uses code from [[Reduced_row_echelon_form#Phix]],
[[Gauss-Jordan_matrix_inversion#Phix]],
[[Matrix_transposition#Phix]], and
[[Matrix_multiplication#Phix]]<br>
See std distro for a complete runnable version.
<lang Phix>-- demo\rosetta\Multidimensional_Newton-Raphson_method.exw
function solve(sequence fs, jacob, guesses)
integer size := length(fs),
maxIter := 12,
iter := 0
sequence gu1, g, t, f, g1,
gu2 := guesses,
jac := repeat(repeat(0,size),size)
atom tol := 1e-8
while true do
gu1 = gu2
g := matrix_transpose({gu1})
t := repeat(0, size)
for i=1 to size do
t[i] = call_func(fs[i],{gu1})
end for
f := matrix_transpose({t})
for i=1 to size do
for j=1 to size do
jac[i][j] = call_func(jacob[i][j],{gu1})
end for
end for
g1 := sq_sub(g,matrix_mul(inverse(jac),f))
gu2 = vslice(g1,1)
iter += 1
bool any := false
for i=1 to length(gu2) do
if abs(gu2[i])-gu1[i] > tol then
any = true
exit
end if
end for
if not any or iter >= maxIter then exit end if
end while
return gu2
end function
function f1(sequence v) atom {x,y} = v return -x*x+x+0.5-y end function
function f2(sequence v) atom {x,y} = v return y+5*x*y-x*x end function
function f3(sequence v) atom {x,y,z} = v return 9*x*x+36*y*y+4*z*z-36 end function
function f4(sequence v) atom {x,y,z} = v return x*x-2*y*y-20*z end function
function f5(sequence v) atom {x,y,z} = v return x*x-y*y+z*z end function
function j1(sequence v) atom {x} = v return -2*x+1 end function
function j2(sequence /*v*/) return -1 end function
function j3(sequence v) atom {x,y} = v return 5*y-2*x end function
function j4(sequence v) atom {x} = v return 1+5*x end function
function j11(sequence v) atom {x} = v return 18*x end function
function j12(sequence v) atom {?,y} = v return 72*y end function
function j13(sequence v) atom {?,?,z} = v return 8*z end function
function j21(sequence v) atom {x} = v return 2*x end function
function j22(sequence v) atom {?,y} = v return -4*y end function
function j23(sequence /*v*/) return -20 end function
function j31(sequence v) atom {x} = v return 2*x end function
function j32(sequence v) atom {?,y} = v return -2*y end function
function j33(sequence v) atom {?,?,z} = v return 2*z end function
procedure main()
sequence fs, jacob, guesses
/*
solve the two non-linear equations:
y = -x^2 + x + 0.5
y + 5xy = x^2
given initial guesses of x = y = 1.2
Example taken from:
http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
Expected results: x = 1.23332, y = 0.2122
*/
fs = {routine_id("f1"),routine_id("f2")}
jacob = {{routine_id("j1"),routine_id("j2")},
{routine_id("j3"),routine_id("j4")}}
guesses := {1.2, 1.2}
printf(1,"Approximate solutions are x = %.7f, y = %.7f\n\n", solve(fs, jacob, guesses))
/*
solve the three non-linear equations:
9x^2 + 36y^2 + 4z^2 - 36 = 0
x^2 - 2y^2 - 20z = 0
x^2 - y^2 + z^2 = 0
given initial guesses of x = y = 1.0 and z = 0.0
Example taken from:
http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
*/
fs = {routine_id("f3"), routine_id("f4"), routine_id("f5")}
jacob = {{routine_id("j11"),routine_id("j12"),routine_id("j13")},
{routine_id("j21"),routine_id("j22"),routine_id("j23")},
{routine_id("j31"),routine_id("j32"),routine_id("j33")}}
guesses = {1, 1, 0}
printf(1,"Approximate solutions are x = %.7f, y = %.7f, z = %.7f\n", solve(fs, jacob, guesses))
end procedure
main()</lang>
{{out}}
<pre>
Approximate solutions are x = 1.2333178, y = 0.2122450
Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.04008929
</pre>
| |||