Find limits of known integral

Dear all,

I would like to find the limits ( -a and +a) of a known integral.That means the integral is known (e.g. I=45) but the limits -a and +a are unkown. How can I find these limits?

This my example:

.......

.......

.......

xx=linspace(x(1);

x(end),8000);

pp=spline(x,y);

yy=ppval(pp,xx);

%%integral over [-a,+a] is known

I=quad(@ppval,-a,+a,[],[],pp)

.....

.....

....

Your help is highly welcome and aprreciated

Cheers

[Star Strider]

Without knowing what the function is, my guess would be ±NaN.

Seriously though, is it defined over ±Inf or some other interval? What values does it take? Is it symmetrical, even, odd, or what?

What your question lacks in detail, it makes up for in ambiguity.

[John D'Errico]

There is no need to know if the function is symmetrical, etc., merely that it is a known function, and that it has some integral we can evaluate.

[Adam]

Hi all,

the function is unkown. I have only data (Point Spread Function: PSF is symmetric)

x=PSF(:,1);

y=PSF(:,2);

data points are interolated using spline. Integrating over all the points is simple and it is the area A under the PSF. But assuming now that the integral between -a and +a is known and it is , for example, 0.35*A. What is the value of a in this case?

Chhers

[John D'Errico]

Again, I fail to see the problem. Do exactly as I did, but now you need to integrate a spline instead of some known function. While that is most efficiently and accurately done analytically, it is still trivial to do using a numerical tool such as quad. WTP? (You have shown in your own example that you know how one can integrate a spline, which while not the best way, is still serviceable.)

Answer 1

[Adam]

Dear John,

Thank you very much for your help. I greatly appreciate it! My PSF is not pur Gaussian and therefore I interpolated the data points usimg Spline instead of fitting a Gaussian function.

My question now is how to integrate Spline using Trapezoidal method. For a known function, eg. a Gaussian function, the integral can be calculated by Trapezoidal method but for Spline it does not work! To integrate a function using trapezoidal method I am calling the following function:

function I = Trapez(f,a,b,n)

x = a;

h = (b - a)/n;

s = feval(f,a);

for i = 1 : n-1

x = x + h;

s = s + 2*feval(f,x);

end

s = s + feval(f,b);

I = (b - a) * s/(2*n);

How can I integrate the Spline using Trap. Method?

Cheers

[John D'Errico]

First of all, PLEASE DON'T add a new answer for every comment you make! Use comments when you make a comment, or when you have a question about a true answer! This is NOT an answer that you have added.

[John D'Errico]

Why would you want to integrate a spline using the trapezoidal method? That is just silly.

Answer 2

[Adam]

Hi John,

I think the problem will remain the same beacuse in your previous posting you were dealing with know function. You can integrate and find the desired limits of the integral. Spline and interpl are unkonw functions the probelm is when I call the function:

I = Trapez(f,a,b,n) the argument f is an unknown function and thus can not be integrated. maybe I am worng!

cheers

[John D'Errico]

Again, DON'T ADD NEW ANSWERS when you want to make a comment! There is a comment button for that. USE IT.

[John D'Errico]

And a spline IS a known function. It is trivial to integrate a spline, because it is a piecewise polynomial. But if that is too difficult, as you showed yourself, integrating a spline using a numerical integration is also trivial.

Again, a spline IS a known function. See my edit to the answer I posed above.

Answer 3

[John D'Errico]

So you have

1. Some known function.

2. A method to integrate the function. If the function is a spline, then don't use quad to integrate it! Splines can be integrated analytically, for a more accurate result that takes less time.

3. A given value for the integral.

The answer is clear. Use a rootfinder, thus fzero. So create a function handle that calls the integrator, then subtracts off the given value of the integral. Fzero will find the value for the limits that yields zero for that difference.

Whats the problem?

I'll show an example...

Given f(x) = cos(x), and assume we are given the integral over [-a,a] is 1. It is the stuff of a basic calculus homework assignment, made trivial with the symbolic toolbox.

syms a x

int(cos(x),-a,a)
ans =
2*sin(a)

So we know that

2*sin(a) == 1

Solving...

solve(2*sin(a) == 1)
ans =
   pi/6
(5*pi)/6

So a == pi/6.

If you wanted to do it numerically...

fun = @(a) integral(@(x) cos(x),-a,a) - 1;

See that fun is a function we can evaluate, for any value of a.

fun(.5)
ans =
    -0.041149

fzero(fun,[0,2])
ans =
       0.5236

Did we get the correct answer?

pi/6
ans =
       0.5236

So next, since it seems I did not make myself clear enough about a spline being a function, and a KNOWN function, lets solve it using a spline.

First, I'll create a spline approximation to cos(x), over the interval [-pi,pi].

x = linspace(-pi,pi,51);
pp = spline(x,cos(x));

I'll use my own SLMPAR to do the integration. Find it as part of my SLM toolbox on the file exchange.

fun = @(a) slmpar(pp,'integral',[-a,a]) - 1;
fzero(fun,[0,2])
ans =
     0.5236

Since I could have done this with ANY set of points where I fit a spline to them, this shows that a spline is easy to work with here.

Just to beat a dead horse, we can use other tools too. Here is an example of the use of chebfun on my cos example, to show that it can generate what is essentially an exact result. (CHEBFUN is another toolbox found on the FEX.)

fun = chebfun(@cos,[-pi,pi]);
funq = @(a) integral(fun,-a,a) - 1;

format long g
A = fzero(funq,[0,2])
A =
         0.523598775598299

pi/6
ans =
         0.523598775598299

So that worked nicely indeed. Of course, it was not a spline. How about if we find the 95% central region of a normal distribution?

x = linspace(-10,10,101)';
y = normpdf(x);
pp = spline(x,y);
fun = chebfun(@(x) ppval(pp,x),[-10,10]);
funq = @(a) integral(fun,-a,a) - 0.95;
A = fzero(funq,[0,10])
A =
           1.95995957360485

Of course, given that it was based on a spline approximation to the points, I don't expect perfection.

norminv(0.975)
ans =
          1.95996398454005

(Note: I imagine I could have done this more elegantly using chebfun, but I am but a novice with that tool.)

[John D'Errico]

Here. I deleted my first answer, and reposted the same text. But please don't keep adding new answers everytime you make a comment. USE THE COMMENTS instead!

[Adam]

Hi john,

I am sorry about it. I think in your first answer you mentionned another method to find the limits of known integral.

Anyway, my problem is still not solved. I fellow your answer. But chebfun and integrate functions do not work. Therefore, I integrate my function using quad:

F=quad(@ppval,x(1),x(end),[],[],pp)

After that I wanted to find the limits of the integral (-a and +a) but still not working!! any additional hepl is highly appreciated.

I want to add an additional information. Why I want to find the limits of the integral?

1- Firstly I calculated the integral, let say Int1, over the whole interval [x(1), x(end)].

2- In the second step I have to take only 60 % (Int2=0.6*int1) of the central area under the PSF and then find the corresponding limits of the this Int2.

is it possible to do it directly in one step in the case of my PSF which not pure Gaussian as you mentioned in the example of normal distribution?

best regards...

cheers

[John D'Errico]

Wanting to do it in one line is silly. Is ONE line truly better than a few short ones? It won't be more efficient, faster, less memory intensive, etc. MATLAB does not charge you by the character. Writing your code in smaller, readable chunks is better code. It is more readable, more easily debugged.

As far as not being able to use CHEBFUN or SLM, you need to download them FIRST. Beyond that, I can't know what your problems were if you give me no information. I showed you at least three ways to solve a general problem. Read what I said, then if you have a specific problem, ask for clarification, SHOWING ME WHAT WAS YOUR PROBLEM. If you got an error message, paste it in, in its entirety.