how to program integration method with a variable step ?

Hello,

Please I nedd your help, I spent a lot of time searching how to make a variable step on an integration method like trapeze or simpson, I tried to variate h the step on the area of the signal. I use this code , but I'm confused !

The purpose is to variate h to not have the same h for N subdivisions. h=(b-a)/N ( fixed sample)

Iexa = (2*(sqrt(9)-sqrt(5)));
alpha = 0;
beta = 4;
%t=linspace(0,4);
%l=rand(N,1)
h =  alpha + ((beta-alpha)./N).*rand(N,1);
% for j=1:N
%     h = (beta-alpha).*rand(0,1);
% %   h=(t(j+1)-t(j));
% %   t=t+j;
%  end 
x = [alpha:h:beta];
f = inline('1./sqrt(5+x)','x'); 
Isim =zeros(N,1);
for i=1:N
    Isim= Isim+h*(1/6*f(x(i))+2/3*f((x(i)+x(i+1))/2)+1/6*f(x(i+1)));
end
Isim 
Esim= abs(Iexa - Isim)

The obtained values of Isim aren't raisonable. Do you have an Idea to make a non uniform sampling on integral method ? Thank you

Answer 1

[John D'Errico]

Note that trapz allows you to provide any spacing your heart desires, because you can give it a vector of values for x and y. So why are you writing your own code here? Never write code to solve a problem when professionally written code already exists.

If you feel you need a more accurate result, could you write a Simpson's rule integral that uses a variable step? Well, yes, in theory, you could do something like that. In practice, it is just a bad idea, because there are better tools out there. What tools?

1. Build a cubic interpolating spline. (A pchip interpolant might be a better choice some of the time.)
2. Integrate the interpolant from step 1. That would be done most easily using fnint (as long as you have the curve fitting toolbox.)

For example...

x = [0,logspace(-1,0,5)]
x =
            0          0.1      0.17783      0.31623      0.56234            1
y = exp(x);

Can we compute the integral of this function from 0 to 1 from the data? Analytically, the integral is exp(1)-1 == 1.71828...

format long g
exp(1) - 1
ans =
  1.71828182845905

Trapz is trivial to use:

trapz(x,y)
ans =
  1.73595871210986

Better is to use a spline integral, as I suggested.

spl = spline(x,y);
fnval(fnint(spl),1)
ans =
          1.71853051923856

Pchip would not have been as accurate for this particular function. But for SOME problems, pchip would have been a better choice.

[Sarah CHOUCHENE]

Hello John D'Errico,

Thank you very mush for your quick answer, I'll search more about Pchip and spline integrals, I'm downloading the curve fitting toolbox.

However, I make my own program to variate the step h, because my main purpose is to comparate between uniform subdivisions and the non-uniform one. So to do that, I have to write my own 'x' right ?

I understand that in theory it is possible but numerically it is not to variate the step of integral in numerical method such simpson or trapeze..

Thanks for your help.

Regards,