matlabmultidimensional-arraymontecarlonumerical-integration

monte carlo integration on a R^5 hypercube in MATLAB


I need to write MATLAB code that will integrate over a R^5 hypercube using Monte Carlo. I have a basic algorithm that works when I have a generic function. But the function I need to integrate is:

∫dA

A is an element of R^5.

If I had ∫f(x)dA then I think my algorithm would work.

Here is the algorithm:

% Writen by Jerome W Lindsey III

clear;
n = 10000;

% Make a matrix of the same dimension
% as the problem.  Each row is a dimension

A = rand(5,n);

% Vector to contain the solution

B = zeros(1,n);


    for k = 1:n
        % insert the integrand here
        % I don't know how to enter a function {f(1,n), f(2,n), … f(5n)} that
        % will give me the proper solution
        % I threw in a function that will spit out 5!
        % because that is the correct solution.
        B(k) = 1 / (2 * 3 * 4 * 5);

    end

mean(B) 

Solution

  • %Question 2, problem set 1
    % Writen by Jerome W Lindsey III
    
    clear;
    n = 10000;
    
    % Make a matrix of the same dimension
    % as the problem.  Each row is a dimension
    A = rand(5,n);
    
    % Vector to contain the solution
    B = zeros(1,n);
    
    
        for k = 1:n
            % insert the integrand here
            % this bit of code works as the integrand
            if sum(A(:,k)) < 1
                B(k) = 1;
            end
    
        end
        clear k;
    
    clear A;
    
        % Begin error estimation calculations
        std_mc = std(B);
        clear n;
        clear B;
    
        % using the error I calculate a new random
        % vector of corect length
        N_new = round(std_mc ^ 2 * 3.291 ^ 2 * 1000000);
        A_new = rand(5, N_new);
        B_new = zeros(1,N_new);
        clear std_mc;
    
            for k = 1:N_new
                if sum(A_new(:,k)) < 1
                    B_new(k) = 1;
                end
            end
            clear k;
    
        clear A_new;
    
        % collect descriptive statisitics
        M_new = mean(B_new);
        std_new = std(B_new);
        MC_new_error_999 = std_new * 3.921 / sqrt(N_new);
        clear N_new; 
        clear B_new;
        clear std_new;
    
    % Display Results
    disp('Integral in question #2 is');
    disp(M_new);
    disp(' ');
    disp('Monte Carlo Error');
    disp(MC_new_error_999);