I am trying to implement Romberg's integration process in MAT LAB (i.e. recursive-trapezoidal to simpson's 1/3 rule(or Richardsons's extrapolation) to finally Boole's rule. After computing my answers for the trapezoidal rule, I then call on the results from my trapezoidal code to then apply Romberg intergration for Simpson's 1/3 rule. And then call on the results from the Rom berg simpson code to finally apply Romberg integration using Boole's rule.
My question is simply around the inability to call the results effectively. My code for Romberg Simpson's 1/3 is able to call on the trapezoidal results correctly, but for some reason my Romberg boole code keeps calling the trapzoidal data instead of the newly extrapolated Romberg Simpson's 1/3 data and applies the equation. And I dont understand why?
n=1; %first value of area (i.e. trapezoidal rule) result(n)=((b-a)./2*(func(b)-func(a))) %start error error_result=1; %while loop that creates error bounds. Will cycle until error is larger %than toler while error_result(n)>toler; it=2.^(n-1); tnm=it; h=(b-a)./tnm; %first x values is from startig from a +0.5 h interval due to split x=a+0.5.*h; %first sum=0 sum=0.0; for j=1:it %sum is sum of sum=sum+feval(func,x); %finding next x value x=x+h; end %Answer/Result of integration at given interval split result(n+1)=0.5.*(result(n)+h.*sum); %error analysis for each interval split error_result(n+1)=abs((result(n+1)-result(n))); %counter n=n+1; end
function[result]=Rombergsimp(func,a,b,toler) %starting counter n=1; %error starter error_result=1; while error_result(n)>toler; %calling result values from trapezoidal code result = trapz(func,a,b,toler); %defining new result values by applying simpsons 1/3 romberg integration. result_romsim(n+1)=result(n+1)+1/3.*(result(n+1)-result(n)); %error analysis of each step error_result(n+1)=abs(((result_romsim(n+1)-result_romsim(n)))); %counter+1 n=n+1; %printing results end
function[result]=Rombergboole(func,a,b,toler) %counter n=1; %error starter error_result=1; %while lopp defining error analysis while error_result(n)>toler; %calling result values from Rombergsimp m-function result = Rombergsimp(func,a,b,toler) %Applying Boole's rule to the already extrapolated simpson data result_romboole(n+1)=result(n)+16/15.*(result(n+1)-result(n)); %displaying results disp(result_romboole); %error check between each new Romber_boole value. error_result(n+1)=abs(((result_romboole(n+1)-result_romboole(n)))); %counter+1 n=n+1; end
I am sorry for the lengthiness of the code and question but I didn't know how else to ask it.
I don’t completely understand how you’re calling your routines, but it seems if you call ‘Rombergboole’ alone, it calls all the others, first ‘Rombergsimp’ that itself calls ‘trapz’ and returns the eventual result.
I don’t understand what ‘func’ is. The trapz function takes a vector (or array) argument — not a function argument — and produces a scalar (or vector) result. There would seem to me to be nothing left to integrate. In other places, it seems that you are calling ‘func’ as a function handle. If you want to integrate a function, the integral (or in earlier versions, quad) function would be appropriate.