multiple integration with non uniform spacing

I need help to find an area intergeral of the attached file(av_val_1) with 10 columns.

where 1st column = R , second column = C, Fifth column = F(R,C).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

R= avg_val_1(:,1);%first columsn = r axis

C = avg_val_1(:,2);%second column = c axis

F = avg_val_1(:,5);%fifth column = F(r,c)

I = trapz(C,trapz(R,F,2)); %area integeral of the total area. This line is wrong.

And

q = integral2(F(r,c),Rmin,Rmax,Cmin,Cmax) %area integeral of the particular area portion. This line is wrong. please see attached figure for reference

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Answer 1

[Ameer Hamza]

You first need to convert the data into a grid before applying the trapz twice

x = load('avg_val_1.txt');

R = x(:,1);
C = x(:,2);
F = x(:,5);

r = linspace(min(R), max(R), 1000);
c = linspace(min(C), max(C), 1000);
[Rg, Cg] = meshgrid(r, c);
Fg = griddata(R, C, F, Rg, Cg);

result = trapz(c, trapz(r, Fg, 2));

Answer 2

[MS]

Hi Ameer,

Thank you very much. I need two modfications in the code.

1, I need to apply the same code through loop and save result as seperate file.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,5);
    r = linspace(min(R), max(R), 1000);
    c = linspace(min(C), max(C), 1000);
    [Rg, Cg] = meshgrid(r, c);
    Fg = griddata(R, C, F, Rg, Cg);
    result = trapz(c, trapz(r, Fg, 2));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

2, I want to choose min(R) and max(R) values for each and every file in the for loop as shown in the atttachment figure.

thanks

[Ameer Hamza]

You can use writematrix to save these three values.

for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,5);
    r = linspace(min(R), max(R), 1000);
    c = linspace(min(C), max(C), 1000);
    [Rg, Cg] = meshgrid(r, c);
    Fg = griddata(R, C, F, Rg, Cg);
    result = trapz(c, trapz(r, Fg, 2));
    writematrix([min(r) max(r) result], ['min_max_integral_val' num2str(i)])
end

[MS]

Thanks. is it possible to save the result in a single file.

Also, Sorry for confusing you for the second part. Basically I need to choose different rectangular area for a different files. Is it possible to choose the limts(rmin,rmax,cmin, cmax) using disp command or something. Kindly modify it.

[Ameer Hamza]

To save in single file:

results = zeros(numl(avg_mat),1);
for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,5);
    r = linspace(min(R), max(R), 1000);
    c = linspace(min(C), max(C), 1000);
    [Rg, Cg] = meshgrid(r, c);
    Fg = griddata(R, C, F, Rg, Cg);
    result(i) = trapz(c, trapz(r, Fg, 2));
end
writematrix(result, 'integral_val.txt')

I am not clear about the 2nd question.

[MS]

Thanks. I am atttaching you an image. In the image you can see two rectangles in two diffrent files. I am trying to do an area integeral of those rectangles with a diffrent orientation.

In that case, the 'r' and 'c' limts changes for every file. is there a way to choose the 'r' and 'c' limts manually for each file? Kindly let me know if you need any calrifications.

Ideally the rectangle area should be chosen automatically based on the masked strucutre(nans) rotation. I know, it is difficult to automate it. It will be great if it is possible to choose the limits automatically.

[Ameer Hamza]

You can use getrect(): https://www.mathworks.com/help/images/ref/getrect.html to manually specify the rectangular coordinate on an image. I am not sure how to handle the area integration for a different orientation.

Also, automating it does not seems to be an easy task. It will require some sophisticated image analysis skills.

[MS]

Many thanks for your help. you helped me alot.

[Ameer Hamza]

Glad to be of help.

[MS]

thanks. We have done the area integeral from the identfied points [Rg, Cg] till now. is it possible to do line integeral around the rectangular line from the identified points[Rg, Cg]?

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,6);
    xy1 = [min(R), max(R)];
    xy2 = [min(C), max(C)];
    t = linspace(0,1,1000)';
    xyseg = (1-t)*xy1 + t*xy2;
    [xmesh,ymesh] = meshgrid(1:1000,1:1000);
    zseg = interp2(xmesh,ymesh,F,xyseg(:,1),xyseg(:,2)); %%this line is wrong, The number of input coordinate arrays does not equal the number of dimensions (NDIMS) of these arrays%%%%
    d = cumsum([0;sqrt(sum(diff(xyseg).^2,2))]);
    line_integeral = trapz(d,zseg)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

please insert the line it if possible.

[Ameer Hamza]

I think that xy1 and xy2 are written incorrect. Also mesh is not used correctly. Try this

for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,6);
    
    rc1 = [min(R), min(C)];
    rc2 = [max(R), max(C)];
    t = linspace(0,1,1000)';
    rcseg = (1-t)*rc1 + t*rc2;
    
    r = linspace(min(R), max(R), 1000);
    c = linspace(min(C), max(C), 1000);
    [Rg, Cg] = meshgrid(r, c);
    Fg = griddata(R, C, F, Rg, Cg);
    
    fseg = interp2(Rg,Cg,Fg,rcseg(:,1),rcseg(:,2)); %%this line is wrong, The number of input coordinate arrays does not equal the number of dimensions (NDIMS) of these arrays%%%%
    d = cumsum([0;sqrt(sum(diff(rcseg).^2,2))]);
    line_integeral = trapz(d,fseg)
