# Improve this relatively simple FindRoot code?

Posted 1 year ago
1865 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
 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).
Posted 1 year ago
 Thanks Jim, but the question here is not about how to find the roots in those particular examples, but rather - if FindRoot is using the Newton Raphson method to search for roots, what kind of tweaks is it doing that result in a different behavior in that particular case. If FindRoot manages to work properly with the initial {10,10} /I know it's far and that's on purpose/ I'd like to be able to do the same manually to replicate the exact intermediate states