Integration of a derivative(arc length formula)

while (sz(1)-1 > 0)
                 x = linspace(points(z-1),points(z));
                 y = X^2+2*z; 
                 dy = diff(y); 
                 func= sqrt(1+dy.^2); %%arc length formula
                 I = I + integral(func,points(z-1),points(z)); %Error here
                 sz(1)=sz(1)-1;
                 z = z+1;
end

Error is

First input argument must be a function handle.

Can someone please help?

Thanks

Answer 1

[John D'Errico]

I have absolutely no idea what it is you are trying to do, because your code makes no sense.

So instead, I'll just show you how to compute arclength, using several methods. I'll pick a simple function - y=x^2. So what is the arclength of the curve (x,x^2), where x varies from 0 to 2?

First, we can compute that result in a symbolic form, and hope that the toolbox can handle it. The sqrt in there might be a problem.

syms x
y(x) = x^2;

int(sqrt(1 + diff(y).^2),x)
ans(x) =
asinh(2*x)/4 + x*(x^2 + 1/4)^(1/2)

int(sqrt(1 + diff(y).^2),x,[0,2])
ans =
log(17^(1/2) + 4)/4 + 17^(1/2)

So over the desired interval, we got a symbolic solution.

vpa(ans)
ans =
4.646783762432935873382616

Next, we can solve this using other schemes. Next, integral:

dy = @(x) 2*x;
integral(@(x) sqrt(1 + dy(x).^2),0,2)
ans =
          4.64678376243294

Integral works quite well. Of course there I computed dy/dx myself.

Next, we can use trapz.

x = linspace(0,2,100);
yfun = @(x) x.^2;

sum(sqrt(diff(x).^2 + diff(yfun(x)).^2))
ans =
          4.64675076773921

Trapz did quite well, considering that we only used 100 points along the curve.

Next, you could download my arclength utility from the file exchange. (I'll add a link in a minute.) Arclength actually passes a spline though those points, then computed the arclength along the spline interpolant itself. It forms an integral as needed internally. As you should expect the arclength computed using the integral of a spline interpolant should be more ccurate than what we get from trapz on the same set of points. But it is not as accurate as what we get from the sym solution or from integral.

arclength(x,yfun(x),'spline')
ans =
          4.64678376533212

You can find arclength here for free download:

https://www.mathworks.com/matlabcentral/fileexchange/34871-arclength

Next, we can use an ODE solver to solve the problem. A nice thing about the ODE solver is we also get the arclength at intermediate points along the curve. We would have gotten something similarly useful before, had we used cumtrapz instead of trapz.

dy = @(x) 2*x;
tfun = @(t,Y) sqrt(1 + dy(t).^2);
[tout,lenout] = ode45(tfun,[0,2],0);
lenout(end)
ans =
          4.64678373566526
          
plot(tout,lenout)
grid on

So ODE45 did quite well too.

[Adham Elkhouly]

y = @(x) D(z-1)./(6.*(points(z)-points(z-1))).*(points(z)-x).^3+D(z)./(6.*(points(z)-points(z-1))).*(x - points(z-1)).^3+(points(z-1,2)./(points(z)-points(z-1))-(D(z-1).*(points(z)-points(z-1)))./6).*(points(z)-x)+(points(z,2)./(points(z)-points(z-1))-(D(z).*(points(z)-points(z-1)))./6).*(x-points(z-1));
 %please ignore what y equal cosider it 2x+2+z or something 
 I = I+  int(sqr(1+diff(y).^2),x, [point(z-1),points(z)]);

Thank you so much for replying

I get the following error: Undefined function 'diff' for input arguments of type 'function_handle'.

(I am trying to construct cubic splines and calculate the total length)

[Adham Elkhouly]

I don't know what the expression of the derivative will be because it changes every iteration. I just want to be able to compute the function's derivative and use it in the integral, which produces an integer.

[Adham Elkhouly]

I solved the issue - Thank you

Just one question what does vpa produce? I tried to convert it to string with num2str but it says that num2str takes only numeric values

[John D'Errico]

I'm tying to figure out what function this represents in y as you wrote it.

[Adham Elkhouly]

This function represented the cubic spline interpolation formula

[haha mark]

Can I ask you for help on this issue：https://www.mathworks.com/matlabcentral/answers/525051-how-to-integrate-discrete-values-over-a-known-x-y-coordinate-image?s_tid=prof_contriblnk