numerical integration

Recently, when I use Gauss-Legendre Formula to handle it. But it is not correct. So I want someone to help me to solve this problem, thanks for your attention!

This equation is: 1/(1-18/9*x+8/9*exp(3*(x-1))), and the integrating range is [1,H], and 0<=H<=1.

Generally speaking, for any H, this result is positive, but for my code, there are some negetive number, and this puzzle me for a very long time, and did not fint the answer.

so, now i give my code. And i do not know how to handle it.

function s=traprl(f,a,b)
% [a,b] is the integrate region, b always equals unity, and f is the function 
h=(b-a)/10000;
s=0;
for k=1:10000-1
    x=a+h*k;
    s=s+feval(f,x);
end
s=h*(feval(f,a)+feval(f,b))/2+h*s;

function y=fh(x)
y=1/(-2/9+(10/9)*x-8/9*exp(4*(x-1)));

Thanks!

[Sargondjani]

can you be (alot) more specific?

what did you try? (your code)

what (do you think) should be the answer?

[Walter Roberson]

If 0 <= H <= 1 and the integrating range is [1,H], then is there a reason you did not describe the integrating range as being [H,1] ?

Answer 1

[Mike Hosea]

You can't just integrate over the singularity near x=0.66. That's why you're getting the wrong answer.

[u]

If I want to integral from zero to one, should I except singularity x=0.66,thank you!

[Mike Hosea]

I don't know what you mean by "except". You can integrate from 1 down to some number larger than r, where x=r is the location of the singularity. Since the numerator is 1, r is a root of the denominator, easily found to good accuracy with FZERO. The singularity is a barrier. It is too strong for QUADGK to handle (I tried what Titus suggested before I responded the first time). The linear term in the denominator troubles me--I don't think it's an integrable singularity.

Answer 2

[Titus Edelhofer]

Hi,

I would suggest to compute the position of the singularity (e.g. using fzero) and try quadgk on the left integral and the right integral separately.

Titus