Muhammed Tuğrul
Muhammed Tuğrul

Reputation: 15

NUMERİCAL METHOD Trapezoidal Rule Fails When N=3 WHY?

clc
clear all
close all
f=@(x) (exp(3*x)+sin(x));
x0=pi/2;
xn=pi;
%x0=input('lower= ');
%xn=input('upper= ');
N=input('N= ');
h=((xn-x0)/N);
area =0;
while (x0<xn)
   area=area+(h/2)*(f(x0)+f(x0+h));
   x0=x0+h;
end
fprintf('area=%f',area);

Everything is okey in this code without entering the N value as 3. There is a fail why?

Upvotes: 0

Views: 37

Answers (1)

Tommaso Belluzzo
Tommaso Belluzzo

Reputation: 23675

This works:

% ...
N = input('N= ');
h = (xn-x0) / N;
area = 0;

for i=1:N
    area = area + ((h/2) * (f(x0) + f(x0+h)));
    x0 = x0 + h;
end

fprintf('area=%f',area);

and now it's time to explain you why.

Since you are dividing the difference between xn and x0 by N, the number of loops needed in order to obtain x0 = xn are equal to N because at the end of every loop you do x0 = x0 + h. If N = 3, you automatically know that after the third loop you will need to stop since x0 = xn.

Using a while and comparing x0 and xn is dangerous, since they are both double precision floating-point numbers and they may loss precision during the iterations. There is the possibility that they will not match even if they are very very very very very very very ... very close to each other.

Using your approach, due to loss of precision, after the third loop you will have x0 and xn almost identical... well, identical if you round them up to like 10 digits or so, but not exactly identical on a higher precision after the decimal dot... so your while condition will fail and the loop will continue.

If you want to perform a good comparison between to floating-point numbers keeping a high precision profile you need to check for:

abs(a - b) > eps(max(abs(a), abs(b)))

like (this while can do the job within your code):

while (abs(xn - x0) > eps(max(abs(xn), abs(x0))))
    %...
end

In other words, you have to define a tolerance level.

Upvotes: 1

Related Questions