How to find area enclosed by points in octave using Quadrature or any other method

Question

I have two sets of coordinates (both positive and negative values, not necessarily in increasing order, and in many cases there are different values of y for the same value of x) which I can load into two row vectors of equal size.

I want to calculate the area enclosed by the curve. How to do it with octave?

I tried this answer but it does not work because it seems that the area printed (204.64) is too high (see picture).

I tried the code:

function showdata(fName)
    M = dlmread(fName);
    H = M(2:end, 1); % starting row number is 2
    B = M(2:end, 2);
    aux = figure();

    plot(H, B,'linewidth',2);
    xlabel ("Auxilary field H (A/m)"); 
    ylabel ("Magnetic Field B (Tesla)");
    area = polyarea(H,B)

    axis([min(H), max(H), min(B), max(B)]);
    grid on;
    grid minor on;
    title (area,"fontsize",20);

Then I am calling showdata('data.txt') in Octave. Picture of Data points:

This is the data file I am using.

Answer 1

There is a function for computing convex hull called "convhull" in Octave. It returns the indices of the points formming convex hull data.

M = dlmread("data.txt"); #I removed the header in data.txt 
x = M(:,1);
y = M(:,2);
k = convhull(x,y);
plot (x(k), y(k), "r-", x, y, "b+");
n = rows(k);
x_prime = vertcat(x(k(n)), x(k(1:n-1)));
y_prime = vertcat(y(k(n)), y(k(1:n-1)));
A = .5*abs(x_prime'*y(k)-y_prime'*x(k)); #80.248
polyarea(x(k), y(k)) == A and true

Maybe convex hull is not good estimate of area because the top left and the down-right lines are a little far away from the points. There are other ways to form a polygon from data
, one of which could be alpha shape. However, alpha shape are more complicated and there is no corresponding pre-built function in Octave.

Update: Each x corresponds to at least one y cordinate. I marked the highest point and lowest point laying on the same x and estimate the area again. There is the code:

[uni, ~] = sort(unique(x));
n = rows(uni);
outline = [];
for i = 1:n
  y_list = y(x==uni(i));
  [y_max, ~] = max(y_list);
  outline(i, :)= [uni(i), y_max];
  [y_min, ~] = min(y_list);
  outline(2*n-i+1,:)= [uni(i), y_min];
endfor
figure;
plot (x(k), y(k), "r-", x, y, "b+", outline(:,1), outline(:,2), "g-", "linewidth", 3);
polyarea(outline(:,1), outline(:,2)) #74.856

By the way, if the arguments of function polyarea do not form a close curve function polyarea would return wrong area. Four point on a unit square: [(0,0), (1,0), (1,1), (0,1)], [(0,0), (1,1), (1,0), (0,1)]

polyarea([0,1,1,0],[0,0,1,1])!==polyarea([0,1,1,0],[0,1,0,1]).