Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
A Strange problem on numerical Integration using Quad

Subject: A Strange problem on numerical Integration using Quad

From: Oriole

Date: 1 Oct, 2009 12:05:05

Message: 1 of 4

Hi, I want to use a complicated formula for calculating some ruin probability. In this formula there is a integral term that needs to be evaluated numerically.

Here is my program to evaluate the integral:

function testquad(eta, b, u)
    Q = quad(@(x)myfun(x,eta,b,u),0,Inf)
    ezplot(@(x)myfun(x,eta, b,u))
    
function y = myfun(x, eta, b, u)
     b1=1/b;
     y=(x.^b1 .*exp(-(x+1).*u.*b1))./((x.^b1.*(1+(1+eta).*(x+1).*b1)-cos(pi.*b1)).^2+ (sin(pi.*b1)).^2);

Then in the command window:
>> eta=0.1;b=10;u=20;
>> testquad(eta, b, u)

Q =

   NaN

Ok, maybe because "quad" can not evaluate with infinite bound. So I tried to replace the Inf in the above code (Q = quad(@(x)myfun(x,eta,b,u),0,Inf)) by numbers like 5, 10, 20 ..and up to 81. The "Q" I get is always
 Q = 0.4965

But when I change into:
Q = quad(@(x)myfun(x,eta,b,u),0,82)

That is, I execute the following program:

function testquad(eta, b, u)
    Q = quad(@(x)myfun(x,eta,b,u),0,82)
    ezplot(@(x)myfun(x,eta, b,u),[81, 82])
    
function y = myfun(x, eta, b, u)
     b1=1/b;
     y=(x.^b1 .*exp(-(x+1).*u.*b1))./((x.^b1.*(1+(1+eta).*(x+1).*b1)-cos (pi.*b1)).^2+ (sin(pi.*b1)).^2);

Then the output is very strange:
>> testquad(eta, b, u)

Q =

  1.4370e-005

However in the ezplot, it shows that the function is smooth and declining on the interval [81 82]. How to explain this strange situation? what is wrong here?

And the last question: what is the best to integrate "myfun" from 0 to infty?

It is almost the first time that I do numerical integration, really confused.. please give some advices. Thank you very very much!
Oriole

Subject: A Strange problem on numerical Integration using Quad

From: Steven Lord

Date: 1 Oct, 2009 13:16:47

Message: 2 of 4


"Oriole " <oriole_ni@hotmail.com> wrote in message
news:ha25th$8tu$1@fred.mathworks.com...
> Hi, I want to use a complicated formula for calculating some ruin
> probability. In this formula there is a integral term that needs to be
> evaluated numerically.
>
> Here is my program to evaluate the integral:
>
> function testquad(eta, b, u)
> Q = quad(@(x)myfun(x,eta,b,u),0,Inf)

This will not work. QUAD subdivides the interval over which you want to
integrate and evaluates your integrand at the subdivision points. If the
interval is infinite, at least one of the subintervals must be infinite.

> ezplot(@(x)myfun(x,eta, b,u))
>
> function y = myfun(x, eta, b, u)
> b1=1/b;
> y=(x.^b1
> .*exp(-(x+1).*u.*b1))./((x.^b1.*(1+(1+eta).*(x+1).*b1)-cos(pi.*b1)).^2+
> (sin(pi.*b1)).^2);
>
> Then in the command window:
>>> eta=0.1;b=10;u=20;
>>> testquad(eta, b, u)
>
> Q =
>
> NaN
>
> Ok, maybe because "quad" can not evaluate with infinite bound.

That's right. The QUADGK function, however, can.

> So I tried to replace the Inf in the above code (Q =
> quad(@(x)myfun(x,eta,b,u),0,Inf)) by numbers like 5, 10, 20 ..and up to
> 81. The "Q" I get is always
> Q = 0.4965
>
> But when I change into:
> Q = quad(@(x)myfun(x,eta,b,u),0,82)
>
> That is, I execute the following program:
>
> function testquad(eta, b, u)
> Q = quad(@(x)myfun(x,eta,b,u),0,82)
> ezplot(@(x)myfun(x,eta, b,u),[81, 82])
>
> function y = myfun(x, eta, b, u)
> b1=1/b;
> y=(x.^b1 .*exp(-(x+1).*u.*b1))./((x.^b1.*(1+(1+eta).*(x+1).*b1)-cos
> (pi.*b1)).^2+ (sin(pi.*b1)).^2);
>
> Then the output is very strange:
>>> testquad(eta, b, u)
>
> Q =
>
> 1.4370e-005

This type of 'issue' has come up in CSSM a number of times in the past:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/254970

http://www.mathworks.com/matlabcentral/newsreader/view_thread/246456

http://www.mathworks.com/matlabcentral/newsreader/view_thread/122969#310160

http://www.mathworks.com/matlabcentral/newsreader/view_thread/151640

You've "zoomed out" of the region of interest too far for QUAD to "see what
the function looks like".

> And the last question: what is the best to integrate "myfun" from 0 to
> infty?

Integrate it over the region where it's sufficiently different from zero
instead of from 0 to Inf and/or use QUADGK.

--
Steve Lord
slord@mathworks.com
comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ

Subject: A Strange problem on numerical Integration using Quad

From: Oriole

Date: 1 Oct, 2009 13:33:03

Message: 3 of 4

"Steven Lord" <slord@mathworks.com> wrote in message <ha2a36$bd5$1@fred.mathworks.com>...
> You've "zoomed out" of the region of interest too far for QUAD to "see what
> the function looks like".
>
> > And the last question: what is the best to integrate "myfun" from 0 to
> > infty?
>
> Integrate it over the region where it's sufficiently different from zero
> instead of from 0 to Inf and/or use QUADGK.
>
> --
> Steve Lord
> slord@mathworks.com
> comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ


Thank you Steve, I have read these posts, I see the point now. I think I will first evaluate the integrand at certain points so I can find the "significnat" region to integrate over.

Oriole

Subject: A Strange problem on numerical Integration using Quad

From: Oriole

Date: 1 Oct, 2009 13:43:04

Message: 4 of 4

"Oriole " <oriole_ni@hotmail.com> wrote in message <ha2b2f$fan$1@fred.mathworks.com>...
> "Steven Lord" <slord@mathworks.com> wrote in message <ha2a36$bd5$1@fred.mathworks.com>...
> > You've "zoomed out" of the region of interest too far for QUAD to "see what
> > the function looks like".
> >
> > > And the last question: what is the best to integrate "myfun" from 0 to
> > > infty?
> >
> > Integrate it over the region where it's sufficiently different from zero
> > instead of from 0 to Inf and/or use QUADGK.
> >
> > --
> > Steve Lord
> > slord@mathworks.com
> > comp.soft-sys.matlab (CSSM) FAQ: http://matlabwiki.mathworks.com/MATLAB_FAQ
>
>
> Thank you Steve, I have read these posts, I see the point now. I think I will first evaluate the integrand at certain points so I can find the "significnat" region to integrate over.
>
> Oriole

Hi, I just replaced quad by quadgk, and it is working well, even for integrating from 0 to Inf. /Oriole

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us