Help me with Matlab code on Newton-Raphson method to solve a system of nonlinear equations.

I want to approximate the zeros of the following nonlinear system, using N-R method:

f(x,y)=x+y^9/3+x^{243}/9+y^{2187}/27=0;

g(x,y)=y+x^{27}/3+y^{243}/9+x^{6561}/27=0.

I have tried with the following Matlab code. But I am getting error. Can you please help me out ?

close all; clc,clear

% Newton Raphson solution of two nonlinear algebraic equations

xy = [-1 -1]; % initial guesses

iter=0;

maxiter=100;

xy_N = 0.5;

TOL = 0.05;

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+y^9/3+x^{243}/9+y^{2187}/27; y+x^{27}/3+y^{243}/9+x^{6561}/27];

% 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, :));

Answer 1

[Ahmad]

why you use { } ?

this is for cell type

[Mabud Sarkar]

Thanks. I have changed the notation and used the following expression:

close all; clc,clear

% Newton Raphson solution of two nonlinear algebraic equations

xy = [-1 -1]; % initial guesses

iter=0;

maxiter=100;

xy_N = 0.5;

TOL = 0.05;

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

Can you please help me to fix the error ?

[Torsten]

The power of x and y makes the solver diverge. Use a smaller initial guess, e.g. xy = [0 0] :-), and a smaller tolerance, e.g. TOL = 1e-8.

[Mabud Sarkar]

Thanks. Yes it converged the code with initial guess (0,0) and smaller tolerance but shows peculiar behavior, it converges to (0,0), which is trivial. I am searching other roots. Someone suggested to use 'fsolves'. There is an answer https://in.mathworks.com/matlabcentral/answers/372827-fsolve-multi-variables-help. Can you please implement it ?

[Torsten]

I suggest you substitute

x~ = x^6561
y~ = y^2187

in your function.

Then the NaN problems will vanish.

[Mabud Sarkar]

Ok, But how can I calculate the Jacobian, then ? Because I can not do derivative of x~ or y~. Or am I mis interpratating your comment ?

[Torsten]

f(iter,:) = [x^(1/6561)+(1/3)*y^(9/2187)+(1/9)*x^(243/6561)+(1/27)*y^(2187/2187); y^(1/2187)+(1/3)*x^(27/6561)+(1/9)*y^(243/2187)+(1/27)*x^(6561/6561)];

Now build the Jacobian.

[Alex Sha]

[x,y]=[0,0] is the exact and only solution

[Mabud Sarkar]

No there are more 3 zeros nearest to (-1,1), (1,-1) and (-1,-1). Please see the answer of @Cesareo here