end

[Ameer Hamza]

I am not sure about this issue. You can to see the specification of your dataset, that which column corresponds to which data. Also, I am not sure why two integrals should be equal.

[MS]

Thanks for the reply. It should be equal according to stokes theorem if we do the calculation correct.

In the data file(avg_mat). is there a way to find the line integeral like u.dr+v.dc for the same limits[Rg, Cg] .It involves two seprate integeration

integeration of u.dr with the limts[min(R), max(R)] + integeration of v.dc with the limts[min(C), max(C)]

can you help me to modify it. .

First column = X direction(R)

Second column =Y direction(C)

Third column = X direction_velocity(u)

Fourth column = Y direction_velocity(v)

[Ameer Hamza]

Sorry, I don't have much expertise in multivariate calculus.

Isn't this just two independent integrals. You can do it like this.

trapz(r,u) + trapz(c,v)

[Ameer Hamza]

For simple integrals, you don't need to use geiddata. You can do something like this

%%%%%%%%%%%%%%%%
for i=1:numel(avg_mat)
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    u = avg_mat{i}(:,3);
    v = avg_mat{i}(:,4);
    
    line_integeral(i)= trapz(R,u) + trapz(C,v)
    
%     fseg = interp2(Rg,Cg,Fg,rcseg(:,1),rcseg(:,2)); %%this line is wrong, The number of input coordinate arrays does not equal the number of dimensions (NDIMS) of these arrays%%%%
%     d = cumsum([0;sqrt(sum(diff(rcseg).^2,2))]);
%     line_integeral(i) = trapz(d,fseg)
end

Answer 3

[MS]

Thanks a bunch. Since my u and v values has got alot of 'NANs' the out put seen as 'NANs' . is there a way to avoid NANs by using any logical operator. Please include a way to add logical in the code.

[Ameer Hamza]

You can set NaN values to 0.

u(isnan(u)) = 0;
v(isnan(v)) = 0;

[MS]

Thank you very much.

[MS]

Hi ameer,

I want to do modfication(nested loop) in the limts[r,c].

I have twelve files(i = 12) to Process . But each file needs to use diffrent r and c values in the code. I need to mutiply r and c with cosine theta for every file.

where theta increases from 0deg to max deg and decreases to 0deg at the end with the increment (max/number of files). Please help.

eg,

for theta = [0:max:0] in degrees

first file limits

r = [min(R)*cos(0), min(C)*cos(0)];%add cos theta for every file

c = [max(R)*cos(0), max(C)*cos(0)];

middle file limts

r = [min(R)*cos(max), min(C)*cos(max)];%add cos theta for every file

c = [max(max)*cos(max), max(C)*cos(max)];

last file limts

r = [min(R)*cos(0), min(C)*cos(0)];%add cos theta for every file

c = [max(R)*cos(0), max(C)*cos(0)];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

results = zeros(numl(avg_mat),1);
for i=1:r and c
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,5);
max = 90;
for theta = [0:(max/numl(avg_mat)):max:-(max/numl(avg_mat)):0]
    r = linspace(min(R), max(R), 1000)*cos(theta);%need modification based the file
    c = linspace(min(C), max(C), 1000)*cos(theta);%need modification based on the file
    [Rg, Cg] = meshgrid(r, c);
    Fg = griddata(R, C, F, Rg, Cg);
    result(i) = trapz(c, trapz(r, Fg, 2));
end
end
writematrix(result, 'integral_val.txt')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[Ameer Hamza]

The code logic is getting quite complicated to understand. Only you can understand what is happening at each line. I suggest you add breakpoints in your code and run it line by line. Then you can see which line is not working as expected.

[MS]

hi ,

I made the code finally with all your help. i need to thank you for helping me with a free will. thank you so much.

regards,

MS

[Ameer Hamza]

I am glad that my comments were of some help in your project.

[MS]

your help is major for my project. now, i need help to plot (r(i),c(i)) as a rectangle for each file. can you help me to add a line.

 navg = numel(avg_mat);
r = cell(navg,1);
c = cell(navg,1);
maxi = 90;

%this code is designed to handle odd navg as well as even
t1 = linspace(0, maxi, ceil(navg/2)+1);
t2 = linspace(maxi, 0, navg - length(t1)+2);
tvals = [t1(2:end), t2(2:end)];

results = zeros(navg,1);
Fg = cell(navg,1);
for i = 1:navg
    writematrix(avg_mat{i}, ['avg_val_' num2str(i)]);
    R = avg_mat{i}(:,1);
    C = avg_mat{i}(:,2);
    F = avg_mat{i}(:,5);
    r{i} = linspace(min(R), max(R), 1000) * cosd(tvals(i));
    c{i} = linspace(min(C), max(C), 1000) * sind(tvals(i));
    [Rg, Cg] = meshgrid(r{i}, c{i});
    tg = griddata(R, C, F, Rg, Cg);
    
end

[Ameer Hamza]

See rectange(): https://www.mathworks.com/help/releases/R2020a/matlab/ref/rectangle.html

Answer 4

[MS]

Thank you for helping to find it.