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
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?
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.)
