Message Boards Message Boards

GROUPS:

Improve this relatively simple FindRoot code?

Posted 7 months ago
775 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

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 7 months 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

Posted 7 months 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

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract