# Improve this relatively simple FindRoot code?

Posted 1 year ago
1817 Views
|
4 Replies
|
3 Total Likes
|
 I suspect the answer to this question is likely to be embarrassingly simple, but right now I'm at a loss as to why the following two examples behave so differently: Trying to replicate the FindRoot function behavior in Mathematica for solving a system of two transcendental equations. I've setup two sample problems and I'm using the FindRoot function to check the output of my code. This is the first example and it works as expected, matching exactly FindRoot for every iteration cycle: U[x_, y_] := x^3 - 3 x y^2 + 1; V[x_, y_] := 3 x^2 y - y^3; Ux[x_, y_] = D[U[x, y], x]; Uy[x_, y_] = D[U[x, y], y]; Vx[x_, y_] = D[V[x, y], x]; Vy[x_, y_] = D[V[x, y], y]; J[x_, y_] = Ux[x, y] Vy[x, y] - Uy[x, y] Vx[x, y] // Simplify; Iterate[x_, y_] := {x - (U[x, y] Vy[x, y] - V[x, y] Uy[x, y])/J[x, y],y - (Ux[x, y] V[x, y] - U[x, y] Vx[x, y])/J[x, y]} NestList[Iterate[#[[1]], #[[2]]] &, {10., 10.}, 10] // TableForm FindRoot[{U[x, y] == 0, V[x, y] == 0}, {{x, 10.}, {y, 10.}},StepMonitor :> Print[x, " ", y],Method -> {"Newton", "UpdateJacobian" -> 1}] This produces the following output: 10. 10. 6.66667 6.66833 4.44445 4.4493 2.96297 2.97463 1.97539 2.002 1.31749 1.3768 0.882366 1.00957 0.613064 0.85679 0.505639 0.855438 0.499855 0.866033 0.5 0.866025 6.66667 6.66833 4.44445 4.4493 2.96297 2.97463 1.97539 2.002 1.31749 1.3768 0.882366 1.00957 0.613064 0.85679 0.505639 0.855438 0.499855 0.866033 0.5 0.866025 0.5 0.866025 Now testing the same code with a different set of equations: SampleParams = {x1 -> 2., y1 -> 5., x2 -> 4., y2 -> 2, x3 -> 8, y3 -> 7}; U[x_, y_] := y1 - y3 - x Cosh[(-y + x1)/x] + x Cosh[(-y + x3)/x] /. SampleParams V[x_, y_] := y2 - y3 - x Cosh[(-y + x2)/x] + x Cosh[(-y + x3)/x] /. SampleParams Ux[x_, y_] = D[U[x, y], x]; Uy[x_, y_] = D[U[x, y], y]; Vx[x_, y_] = D[V[x, y], x]; Vy[x_, y_] = D[V[x, y], y]; J[x_, y_] = Ux[x, y] Vy[x, y] - Uy[x, y] Vx[x, y] // Simplify; Iterate[x_, y_] := {x - (U[x, y] Vy[x, y] - V[x, y] Uy[x, y])/J[x, y],y - (Ux[x, y] V[x, y] - U[x, y] Vx[x, y])/J[x, y]} NestList[Iterate[#[[1]], #[[2]]] &, {10., 10.}, 3] // TableForm FindRoot[{U[x, y] == 0, V[x, y] == 0}, {{x, 10.}, {y, 10.}},StepMonitor :> Print[x, " ", y], Method -> {"Newton", "UpdateJacobian" -> 1}] This time I'm not getting the desired results: 10. 10. -61.2969 -34.1643 -2953.96 -1832.23 -6.67248*10^6 -4.14752*10^6 7.00846 8.14691 5.76637 7.37079 5.04881 6.91975 4.58046 6.62407 1.31254 4.55446 1.46434 4.67168 1.51128 4.71681 1.51456 4.72029 1.51457 4.72031 1.51457 4.72031 {x->1.51457,y->4.72031} As you can see my code simply does not work in this case, even though a root clearly exists @ {x->1.51457,y->4.72031} Any thoughts?
4 Replies
Sort By:
Posted 1 year ago
 I think I may have found a possible answer to my question: My code is a pure Newton Rapshon method whereas Mathematica actually uses a dampened version of this method with some builtin intelligence as to how to pick a damping factor to prevent the method from diverging in certain cases Now I need to figure out what that extra intelligence is
Posted 1 year ago
 You can't beat good starting values. Using {10.,10.} is a bit too far away. Try {2.,6.}. (And your first example has 3 roots.) This doesn't match with each iteration but both do end up with the same result (when more iterations are added to the NestList command).