I am trying to solve Nonlinear systems of equations using Newton-Raphson method, having problem while using newtmult function. Please guide me to correct my code.

function [J,f]=jfreact(x,varargin)

del=0.000001;

df1dx1=(u(x(1)+del*x(1),x(2))-u(x(1),x(2)))/(del*x(1));

df1dx2=(u(x(1),x(2)+del*x(2))-u(x(1),x(2)))/(del*x(2));

df2dx1=(v(x(1)+del*x(1),x(2))-v(x(1),x(2)))/(del*x(1));

df2dx2=(v(x(1),x(2)+del*x(2))-v(x(1),x(2)))/(del*x(2));

J=[df1dx1 df1dx2;df2dx1 df2dx2];

f1=u(x(1),x(2));

f2=v(x(1),x(2));

f=[f1;f2];

function f=u(x,y)

f = (5 + x + y) / (50 - 2 * x - y) ^ 2 / (20 - x) - 0.0004;

function f=v(x,y)

f = (5 + x + y) / (50 - 2 * x - y) / (10 - y) - 0.037;

solver function used:- (Initial guesses are x1=x2=3)

function [x,f,ea,iter]= newtmult(@jfreact,[3;3],es,maxit,varargin);

es=50;

if nargin<2,error('at least 2 input arguments required'),end

if nargin<3||isempty(es),es=0.0001;end

if nargin<4||isempty(maxit),maxit=50;end

iter = 0;

x=x0;

while (1)

[J,f]=jfreact(x,varargin{:});

dx=J\f;

x=x-dx;

iter = iter + 1;

ea=100*max(abs(dx./x));

if iter>=maxit||ea<=es, break, end

end

end

Answer 1

[arushi]

Hi Dhiraj,

Here's the corrected and complete code:

function newton_solver()
    % Initial guesses
    x0 = [3; 3];
    es = 0.0001;  % Desired relative error
    maxit = 50;   % Maximum number of iterations

    % Call the Newton-Raphson solver
    [x, f, ea, iter] = newtmult(@jfreact, x0, es, maxit);

    % Display results
    fprintf('Solution found after %d iterations:\n', iter);
    fprintf('x1 = %f, x2 = %f\n', x(1), x(2));
    fprintf('Function values at solution: f1 = %f, f2 = %f\n', f(1), f(2));
    fprintf('Approximate relative error: %f\n', ea);
end

function [J, f] = jfreact(x, varargin)
    del = 0.000001;
    % Calculate partial derivatives for the Jacobian matrix
    df1dx1 = (u(x(1) + del * max(x(1), 1e-6), x(2)) - u(x(1), x(2))) / (del * max(x(1), 1e-6));
    df1dx2 = (u(x(1), x(2) + del * max(x(2), 1e-6)) - u(x(1), x(2))) / (del * max(x(2), 1e-6));
    df2dx1 = (v(x(1) + del * max(x(1), 1e-6), x(2)) - v(x(1), x(2))) / (del * max(x(1), 1e-6));
    df2dx2 = (v(x(1), x(2) + del * max(x(2), 1e-6)) - v(x(1), x(2))) / (del * max(x(2), 1e-6));
    J = [df1dx1 df1dx2; df2dx1 df2dx2];

    % Evaluate functions
    f1 = u(x(1), x(2));
    f2 = v(x(1), x(2));
    f = [f1; f2];
end

function f = u(x, y)
    f = (5 + x + y) / (50 - 2 * x - y)^2 / (20 - x) - 0.0004;
end

function f = v(x, y)
    f = (5 + x + y) / (50 - 2 * x - y) / (10 - y) - 0.037;
end

function [x, f, ea, iter] = newtmult(jfhandle, x0, es, maxit, varargin)
    if nargin < 2
        error('At least 2 input arguments required');
    end
    if nargin < 3 || isempty(es)
        es = 0.0001;
    end
    if nargin < 4 || isempty(maxit)
        maxit = 50;
    end

    iter = 0;
    x = x0;
    while true
        [J, f] = feval(jfhandle, x, varargin{:});
        dx = J \ f;
        x = x - dx;
        iter = iter + 1;
        ea = 100 * max(abs(dx ./ x));
        if iter >= maxit || ea <= es
            break;
        end
    end
end

Hope this helps.