# How to extend a root finding method from 1D to multi dimensions

The number (0.25 Pi) is a third order root of f[x] below.

Clear[f];
f[x_] := Sin[x] - Cos[x] - Sqrt[2] (x - 0.25 Pi);
Chop@Normal@Series[f[x], {x, 0.25 Pi, 4}]


-0.235702 (-0.785398 + x)^3

Newton's method converges slowly to a multiple root such as the third order root above, but if f[x] is sufficiently well behaved, g[x_]:=If[ f[x]==0, 0, f[x]/f'[x] ] only has simple roots. I can't understand why, but very few books on numerical methods mention that. Anyway we can implement Newton's method to find a root of g[x] above and have quadratic convergence on the root. That can be done as follows:

rt=NestList[With[{y=f[#],yp=f'[#]},
If[y==0,#,#+y*yp/(y*f''[#]-yp^2)]]&,
5.0,5]


{5.,3.86977,2.28629,0.996825,0.786027,0.785398}

0.25 Pi - %


{-4.2146, -3.08437, -1.50089, -0.211427, -0.000629266, -7.64783*10^-10}

How do we extend the code above to find a root of a system of nonlinear equations in two or more dimensions? Example given:

f[x_,y_]:={Sin[x+y]+Cos[x-y]+x,Sin[x+y]+y-Cos[x-y]};


find the root of f[x,y] at { x-> -0.514933, y-> 0.514933 } . I know how to find that root with the built-in FindRoot. I want to know how to extend the method above.

• It is certainly unfortunate that Schröder's method is not more well-known; nevertheless, this is more a math question than a Mathematica question, no? Jan 15 '17 at 17:38
• I think I will put a similar question on the stack-exchange math site. Jan 15 '17 at 17:51
• Really not a Mathematica question. There is a substantial literature on root deflation, though. Macaulay matrices get used in many such, at least in the case of polynomial systems. I should add that in general this is more complicated than the univariate case, because the "multiplicity structure" can be more complicated: in the univariate case multiplicity n implies n-1 derivatives vanish at the root and these provide a basis for what is called the "dual space. Dual to what? I won't get into that. But the idea is that in the multivariate setting things need not be so straightforward. Jan 15 '17 at 20:22
• As for f[x]/f'[x] having simple roots, just factor out the root, that is, write f as f[x] = (x-r)^n*g[x] where g does not vanish at the root r. Then the quotient is (x-r) * g/(n*g+(x-r)*g' and the second factor is simply 1/n at x=r, hence this quotient vanishes to first order. Jan 15 '17 at 20:29

2D Schröder's method, which gives the idea how it can be extended to whatever dimensions. The first step is to construct the derivative tensors {zz, zp, zpp} and the iteration step.

ff[x_, y_] := {Sin[x + y] + Cos[x - y] + x,
Sin[x + y] + y - Cos[x - y]};

{zz, zp, zpp} = NestList[D[#, {{x, y}}] &, ff[x, y], 2];

Block[{LinearSolve, x, y, p},
schroederStep[p : {x_?NumericQ, y_?NumericQ}] =
p - LinearSolve[zp.zp - zpp.zz, zp.zz];
]


Try it out:

root = {x, y} /. FindRoot[ff[x, y] == 0, {x, 0}, {y, 0}];
seq = NestWhile[schroederStep, {0., 0.},
Norm[Subtract[##]] > 100 \$MachineEpsilon &,
2, 6]
seq - root
(*
{-0.514933, 0.514933}
{0., 0.}               <-- cool: no error
*)


wp = 1000;
root = {x, y} /.
FindRoot[ff[x, y] == 0, {x, 0}, {y, 0}, WorkingPrecision -> wp];
Block[{$$MinPrecision = wp,$$MaxPrecision = wp},
seq = NestList[schroederStep, {0, N[1/10, wp]}, 10]
];
Log[
Norm /@ Transpose[Transpose@seq - root]
] // Ratios // N
(*
{2.54409, 2.48968, 2.41174, 2.20723, 2.09436,
2.04506, 2.02203, 2.0109,  2.00542, 2.0027}  <-- approaches 2
*)
`