1
$\begingroup$

How do I use second order Taylor method to solve a system of non-linear equations? Is there a good reference that gives details? I found mentions of it in a dozen of numerical analysis books, but no examples

Specifically, $f:\mathbb{R}^n \to \mathbb{R}^m$, solve $f(\mathbf{x})=\mathbf{0}$ using second order Taylor expansion of $f$ around initial guess $\mathbf{x_0}$

  • 0
    Which NA books have you been looking at? Note that the multidimensional version of Newton-Raphson involves expanding $f(\mathbf{x})$ up to the Jacobian-containing term (first-order).2010-09-30
  • 0
    the ones that come up in google books when I search for "higher order Taylor"2010-09-30
  • 0
    I don't see any practical generalizations of Halley's method to multidimensional equations, if that's what you're getting at; the quadratic term involves a rank-3 tensor, and it looks unwieldy to manipulate in manner of how one would derive Halley's method from the Taylor expansion.2010-09-30
  • 0
    It looks unwieldy, which is why I'm looking for some reference that goes through the details2010-09-30
  • 0
    Actually, what I was getting at is that there's one question you have to ask first: "how does one 'invert' a rank-3 tensor?"2010-09-30

2 Answers 2

1

One way to solve $f(x) = y$ is to minimize $g(x) = (f(x) - y)^2$. You could do this by taking the quadratic approximation to $g(x)$ that comes from the second order Taylor series centered at your initial guess, finding its minimum, and letting that minimum be the starting guess for the next iteration. This isn't necessarily the best algorithm, but it's easy to understand and implement.

For high-dimensional problems, the work is in solving the system of equations to minimize the quadratic approximation. A great deal of research has gone into doing that step cleverly, taking advantage of sparse matrix structure etc.

  • 0
    Well, I don't need a high-dimensional problem, just an example for n=2, m=2. As J.M. notes, it involves rank-3 tensor2010-09-30
  • 0
    I think you and I may have different meanings of "second order" in mind. By second order, I mean up to and including 2nd derivatives. The quadratic approximation involves a linear system of equations, a rank-2 tensor. Sounds like you have in mind including third derivatives and hence encountering a rank-3 tensor.2010-09-30
  • 0
    Second derivative of f in my example is an m-by-n-by-n tensor2010-09-30
  • 0
    I still consider converting simultaneous nonlinear equations into a nonlinear least squares problem as a bit of a cheat, but... :D Anyway, what Yaroslav has originally is a $m$-dimensional vector function with $n$ independent variables, so the first term of the Taylor expansion has a vector, the second (linear) term has a matrix... you get the idea. Rather unfortunate "order" is a very overloaded mathematical term, I must say.2010-09-30
0

What you probably want is the Euler-Chebychev method for functions $F:\Bbb R^n\to\Bbb R^n$, a multi-dimensional generalization of the Halley method.

You know that $s=-F'(x)^{-1}F(x)$ gives $F(x+s)=O(s^2)$ and thus $(x+s)-x^*=O(s^2)$. Now add a further refinement $v=O(s^2)$ based on the values at the point $x$ and disregard any terms $O(s^3)$, \begin{align} F(x+s+v)&=F(x)+F'(x)(s+v)+\frac12F''(x)[s+v,s+v]+O((s+v)^3)\\ &=F'(x)v+\frac12F''(x)[s,s]+O(s^3) \end{align} Thus with $v=-\frac12F'(x)^{-1}F''(x)[s,s]$ (which is indeed $O(s^2)$) you obtain a third order method, \begin{align} s_n&=-F'(x_n)^{-1}F(x_n)\\ x_{n+1}&=x_n+s_n-\frac12F'(x_n)^{-1}F''(x_n)[s_n,s_n] \end{align}

For non-square systems replace the inverse matrix by the Moore-Penrose pseudo-inverse of $F'(x_n)$ to get cubic convergence to a point of minimal error.