I want to approximate the zeros of the following system of polynomial equations using the Newton-Raphson method:
\begin{align} f(x,y) &= x + \frac13 y^9 + \frac19 x^{243} + \frac{1}{27} y^{2187} =0\\ g(x,y) &= y + \frac13 x^{27} + \frac19 y^{243} + \frac{1}{27} x^{6561} =0 \end{align}
I have used the following Matlab code with initial guess $(-1,-1)$ and with error tolerance $0.01$:
close all; clc,clear
% Newton Raphson solution of two nonlinear algebraic equations
xy = [-1 -1]; % initial guesses
iter=0;
maxiter=100;
xy_N = 0.1;
TOL = 0.01;
error1 = [1, 1];
f = zeros(iter, 2);
error = zeros(iter, 2);
% begin iteration
while xy_N>TOL
iter=iter+1;
x = xy(1);
y = xy(2);
% calculate the functions
f(iter,:) = [x+(1/3)*y^9+(1/9)*x^243+(1/27)*y^2187;
y+(1/3)*x^27+(1/9)*y^243+(1/27)*x^6561];
% calculate the Jacobian
J= [1+27*x^242, 3*y^8+81*y^2186; 9*x^26+24*x^6560, 1+27*y^242];
xy_new = xy - ((J)\f(iter,:)')';
error1=abs(xy_new-xy);
error(iter,:)=error1;
% Calculate norm:
xy_N = norm(f(iter,:));
XYN(iter) = xy_N;
xy=xy_new;
%iter=iter+1;
if (iter > maxiter)
fprintf('Did not converge! \n', maxiter);
break
end
end
if iter<maxiter
f(iter+1:end,:)=[];
end
% Print results
fprintf('Number of iterations is: %f \n ', iter)
fprintf('Final values of [x, y] = %f %f \n', xy(end, :));
fprintf('Final values of EQN: %f %f \n', f(end, :));
fprintf('Final values of ERROR: %f %f \n', error(end, :));
But it is giving error. The following display:
Warning: Matrix is singular to working precision. Warning: Matrix is singular, close to singular or badly scaled.
Results may be inaccurate. RCOND = NaN.
Number of iterations is: 4.000000
Final values of [x, y] = NaN NaN
Final values of EQN: NaN NaN
Final values of ERROR: NaN NaN
Could you please help me fix the code to ensure that convergence is attained? You can suggest more flexible code to solve my system.
Edit: Since the Jacobian is not invertible, I am adding one more term to each function in order to check my program as follows:
\begin{align} f(x,y) &= x + \frac13 y^9 + \frac19 x^{243} + \frac{1}{27} y^{2187} +\frac{1}{81}x^{59049}=0\\ g(x,y) &= y + \frac13 x^{27} + \frac19 y^{243} + \frac{1}{27} x^{6561} +\frac{1}{81}y^{59049}=0 \end{align} I hope now the Jacobian will be non-singular.
Problem seems to be that the Jacobian is not invertible in the fourth iteration. This is because I think the NR-method wants to converge to the solution $x:= (-1.0005, 1.0014)^\top$ but $J(x)$ is not invertible. The convergence results for the NR-method usually all apply to non-degenerate solutions $y$, i.e. $\det{J(y)} \neq 0$ (of course you cannot check that in advance).
One suggestion that comes to mind is just to store the iterations in which $\texttt{J}$ is invertible and use that old $\texttt{J}$ when it is not invertible in the current iteration (so called frozen NR-method). Even though I do not know for certain if this converges, chance is certainly higher.
(Fun fact: This would actually converge locally - though not as fast as normal NR-method - if $J(x)$ was invertible)
Another one is using the pseudo-inverse of $J$ instead of the $\texttt{\\}$ operator.