MATLAB-output values matrix

[Bostonpancake0]

Hi,

I'm trying to write a code for the recursive trapezoidal rule. The code must have the bounds of integration and tolerance. The code is meant to stop when the tolerance between the next and previous value is below 10^-4, given n number of times to sub split the intervals.

My progress on the code is below;

function [s1] = trapzd(func,a,b,toler)
%
% function s1 = trapzd(func,a,b,s,n)
%
% This routine uses the trapezoidal
%rule to evaluate
% the integral between a and b of the
%integrand defined
% by the function func. n is the number
%of times the interval is subdivided.
%
% The new approximation (s1) is
%adjusted from the previous
%approximation
% (s) from a prior call.
%

s=0

for n=1:10

it=2.^(n-2);
tnm=it;
h=(b-a)./tnm;
x=a+0.5.*h;
sum=0.0;


for j=1:it
sum=sum+feval(func,x);
x=x+h;
end
s=0.5.*(s+h.*sum);

s1 = s

end


end

The code seems to be working fine but I'm wondering if there is anyway to put the out put values for s1 into a matrix immediately after one loop. What I am hoping to achieve is, is to put a while loop around all other for loops stating that the code will stop until s1(n+1)-s1(n) is less than 10^-4. Its probably an easy answer to this question but I just cannot remember how to do it.

Also if there is any other tips that can be given will be much appreciated.


Thanks

 

[kreil]

I'm wondering if there is anyway to put the out put values for s1 into a matrix immediately after one loop.

Can you be more specific? Are you saying you think your code can be written with 1 loop rather than two, or that you don't like the way s1 is defined?

What I am hoping to achieve is, is to put a while loop around all other for loops stating that the code will stop until s1(n+1)-s1(n) is less than 10^-4. Its probably an easy answer to this question but I just cannot remember how to do it.

You'll have to store the s1(n) value to later compare it to the s1(n+1) value. Something like this should work, but I haven't tested it much:

s=0;
ds=1;

while ds > toler
    for n=1:10
        sN = s;
        it=2.^(n-2);
        tnm=it;
        h=(b-a)./tnm;
        x=a+0.5.*h;
        sum=0.0;
        
        for j=1:it
            sum=sum+feval(func,x);
            x=x+h;
        end
        s=0.5.*(s+h.*sum);
        
        s1 = s;
        ds = abs(s1-sN);
    end
end

You might also wish to add some maximum number of iterations to an implementation like this, to ensure it doesn't get stuck trying to reach an unachievable accuracy. Just add a ticker to count the iterations and error out when it goes above some threshold.