Multivariate Newthon's Method

Dear All,

I am trying to solve a multivariate function by using the Newthon's Method. However, I get an error message in the lwhen I run my code (exactly in newtons_method file), and I could not figure out what is going on. Here is the files that I am using:

%This file is the one that appears here: Multivariate Newton’s Method (newtons_method_n)
function [x,k,x_all] = newtons_method_n(f,J,x0,opts)

    if (nargin < 4) || isempty(opts) || ~isfield(opts,'k_max')
        k_max = 200;
    else
        k_max = opts.k_max;
    end

    if (nargin < 4) || isempty(opts) || ~isfield(opts,'return_all')
        return_all = false;
    else
        return_all = opts.return_all;
    end

    if (nargin < 4) || isempty(opts) || ~isfield(opts,'TOL')
        TOL = 1e-10;
    else
        TOL = opts.TOL;
    end

    n = length(x0);
  
    if f(x0) == zeros(n,1)
        x = x0;
        return
    end
    
  
    x_curr = x0;
    x_next = zeros(n,1);
 
    if return_all
        x_all = zeros(n,k_max+1);
    end
  
    for k = 1:k_max
 
        if return_all
            x_all(:,k) = x_curr;
        end

        y = J(x_curr)\(-f(x_curr));
        
        x_next = x_curr+y;
        
        if (norm(y) < TOL)
        break;
        end
      
        x_curr = x_next;
        
    end
    
    x = x_next;
 
    if return_all
        x_all(:,k+1) = x;
        x_all = x_all(:,1:(k+1));
    end
    
end

This is the main.n file:

clear all
clc

% Global variables are only parameters which do not change
global a b theta t d n L Y g 

load DATA

n = 19;                              % Number of countries
a = 0.134813590592666;                
b  = 0.21221;
theta = 8.28;                        % Shape of Pareto dist
relaw = aw./aw(n,1);                 % Relative wage to US
t = exp( b*(S+theta*(log(relaw))) ); % Technology: see Table VI
d = D;                               % Geographic barriers dni^(-theta)
L = l;                               % Labor
Y = y;                               % Income
elast = 0.1;                         % Elasticity of wage adjustment

% Gamma constant (as in the original code)
g = (b^(-b))*((1-b)^(-(1-b)));                    
g = g^(theta);                                    

disp('Newton iterative procedure starts here')

w = aw;                     % Set initial wages at data
    
    % Set bounds (implied by the model) for prices
     pbounds = bounds(w);                   
    p0 = 0.5*(pbounds(:,1)+pbounds(:,2));  % Half way between min and max

   
  [pval,fval] = newtons_method(@(p) prices(p,w),p0);

And finally the prices.m file:

  function [fval,fjac] = prices(p,w)

global b theta t d n g 


aux1 = g*t.*(w.^(-theta*b));
aux2 = repmat(aux1,1,n)';
G = d.*aux2;
fval = G\p -p.^(1-b);
fjac = inv(G) - (1-b)*diag(p.^(-b));

With all this on hand, any idea about why is there an error?

Thank you in advanced!

Answer 1

[Nipun]

Hi Ibai,

I understand that you are encountering an error when running your Newton's Method code for solving a multivariate function. The error likely occurs due to a mismatch between the function signature in your main script and the actual function definition for newtons_method_n. Specifically, newtons_method_n requires both the function and its Jacobian as inputs, but you are passing only the function handle prices.

Here is the corrected code for your main script:

clear all
clc
% Global variables are only parameters which do not change
global a b theta t d n L Y g 
load DATA
n = 19;                              % Number of countries
a = 0.134813590592666;                
b  = 0.21221;
theta = 8.28;                        % Shape of Pareto dist
relaw = aw./aw(n,1);                 % Relative wage to US
t = exp( b*(S+theta*(log(relaw))) ); % Technology: see Table VI
d = D;                               % Geographic barriers dni^(-theta)
L = l;                               % Labor
Y = y;                               % Income
elast = 0.1;                         % Elasticity of wage adjustment
% Gamma constant (as in the original code)
g = (b^(-b))*((1-b)^(-(1-b)));                    
g = g^(theta);                                    
disp('Newton iterative procedure starts here')
w = aw;                     % Set initial wages at data
    
    % Set bounds (implied by the model) for prices
pbounds = bounds(w);                   
p0 = 0.5*(pbounds(:,1)+pbounds(:,2));  % Half way between min and max

opts.k_max = 200;
opts.return_all = false;
opts.TOL = 1e-10;

% Define the function and its Jacobian
f = @(p) prices(p,w);
J = @(p) jacobian_prices(p,w);

[pval, k, x_all] = newtons_method_n(f, J, p0, opts);

For more information on functions and arguments in MATLAB, refer to the following MathWorks documentation: https://www.mathworks.com/help/matlab/ref/function.html

Hope this helps.

Regards,

Nipun

[Ibai Ostolozaga Falcon]

Thank so much for your reply. Finally, I could manage to solve the problem but your answer is more than welcome! As I can observe from your respond, I've arrived to the same solution, because as you said "newtons_method_n requires both the function and its Jacobian as inputs, but you are passing only the function handle prices".

Thank you!

Ibai