# Working Iterative formula for a system of equations

I have been given a project, I need to show the use of a version of Newton's method to solve these non-linear equations. The version of Newton's method I am required to use is: $$X_{n+1} = X_n - J^{(-1)} F(X_n)$$. I have all the values required here, and this works to find the point $$X_1$$. However I need to find a code that inputs the following points $$X_2, X_3, ..., X_n$$ automatically. Here are the given values:

f[x_, y_] := x^2 + y^2 - 5;
g[x_, y_] := x^3 - y^3 - 7;

x0 = 2.1;
y0 = 0.9;

f[x0, y0]
g[x0, y0]


0.22

1.532

M = {{2*x0, 2*y0}, {3*x0^2, -3*y0^2}}


{{4.2, 1.8}, {13.23, -2.43}}

J = Inverse[M]


{{0.0714286, 0.0529101}, {0.388889, -0.123457}}

F0 = {{f[x0, y0]}, {g[x0, y0]}}
X0 = {{x0}, {y0}}
X1 = X0 - J.F0


{{0.22}, {1.532}}
{{2.1}, {0.9}}
{{2.00323}, {1.00358}}

• Have a look at FixedPoint and FixedPointList and their option SameTest. Nest(List) and NestWhile(List) also come to mind. You can also use FindRoot which has all this built-in. – Henrik Schumacher Dec 11 '18 at 12:32
• Ok Thank you, I will try some of the commands you have listed. Unfortunately, FindRoot does not help me in my situation as I need to show all the iterations leading up to the point, rather than just find the point itself. – Andrew Bradley Dec 11 '18 at 13:02
• Oh, this can also be done with FindRoot: Try Reap[FindRoot[Sin[x] == 0.2, {x, Pi/42}, EvaluationMonitor :> Sow[x]]]. Btw.: Using LinearSolve instead of Inverse should be faster for larger systems and should prevent certain problems with precision loss. – Henrik Schumacher Dec 11 '18 at 14:20

jac = D[{f[x, y], g[x, y]}, {{x, y}, 1}];
Xlist = NestList[# - Inverse[jac /. Thread[{x, y} -> #]].{f @@ #, g @@ #} &, {x0, y0}, 5]


{{2.1, 0.9}, {2.0032275, 1.0035802}, {2.0000033, 1.0000051}, {2., 1.}, {2., 1.}, {2., 1.}}

You can get the same from FindRoot:

{res, {steps}} = Reap[FindRoot[{f[x, y], g[x, y]}, {{x, x0}, {y, y0}},
Method -> "Newton", StepMonitor :> Sow[{x, y}]]]


{{x -> 2., y -> 1.}, {{{2.0032275, 1.0035802}, {2.0000033, 1.0000051}, {2., 1.}, {2., 1.}}}}

• Thank you so much, you are my hero! – Andrew Bradley Dec 11 '18 at 13:48

I'll show it with FixedPointList

f[x_, y_] = {x^2 + y^2 - 5, x^3 - y^3 - 7};
j[x_, y_] = Grad[f[x, y], {x, y}];


with LinearSolve:

FixedPointList[(# - LinearSolve[j[Sequence @@ #], f[Sequence @@ #]]) &, {2.1, 0.9}, 3]
{{2.1, 0.9}, {2.00323, 1.00358}, {2., 1.00001}, {2., 1.}}


with Inverse:

FixedPointList[(# - Inverse[j[Sequence @@ #]].f[Sequence @@ #]) &, {2.1, 0.9}, 3]
{{2.1, 0.9}, {2.00323, 1.00358}, {2., 1.00001}, {2., 1.}}