Find limits of known integral

11 visualizzazioni (ultimi 30 giorni)
Adam
Adam il 30 Lug 2014
Commentato: John D'Errico il 2 Ago 2014
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
  4 Commenti
Adam
Adam il 30 Lug 2014
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
John D'Errico il 30 Lug 2014
Modificato: John D'Errico il 30 Lug 2014
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.)

Accedi per commentare.

Risposte (3)

Adam
Adam il 30 Lug 2014
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
  2 Commenti
John D'Errico
John D'Errico il 30 Lug 2014
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
John D'Errico il 30 Lug 2014
Why would you want to integrate a spline using the trapezoidal method? That is just silly.

Accedi per commentare.


Adam
Adam il 30 Lug 2014
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
  2 Commenti
John D'Errico
John D'Errico il 30 Lug 2014
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
John D'Errico il 30 Lug 2014
Modificato: John D'Errico il 30 Lug 2014
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.

Accedi per commentare.


John D'Errico
John D'Errico il 30 Lug 2014
Modificato: John D'Errico il 30 Lug 2014
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.)
  3 Commenti
Adam
Adam il 31 Lug 2014
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
John D'Errico il 2 Ago 2014
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.

Accedi per commentare.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by