Calculation of Arc length of grade 2 and grade 4 polynomial function

Hello,

I am a new user of Matlab with no experience, but I try to analyse some of the data of a research project using it. The task is to calculate the arc length of a detected surface over the time. The surface at the beginning has the shape of a squared function, after a certain period of time a deformation starts and the surface is following a polynomial function 4th grade. I tried to solve the problem with Excel but I can not integrate the square root of a 4th-grade function. Some research on the internet showed that this would be possible with Matlab, but unfortunately, I can not find a solution by myself.

I am very thankful for any idea! Thank you in advance.

Answer 1

[Star Strider]

I would use the trapz or cumtrapz function for a vector of data.

For example (from the Wikipedia arc length article):

t = linspace(-1, 1);                                % Create Data
dt = t(2)-t(1);
x = t.^3;
y = t.^5;

dx = gradient(x,dt);                                % Derivatives
dy = gradient(y,dt);
ds = hypot(dx, dy);                                 % Arc Segment Vector
cum_arclen = cumtrapz(t, ds);                       % Cumulative Integration
tot_arclen = cum_arclen(end);                       % Total Arc Length (Equivalent To Using ‘trapz’)

figure(1)
plot(x,y,   x,cum_arclen)
grid
axis equal
legend('Original Data', 'Arc Length', 'Location','NW')

EDIT — Added ‘example’ code.

[doehr001]

Thank you very much for the answer. Would it be possible to explain it a little bit more in detail? I tried it now, but, unfortunately, it is not working. The table of data has the format: the first column is the time the second column is the outermost point of the profile (xmin) all next columns are along the x-axis from -4to4 (572 columns) and each row contains the position of the surface on the y-axis.

I would need a curve fitting and based on that the arc length of the 572 x-columns (y(t) at each x).

Hope this doesn't confuse more, than it helps.

[doehr001]

Thank you for the code. I tried it by using my own data. Unfortunately, I get an error when I run it.

I deleted the create data code and use this:

dt = t(2)-t(1);                                     % Derivatives
dx = gradient(x,dt);                                
dy = gradient(y,dt);

ds = hypot(dx, dy);                                 % Arc Segment Vector

cum_arclen = cumtrapz(t, ds);                       % Cumulative Integration

tot_arclen = cum_arclen(end);                       % Total Arc Length (Equivalent To Using ‘trapz’)

figure(1)
plot(x,y,   x,cum_arclen)
grid
axis equal
legend('Original Data', 'Arc Length', 'Location','NW')

    if true
      % code
    end

The errors are:

>> createFit Error using cumtrapz (line 71) LENGTH(X) must equal the length of Y in dim 2.

Error in createFit (line 8) cum_arclen = cumtrapz(t, ds); % Cumulative Integration

>>

The length of x and y is equal (please see attached picture).

[Star Strider]

Attach your data, preferably as a .mat file (use the save function), and describe what the variables are, and what you want to do. Also, post your code. (My code works with the data I created to test it.) If your data are both functions of time, also provide the time vector.

I cannot troubleshoot my code with your data without having your data to work with.

[doehr001]

Sorry. Attached you will find the file. I use the following code:

dt = t(2)-t(1);
dx = gradient(x,dt);                                % Derivatives
dy = gradient(y,dt);
ds = hypot(dx, dy);                                 % Arc Segment Vector
cum_arclen = cumtrapz(t, ds);                       % Cumulative Integration
tot_arclen = cum_arclen(end);                       % Total Arc Length (Equivalent To Using ?trapz?)
figure(1)
plot(x,y,   x,cum_arclen)
grid
axis equal
legend('Original Data', 'Arc Length', 'Location','NW')

It is the code you posted before just deleted the creation data bit.

t is the time from 0 to 32ms,

x is the profile of the surface from -4 to 4 mm,

y is the position of the surface across the profile

I want to calculate the arc length of the surface at each time point.

Thank you very much for your help!

[Star Strider]

I can’t load the file. It’s not a binary file, and when I attempt to load it as an ASCII file, I get this error:

Error using load
Number of columns on line 3 of ASCII file ArcLength
StarStrider.mat must be the same as previous lines.

Error in ...
D = load('ArcLength StarStrider.mat','-ascii');

Please save it as a binary file.

EDIT — I cannot help unless I have your data, and as I mentioned, i cannot load or read the file you attached. Please save your data as a binary .mat file, as I requested. See the documentation on the save function for details.

[doehr001]

I can only save a .mat by adding .mat manually. Mathlab safe it automatically as a .m file. Is this working for you?

I tried it again with the code and attached another .mat file. Please try this one. I hope it works. Thank you very much for your patience!

[Star Strider]

That seems to have worked!

What is the ‘yAll’ matrix? When I plot it, it looks like a (140x572) parabolic surface. What do you want to do with it?

What axis is time? What is the sampling frequency or sampling interval? Please provide the time vector if the sampling interval is not constant.

[doehr001]

Great! Thank you for looking at it.

the yall matrix is the complete recorded matrix for one sample, each row is one time point. The time axis is in t and starts from 0 to 139, increment 1. The x axis start from -4 to 3.9940 in 0.0140 steps.

I copied the data in the particular table.

[Star Strider]

Apologise for the delay.

This works for me:

D = load('ArcLengthStarStrider.mat');
yAll = D.yAll;

t = 0:size(yAll,1)-1;
xax = linspace(-4, 3.9940, size(yAll,2));
dx = gradient(xax,1);                               % X-Derivative

for k1 = 1:length(t)
    dy = gradient(yAll(k1,:),dx(1));                % Y-Derivatives
    ds = hypot(dx, dy);                             % Arc Length Segment
    arclen(k1) = trapz(xax, ds);                    % Total Arc Length
end

figure(1)
meshc(yAll)
xt = get(gca,'XTick');
xta = [1 xt size(yAll,2)];
set(gca, 'XTick',xta, 'XTickLabel',linspace(-4, 3.9940, length(xta)))
xlabel('x')
ylabel('Time')

figure(2)
plot(t, arclen)
grid
xlabel('Time')
ylabel('Arc Length')

The ‘arclen’ variable is the vector of arclengths for each value of ‘t’, as plotted in figure(2). The values seem to me to be appropriate, given the magnitudes of ‘x’ and ‘y’ in figure(1).

[doehr001]

Brilliant thank you so much! This is exactly what I needed! Many times thank you.

[Star Strider]

As always, my pleasure!