DOUBLE cannot convert the input expression into a double array.Why?

I am using Matlab to calculate a double integral. The error is:

    Warning: Explicit integral could not be found. 
    Error using mupadmex
        Error in MuPAD command: DOUBLE cannot convert the input expression into a double
        array. If the input expression contains a symbolic variable, use the VPA function instead.

    Error in sym/double (line 827)
        Xstr = mupadmex('symobj::double', S.s, 0);

    Error in MeEqu7 (line 20)
         out2 =out2 (int(int(yi,u1,0,u2),u2, 0 ,1000))

and my code is as follows:

    clear all; 

    a=4; 
    la1=1/(pi*500^2); la2= la1*5;
    p1=25; p2=p1/25;
    sgma2=10^(-11);
    index=1;
    g=2./a;

    syms r u1 u2  u3 aa
    powe= 0;
    seta= 10^powe;

    xNor =  (u2./u1).^(-1./g) + 1  ;
    aa= -u2.^(-1./g) .*seta.*( pi.*(la1.*p1.^(g) + la2.*p2.^(g)) ).^(1./g)  ;

    yi =  exp(aa./xNor - u2 );
    out2 = double (int(int(yi,u1,0,u2),u2, 0 ,1000))

Thanks alot for your help!! I also tried vpa. but the out2 is not a number. Actually, I am expecting the out2 is a number, sothat I can use this number to plot a figure. I did some search here. but their case is not helpful for my case. Could you take a deeper look?

Answer 1

[John D'Errico]

Look at what you have created.

(int(yi,u1,0,u2)
ans =
int(exp(- u2 - 1934281311383407/(1208925819614629174706176*u2^2*(u1^2/u2^2 + 1))), u1, 0, u2)

This expression has the symbolic variable u2 in it. The error says the integral could not be found, but even if it could have been found, what you are trying to stuff into out2 is a function of u2.

Perhaps you were planning on trying to integrate u2 out given your initial comment, but regardless, you cannot convert a symbolic expression with a variable in it, into a double precision number.

Edit: (in response to the change made to your question.)

You changed one line of your code. The last line of code is now

out2 = double (int(int(yi,u1,0,u2),u2, 0 ,1000))

So, where before, you had only a single integral, you now do have a double nested integral. But, really, nothing material has changed. The double integral has no analytical solution. The inner integral fails, so the outer one fails too.

You cannot convert that result into a double precision number. Wanting something to happen in mathematics is not sufficient. Demanding that the integral exist is not sufficient, by trying to convert it to a double.

int(int(yi,u1,0,u2),u2, 0 ,1000)
ans =
int(int(exp(- u2 - 1934281311383407/(1208925819614629174706176*u2^2*(u1^2/u2^2 + 1))), u1, 0, u2), u2, 0, 1000)

Maybe your goal is to do a numerical integration, and you were thinking that telling MATLAB to convert to double precision that it would magically know what you wanted. MATLAB does not read your mind. If the input to double is not a number or an array of numbers, it cannot somehow know to do a numerical integration instead of the symbolic integration you asked for.

(Note that IF you do decide to solve this using a numerical integration, that integrals of exponentials with extremely wide limits of integration often tend to fail, so you will need to do that integration carefully.) So now lets look at the integrand of this double integral.

exp(- u2 - 1934281311383407./(1208925819614629174706176*u2.^2.*(u1.^2./u2.^2 + 1)));

First of all, the ratio

vpa('1934281311383407/1208925819614629174706176')
ans =
0.0000000016000000000000002650866688430533

is essentially 1.6e-9 in double precision. Next, basic algebra lets us simplify that mess of a kernel a bit more.

exp(- u2 - 1.6e-9./(u1.^2 + u2.^2));

You want to convert this into a double precision result, with the integral on u1 having limits of 0 and u2, and the integral on u2 having limits of 0 and 1000.

So now compare the two parts of what we find inside the exponential. We have a minus u2 term, which varies from 0 to 1000, and a second term which is essentially zero in double precision over the ENTIRE DOMAIN of your integration, compared to a number that varies from 0 to 1000. The result of that double integral will clearly be just slightly less than 1, since we know that

int(exp(-u2),0,1000) == 1 - exp(-1000)

The double integration seems to be doable though.

u1max = @(u2) u2;
fun = @(u2,u1) exp(- u2 - 1.6e-9./(u1.^2 + u2.^2));

integral2(fun,0,1000,0,u1max,'reltol',1e-13,'abstol',1e-13)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test. 
> In funfun/private/integral2Calc>integral2t at 129
  In funfun/private/integral2Calc at 9
  In integral2 at 106 
ans =
       0.999999986764709

So it looks like I cranked down on the tolerances a bit too much, but the result seems reasonable. In the end, if you want to do a numerical integration, use a numerical integration tool, not a symbolic one.

[sun]

@John D'Errico Thanks for your help here. I just re-edit my question. Yes, out2 was a function of u2. But now, it should not be. I just re edit. sorry about that.

There is double int. I am expecting the out2 is a number ,sothat I will use this number to plot a figure. Thanks again. Could you take a deeper look ?

[sun]

Dear John. First of all. thank you so much for your reply. Allow me to think about it. I will be back soon. oh, btw. Merry Christmas to you.

[sun]

Dear John, I just re-write my equation, which is the original equation. The first int out is not fails. It gives me

int((1934281311383407*exp(-1/(1/r1 + 1/r2))*exp(-(1878969922100481*pi*r2^(1/2))/147573952589676412928))/(4835703278458516698824704*r1^(1/2)*r2^(1/2)), r1 == 0..r2)

Hmm, as you see, the inside int output is still function of r1. Maybe that's the reason outside int output can not convert into a number.

In the line,(int(int(yi,u1,0,u2),u2, 0 ,1000)), 1000 means Inf, I am not sure the Inf for matlab, that's why I just give some "big" number here.

Here is the re-write equation.

a=4; 
la1=1./(pi*500^2); la2= la1.*5;
p1=25; p2=p1./25;
sgma2=10.^(-11);
index=1;
g=2./a;

syms r1 r2 aa bb cc dd
powe= 0;
seta= 10^powe;

aa = exp( -seta./(r1.^-1 + r2.^-1));
bb = (pi.*g.*( la1.*p1.^g +  la2.*p2.^g  )).^2; 
cc = exp( -pi.*(r2.^g).* ( la1.*p1.^g +  la2.*p2.^g )    );
dd = (r1.^(g-1) ) .* (r2.^(g-1) );

yi=aa.*bb.*cc.*dd;
out2 = double (int(int(yi,r1,0,r2),r2, 0 ,1000))

Thanks for your hint, John. The hint of range of int was helpful.