Adaptive step size trapezoidal integration

[gfd43tg]

https://www.physicsforums.com/attachments/72086

https://www.physicsforums.com/attachment.php?attachmentid=72087&d=1407804589

Hello

Here is the code for the adaptive stepsize function

function I = arttrap(fh,a,b,tol,fa,fb)
if nargin == 4
    fa = fh(a); fb = fh(b);
end
m = (a+b)/2;
fm = fh(m);
h = b-a;
% Compute I based on 1 subdivision
I1 = 0.5*(fa+fb)*h;
% Compute I based on 2 subdivisions
I2 = (fa + 2*fm + fb)*h/4;
% Compare both estimates
if abs(I1-I2)<tol
    I = I2;
else
    ILeft = artrap(fh,a,m,tol/2,fa,fm);
    IRight = artrap(fh,m,b,tol/2,fm,fb);
    I = ILeft + IRight;
end
disp(['a = ' num2str(a) ', b=' num2str(b)]);

I am wondering if I did it correctly, or better yet if just by looking at the graph I would be able to tell the values of $x$ which are evaluated. This can't be plotted to test, so it's clear this was a contrived graph to test the student's understanding of how adaptive stepsize integration works. I am just wondering what would indicate the places that would be checked.

From what I see, 2 and 10 would be looked at. Then also 6 and 8, since there is a sharp turn at those locations, and a step could be placed there. But my calculations are saying that it's also evaluated at 7 and 9.

 

[AlephZero]

Your two links don't work.

 

[gfd43tg]

Okay my bad, there they are

 

[AlephZero]

From what I see, 2 and 10 would be looked at. Then also 6 and 8, since there is a sharp turn at those locations, and a step could be placed there. But my calculations are saying that it's also evaluated at 7 and 9.

The computer doesn't magically "know" where the kinks in the graph are and put the steps and the "best" places. It just does what the program tells it to do.

I think you are right about evaluating at 7 and 9. The interval 6 to 10 doesn't converge so it is split into 6 ro 8, and 8 to 10. The computer doesn't "know" what the graph looks like, so it has to evaluate at the mid point of those intervals. Both intervals converge so they are not split any further.

The marks on the graph, and the answer on the second sheet, forgot about checking the interval 2 to 6.

BTW this is a poor algorithm for doing integration (but it's a nice exercise related to trees and recursive functions). See why you get the wrong answer, if you use it to evaluate $\int_0^{2\pi} \sin^2\!\! x\, dx$.

 

[.Scott]

To be clear, when AlephOne said "The marks on the graph, and the answer on the second sheet, forgot about checking the interval 2 to 6.", he's saying you've missed a value (which you have).

 

[gfd43tg]

So that means my answers are right? Even the 2nd question?

 

[.Scott]

So that means my answers are right? Even the 2nd question?

You're list (2,10,6,7,8,9) is missing a number.
You're estimation of what will be printed is also wrong.

Here's what will happen:

artrap(2,10) called.
--fh evaluated at 2, 10, and then 6
--artrap(2,6) called.
----fh evaluated at <you fill this in>
----display 2, 6
--artrap(6,10) called.
----fh evaluated at <you fill this in>
----artrap(6,8) called.
...

Notice that the first display is "a=2, b=6", not "a=2, b=10" which is what you have.

So in the attachments, both H.1 and H.2 have mistakes.

 

[gfd43tg]

So it evaluates at 2,10,6,4,8,7,9.

I am confused about where it has its end points. This code gets confusing because when it splits up the step, it doesn't meet tolerance requirements, so it splits itself again, then I don't understand how to get it to go back.

Does it display
1. a=2 b=6
2. a=6 b=8
3. a=8 b=10
?
If so, is that just coincidentally the places on the graph with the sharp turnsartrap(2,10) called.
--fh evaluated at 2, 10, and then 6
--artrap(2,6) called.
----fh evaluated at <2,4,6>
----display 2, 6
--artrap(6,10) called.
----fh evaluated at <6,8,10>
----artrap(6,8) called.
---- fh evaluated at <6,7,8>
----artrap(8,10) called
---- fh evaluated at <8,9,10>

 

[.Scott]

So it evaluates at 2,10,6,4,8,7,9.

Yes.

I am confused about where it has its end points. This code gets confusing because when it splits up the step, it doesn't meet tolerance requirements, so it splits itself again, then I don't understand how to get it to go back.

Does it display
1. a=2 b=6
2. a=6 b=8
3. a=8 b=10
?
If so, is that just coincidentally the places on the graph with the sharp turns

Those are the first three lines displayed. But there are two more.
The first call to artrap is with a=2, b=10. Since the last line in artrap reports its input, you should know that the last thing that will be reported is "a=2, b=10". You might want to work backwards through the list.

artrap(2,10) called.
--fh evaluated at 2, 10, and then 6
--artrap(2,6) called.
----fh evaluated at <2,4,6>
----display 2, 6
--artrap(6,10) called.
----fh evaluated at <6,8,10>
----artrap(6,8) called.
---- fh evaluated at <6,7,8>
----artrap(8,10) called
---- fh evaluated at <8,9,10>

Not quite.
artrap(2,10) called.
--fh evaluated at 2, 10, and then 6
--artrap(2,6) called.
----fh evaluated at 4 (and only 4)
----display 2, 6
--artrap(6,10) called.
----fh evaluated at 8 (and only 8)
----artrap(6,8) called.
------fh evaluated at 7 (and only 7)
------display <you fill this in>
----artrap(8,10) called.
------fh evaluated at 9 (and only 9)
------display <you fill this in>
----display <you fill this in>
--display <you fill this in>


Notice that I am using two dashes for each level in the stack. It's important to keep track of this.