# 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 (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? – J. M. will be back soon Jan 15 '17 at 17:38
• I think I will put a similar question on the stack-exchange math site. – Ted Ersek 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. – Daniel Lichtblau 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. – Daniel Lichtblau Jan 15 '17 at 20:29