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.
