# Steffensen's Method Implementation in Mathematica

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 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]];


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?

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.

• 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?
– Dan
Sep 18 at 16:35
• 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. Sep 18 at 16:43
• 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. Sep 18 at 16:45
• 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. Sep 18 at 16:46
• You're welcome, Sep 19 at 6:06