Find Indefinite Integral from Array Values

Hi all,

I am trying to obtain indefinite integral of a function described by numerical values inside an array. To build such function i stardet trying to solve a simple integral of sin(x).

The problem i am facing is that i can't get the array containing the points representig the indefinite integral of the function (with constant of integration C=0). In other words vint_theor and vint_numeric must coincide.

a=0;
b=10;
N=10^6;
x=linspace(a,b,N);


v=sin(x);
vint_theor=-cos(x);
vq=(b-a)/N*cumtrapz(v);
vint_numeric=(vq-mean(vq));

plot(x,v,'DisplayName','sin(x)')
hold on
plot(x,vint_theor,'DisplayName','-cos(x)')
plot(x,vint_numeric,'DisplayName','numerical integration')
xlabel('n') 
ylabel('Amplitude')
hold off

legend

[Dyuman Joshi]

"In other words vint_theor and vint_numeric must coincide."

cumptraz() computes the approximate cumulative numerical integral (which is itself an approximation of analytical/symbolic integration).

So, no, they will not coincide.

If you need the exact result, you will have to use symoblic variables/expressions/functions.

Edit - John D'Errico has highlighted and resolved the issue.

[Mattia Deriu]

I'm trying to build a code that works on ARRAY values whose analytic expression is not known.

Answer 1

[John D'Errico]

Sorry. there is no need for them to be exactly the same.

Cumtrapz is a simple variation of a fixed stepsize Euler integration, as if you were solving a simple first order ODE, using Euler's method. It is an APPROXIMATION. It is not the exact solution to an integral. Take more terms, the error may decrease. But even then, Euler's method need not always be a viable solution.

Having said that, WHY ARE YOU DOING WHAT YOU DID???????? Why did you subtract off the mean?

a=0;
b=10;
N=10^6;
x=linspace(a,b,N);
v=sin(x);
vint_theor=-cos(x);

Your first error was here:

vq=(b-a)/N*cumtrapz(v);

What is the step size? linspace genrates a vector of length N. How many intervals are there?

(Answer: N-1)

Therefore, what is the step size?

(b-a)/(N-1)

Next, what is the purpose of this next line?

vq is an APPROXIMATE estimate of the integral. Well, it would be if you had divided by N-1.

vint_numeric=(vq-mean(vq));

But then why subtract of the mean? What would the mean of that expression be?

a=0;
b=10;
N=10^6;
x=linspace(a,b,N);
v=sin(x);
vint_theor = 1 - cos(x); % The actual integral of sin(x), from 0 to x
vint_numeric = cumtrapz(v)*(b-a)/(N-1);

plot(x,v,'DisplayName','sin(x)')
hold on
plot(x,vint_theor,'DisplayName','-cos(x)')
plot(x,vint_numeric,'DisplayName','cumtrapz')
xlabel('n') 
ylabel('Amplitude')
hold off
legend

As you can see, the yellow and red curves now overlay on top of each other.

[Mattia Deriu]

On your plot i see 1-cos(x) and not -cos(x) so it is not the output i was looking for. Many thanks for the N-1 hint anyway.

I'm using the analytic expression as benchmark so don't focus on that. This code shall work with complex oscillating functions having mean(f)=0 whose analytic expression is not known at all and i'm not intrested in it.

I have an array and i want to compute the best APPROXIMATION of its INdefinite integral plot without knowing its analytic expression considering integral constant C=0.

I later notice that my code shall work if the extension of the integrating interval is an interger multpile of 2*pi.

[Mattia Deriu]

a=0;
%% b=10;
b=4*pi;
N=10^6;
x=linspace(a,b,N);
v=sin(x);
vint_theor=-cos(x);
vq=(b-a)/N*cumtrapz(v);
vint_numeric=(vq-mean(vq));
plot(x,v,'DisplayName','sin(x)')
hold on
plot(x,vint_theor,'DisplayName','-cos(x)')
plot(x,vint_numeric,'DisplayName','numerical integration')
xlabel('n') 
ylabel('Amplitude')
hold off
legend

[Torsten]

I later notice that my code shall work if the extension of the integrating interval is an interger multpile of 2*pi.

If you wouldn't have subtracted mean(vq), your code had also worked if the integration interval were not a multiple of 2*pi. With 10^6 discretization points, you can approximate the integral of almost every function :-)

[Mattia Deriu]

If i wouldn't have subtracted mean(vq) my code wound end up with this (which again is not the INDEFINITE integral with constant C=0 and definitely not i was looking for). I'm looking toward getting the plot of -cos(x) which is not what i get from the code.

a=0;
b=10;
N=10^6;
x=linspace(a,b,N);
v=sin(x);
vint_theor=-cos(x);
vq=(b-a)/N*cumtrapz(v);
vint_numeric=(vq);
plot(x,v,'DisplayName','sin(x)')
hold on
plot(x,vint_theor,'DisplayName','-cos(x)')
plot(x,vint_numeric,'DisplayName','numerical integration')
xlabel('n') 
ylabel('Amplitude')
hold off
legend

[Torsten]

You get -cos(x) + cos(0), as expected which is integral_{t=a}^{t=b} sin(t) dt.

I don't know what you mean with "constant C = 0".