This is sort of a math and Mathematica Language question since I do not know which one is going wrong.

I am trying to implement Steffensen's Method for Nonlinear Systems of Equations and the first step is to approximate the Jacobian Matrix. According to this paper, https://downloads.hindawi.com/journals/jam/2014/705375.pdf, the formula is Jacobian Approximation Now here is my code:

a = 6;
g = 1;
V[x_, y_] = 
 Expand[1/6*(3*x^2 + 3*y^2 + 6*x^2*y - 2*y^3) + 
   a*(2*x^6*(y + 1) + x^4*y^2*(2*y + 1) + x^2*y^4*(2*y + 1) + 
      2*y^6*(y + 1)) + 
   g*(x^4*(y - 1) + x^2*(y - 2)*y^2 - y^4*(y + 1))];
Vx[x_, y_] = Expand[Derivative[1, 0][V][x, y]];
Vy[x_, y_] = Expand[Derivative[0, 1][V][x, y]];
J[x_, y_] = ( {
    {Vx[x + Vx[x, y], y] - Vx[x, y], Vx[x, y + Vy[x, y]] - Vx[x, y]},
    {Vy[x + Vx[x, y], y] - Vy[x, y], Vy[x, y + Vy[x, y]] - Vy[x, y]}
   } ) . ( {
    {1/Vx[x, y], 0},
    {0, 1/Vy[x, y]}
   } ); (*Approximated Jacobian*)
Jacob[x, y] = ( {
   {Derivative[1, 0][Vx][x, y], Derivative[0, 1][Vx][x, y]},
   {Derivative[1, 0][Vy][x, y], Derivative[0, 1][Vy][x, y]}
  } );(*Real Jacobian*)
MatrixForm[Jacob[1, 1]];
MatrixForm[J[1, 1]];

My Code Vx = 0 and Vy = 0 are the two nonlinear equations I am trying to solve. At the bottom of my code is my confusion. The approximated Jacobian is way off and I can't figure out why. For other values of Vx and Vy, the approximated Jacobian is not that off but for these functions, it is. Am I doing the Mathematica code wrong?

  • 1
    $\begingroup$ Please edit your question, and copy the code from your notebook to your question. Posting the code instead of an image makes it easier for us to help you. $\endgroup$
    – creidhne
    Sep 18, 2021 at 6:09

1 Answer 1


I think you did everything right. The problem is that {Vx[1,1],Vy[1,1]} equals {253, 292} which is very far from {0,0}. See, the entries of J are finite difference approximations of derivates. For example, the upper left entry of J is

(Vx[x + Vx[x, y], y] - Vx[x, y]) / Vx[x, y]

Replacing the Vx[x, y] that originated from the diagonal matrix by h, this becomes

(Vx[x + h, y] - Vx[x, y]) / h

which approximates

D[ Vx[x, y], x]

only for $h \to 0$. The idea of Steffensen's method seems to be that in the course of the Newton iterations with this fake Jacobian J, Vx and Vy will converge to 0, so that the fake Jacobian J approximates the true Jacobian better and better over time. This suffices as it is well-known that Newton's method converges quadratically, even if we perturb the Jacobian by something that converges sufficiently fast to the true Jacobian at the solution.

But this may also mean, that Steffensen's method won't work if one starts too far away from any solution. In that respect, it seems to be more of a diva as Newton's method already is...

Some remedy could be to Clip h to smaller values like this:

eps = 10^-8;

J[x_, y_] := Module[{Valx, Valy, hx, hy},
  Valx = Vx[x, y];
  Valy = Vy[x, y];
  hx = Clip[Valx, {-eps, eps}];
  hy = Clip[Valy, {-eps, eps}];
   {(Vx[x + hx, y] - Valx)/hx, (Vx[x, y + hy] - Valx)/hy},
   {(Vy[x + hx, y] - Valy)/hx, (Vy[x, y + hy] - Valy)/hy}

But have not tried it, yet. I also don't know which value of eps would be most suitable.

  • $\begingroup$ Thanks for the answer. I checked some values and saw that there are only a small range of values where the method converges. Are there any other methods that are more versatile? $\endgroup$
    – Dan
    Sep 18, 2021 at 16:35
  • $\begingroup$ Hard to tell. Typically it depends very much on the actual problem that you want to solve. A good idea in general is to use a line search to what people from optimization call "globalize convergence". Also Broyden's method, (a quasi-Newton method) comes to mind. $\endgroup$ Sep 18, 2021 at 16:43
  • 1
    $\begingroup$ In your special situation you want to find a the zero of a derivative of the scalar function V. So I guess you want to minimize or maximize V. In that case, I'd suggest L-BFGS. It is a special case of Broyden's method -- that works by far better than Broyden's method for general nonlinear equations. And the best of it: You do not require any second derivatives of V, nor do you have to solve any linear equation. L-BFGS maintains an appropriate approximation of the inverse of the Hessian and updates it in every iteration. $\endgroup$ Sep 18, 2021 at 16:45
  • 2
    $\begingroup$ L-BFGS does not converge quadratically, but still superlinearly. Because one spares the $O(N^3)$ costs for the linear solve, it is usually still considerably faster than Newton's method. $\endgroup$ Sep 18, 2021 at 16:46
  • $\begingroup$ You're welcome, $\endgroup$ Sep 19, 2021 at 6:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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