Find limits of known integral
1 回表示 (過去 30 日間)
古いコメントを表示
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 件のコメント
John D'Errico
2014 年 7 月 30 日
編集済み: John D'Errico
2014 年 7 月 30 日
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.)
回答 (3 件)
Adam
2014 年 7 月 30 日
2 件のコメント
John D'Errico
2014 年 7 月 30 日
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
2014 年 7 月 30 日
Why would you want to integrate a spline using the trapezoidal method? That is just silly.
Adam
2014 年 7 月 30 日
2 件のコメント
John D'Errico
2014 年 7 月 30 日
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
2014 年 7 月 30 日
編集済み: John D'Errico
2014 年 7 月 30 日
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.
John D'Errico
2014 年 7 月 30 日
編集済み: John D'Errico
2014 年 7 月 30 日
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 件のコメント
John D'Errico
2014 年 8 月 2 日
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.
参考
カテゴリ
Help Center および File Exchange で Boundary Conditions についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!