Prev: Variable names to struct fields and back without file I/O?
Next: uigetdir default position
From: Angie on 25 Jun 2010 18:24
Hi,

Is there a way to solve the following definite integral? My integrand is symbolic:

myintegrand = @(x) x*normcdf(x,1,0);
int(myintegrand,x,0,1)

Matlab says:
Warning: Explicit integral could not be found.
ans =
int(x*normcdf(x, 1, 0), x = 0..1)

Basically, it does not give me an answer. I also tried "erfc" with no luck.

Thank you,

A.
From: Angie on 25 Jun 2010 19:14
"Angie" wrote in message <i03aa5$t52$1(a)fred.mathworks.com>...
> Hi,
>
> Is there a way to solve the following definite integral? My integrand is symbolic:
>
> myintegrand = @(x) x*normcdf(x,1,0);
> int(myintegrand,x,0,1)
>
> Matlab says:
> Warning: Explicit integral could not be found.
> ans =
> int(x*normcdf(x, 1, 0), x = 0..1)
>
> Basically, it does not give me an answer. I also tried "erfc" with no luck.
>
> Thank you,
>
> A.

First: In the above post, I should have normcdf(x,0,1) not normcdf(x,1,0)

Second: I tried this approach in another post, yet I am not sure if this is correct. Can anyone tell me if this is correct?

fun1 = @(x) x;
fun2 = @(x) normcdf(x,0,1);
myintegrand = @(x) fun1(x).*fun2(x);
quad(myintegrand,0,1)

ans =
0.3710
From: Roger Stafford on 25 Jun 2010 20:19
"Angie" <angie11tr(a)yahoo.com> wrote in message <i03d7t$3tu$1(a)fred.mathworks.com>...
> First: In the above post, I should have normcdf(x,0,1) not normcdf(x,1,0)
>
> Second: I tried this approach in another post, yet I am not sure if this is correct. Can anyone tell me if this is correct?
>
> fun1 = @(x) x;
> fun2 = @(x) normcdf(x,0,1);
> myintegrand = @(x) fun1(x).*fun2(x);
> quad(myintegrand,0,1)
>
> ans =
> 0.3710
- - - - - - - - -
Yes, that agrees with the answer I got using integration by parts with old-fashioned pen and paper, namely:

1/4 + 1/2/sqrt(2*pi*exp(1)) = 0.370985362

I would have thought the symbolic toolbox could also solve this.

Roger Stafford
From: Angie on 26 Jun 2010 16:36
"Roger Stafford" <ellieandrogerxyzzy(a)mindspring.com.invalid> wrote in message <i03h1r$m51$1(a)fred.mathworks.com>...
> "Angie" wrote in message <i03d7t$3tu$1(a)fred.mathworks.com>...
> > First: In the above post, I should have normcdf(x,0,1) not normcdf(x,1,0)
> >
> > Second: I tried this approach in another post, yet I am not sure if this is correct. Can anyone tell me if this is correct?
> >
> > fun1 = @(x) x;
> > fun2 = @(x) normcdf(x,0,1);
> > myintegrand = @(x) fun1(x).*fun2(x);
> > quad(myintegrand,0,1)
> >
> > ans =
> > 0.3710
> - - - - - - - - -
> Yes, that agrees with the answer I got using integration by parts with old-fashioned pen and paper, namely:
>
> 1/4 + 1/2/sqrt(2*pi*exp(1)) = 0.370985362
>
> I would have thought the symbolic toolbox could also solve this.
>
> Roger Stafford

Thank you. Yes, I thought symbolic integration could be possible, but it was not. I also checked it using Riemann sum, it seems to be correct. Thanks again for confirming it.
From: Roger Stafford on 26 Jun 2010 21:00
"Angie" <angie11tr(a)yahoo.com> wrote in message <i05obl$2pf$1(a)fred.mathworks.com>...
> Thank you. Yes, I thought symbolic integration could be possible, but it was not. I also checked it using Riemann sum, it seems to be correct. Thanks again for confirming it.
- - - - - - - - - - - -
With a Riemann sum or any of matlab's numerical quadrature routines you can only get an approximate value for your integral. However, as I mentioned earlier, using the integration by parts technique this iterated integral can be reduced to the exact expression

1/4 + 1/2/sqrt(2*pi*exp(1))

Below is a very brief sketch of the reasoning involved in case it is of interest to you.

We start with your expression

int(x*int(k*exp(-y^2/2),'y',-inf,x),'x',0,1)

where k = 1/sqrt(2*pi). Using integration by parts, integrate x which gives x^2/2 and differentiate

int(k*exp(-y^2/2),'y',-inf,x)

which gives

k*exp(-x^2/2).

Their product can be written

(x/2)*(k*x*exp(-x^2/2)),

and next differentiate the x/2 which gives 1/2 and integrate the

k*x*exp(-x^2/2)

which gives

-k*exp(-x^2/2).

After all the smoke clears away and everything possible is simplified, the above exact expression remains. At one point we make use of the fact that

int(k*exp(-x^2/2),'x',-inf,0) = 1/2

Roger Stafford
 |  Next  |  Last
Pages: 1 2
Prev: Variable names to struct fields and back without file I/O?
Next: uigetdir default position