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?