# 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[#[], #[]] &, {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[#[], #[]] &, {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? Answer
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). Answer
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 Answer
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 Answer
Posted 1 year ago
 Perhaps you've found these tutorials, but in case you haven't: Answer