5
$\begingroup$

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.

$\endgroup$
4
  • 2
    $\begingroup$ 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? $\endgroup$
    – J. M.'s torpor
    Jan 15 '17 at 17:38
  • $\begingroup$ I think I will put a similar question on the stack-exchange math site. $\endgroup$
    – Ted Ersek
    Jan 15 '17 at 17:51
  • $\begingroup$ 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. $\endgroup$ Jan 15 '17 at 20:22
  • $\begingroup$ 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. $\endgroup$ Jan 15 '17 at 20:29
2
$\begingroup$

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
*)

Verify quadratic convergence:

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
*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.