Error 'DOUBLE cannot convert the input expression into a double array'

I have symbolically integrated as following

function[flux]= superquad2(a1,a2,a3, epsilon1,epsilon2,epsilon3)
 syms eta w 
 wmin=-pi/2;
 wmax=pi/2;
 etamin=-pi;
 etamax=pi;
 elp_x = a1.*abs(cos(eta)).^epsilon1.*abs(cos(w)).^epsilon1.*sign(cos(eta)).*sign(cos(w));
 elp_y = a2.*abs(cos(eta)).^epsilon2.*abs(sin(w)).^epsilon2.*sign(cos(eta)).*sign(sin(w));
 elp_z = a3.*abs(sin(eta)).^epsilon3.*sign(sin(eta));
 ellipsoid=[elp_x,elp_y,elp_z];
 ndseta=diff(ellipsoid,eta);
 ndsw=diff(ellipsoid,w);
 nds=cross(ndseta,ndsw);
 F=[(a1.*abs(cos(eta)).*abs(cos(w))),(a2.*abs(cos(eta)).*abs(sin(w))),(a3.*abs(sin(eta)))];
 F_par=subs(F,[(a1.*sign(cos(eta)).*abs(cos(eta)).*sign(cos(w)).*abs(cos(w))),(a2.*sign(cos(eta)).*abs(cos(eta)).*sign(sin(w)).*abs(sin(w))),(a3.*sign(sin(eta)).*abs(sin(eta)))],ellipsoid);
 realdot=@(u,v)u*transpose(v);
 fflux=realdot(F_par,nds);
 fluxxx=symint2(fflux,w,wmin,wmax,eta,etamin,etamax);
 fluxxxx=eval(fluxxx);
 flux=double(fluxxxx);
end

To test the stability of the code, when I call this function with integer exponents it converges to a solution, but when I try with this

superquad2(25,25,25,.33,.66,.33)

the function returns the following error

Error using symengine
DOUBLE cannot convert the input expression into a double array.

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

Error in superquad2 (line 21)
 flux=double(fluxxxx);

I could not understand it. I am starting to debug it but it doesn't help much.

Please help me to resolve the issue.

[Adam]

What is symint2? And would you expect whatever it outputs, having gone through eval, which, frankly could produce anything, to be convertable to double?

[Hirak]

I have attached the function herewith.

It should return some numerical values as the solution.

[Adam]

So why are you running the result through eval if it returns numerical results? I have no idea what eval does to numerical data as I never use it, but I don't see why it is used here.

[Hirak]

It was used to evaluate the functional form obtained symbolically.

[Stephen23]

"It was used to evaluate the functional form obtained symbolically."

eval is not the right tool for that. It is not even listed as a symbolic toolbox function:

https://www.mathworks.com/help/symbolic/functionlist.html

But that list does have a section entitled "Conversion Between Symbolic and Numeric", which includes double and vpa. Did you try either of the fucntions that are actually listed in the symbolic toolbox documenation?

[Walter Roberson]

you should never eval something that is symbolic . The expression language used by the symbolic toolbox is not exactly the same as MATLAB .

[Hirak]

even I desect that portion though the error continues to return.

Answer 1

[Walter Roberson]

you should be using vpaintegral instead of int for that code.

double() of an int() forces int to try numeric integration, but with a different strategy than vpaintegral .

numeric integration through double() of int() or through vpaintegral() can fail such as if the function diverges or does not converge like int of sin to infinity . When you use double the only way to deal with that is try/catch. When you use vpaintegral it will reply back with aa vpaintegral when it is not able to calculate and you can look at the length of children() to determine whether you got back aa numeric result .

[Hirak]

I tried with vpaintegral.The code returns an error

>> superquad2(25,25,25,.33,.66,.33)
Struct contents reference from a non-struct array object.

Error in sym/vpaintegral (line 196)
                   [x.s '=' a.s '..' b.s],...

Error in superquad2 (line 20)
 flux=vpaintegral(ffluxx,w,wmin,wmax,eta,etamin,etamax);

[Walter Roberson]

vpaintegral is not a two dimensional .you need nested vpaintegral just like you use nested int()

Answer 2

[Hirak]

I have tried with this

syms x y z eta w
 %surface
 a1=25; a2=25; a3=25; epsilon1=.33; epsilon2=.66; epsilon3=.33;
 elp=[a1.*abs(cos(eta)).^epsilon1.*abs(cos(w)).^epsilon1.*sign(cos(eta)).*sign(cos(w)),a2.*abs(cos(eta)).^epsilon2.*abs(sin(w)).^epsilon2.*sign(cos(eta)).*sign(sin(w)),a3.*abs(sin(eta)).^epsilon3.*sign(sin(eta))];
 ndseta=diff(elp,eta);
 ndsw=diff(elp,w);
 ndS=cross(ndseta,ndsw);
 
 %force
 F=[x,y,z];
 
 %vector
 ndseta=diff(elp,eta);
 ndsw=diff(elp,w);
 nds=cross(ndseta,ndsw);
 
 F_par=subs(F,[x,y,z],elp);
 realdot=@(u,v)u*transpose(v);
 
 flux1=realdot(F_par,nds);
 flux2=symint2(flux1,w,0,2*pi,eta,0,pi);
 flux3=double(flux2);

it returns again the error

Error using symengine
DOUBLE cannot convert the input expression into a double array.

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

Error in Untitled7 (line 15)
 flux3=double(flux2);
 

although I am obtaining a functional form, the numericall result could not be obtained.

What is the point? Why Double operation retun values for integer exponents (when we have used the values 2,2 and 2 for epsilon1, epsilon2 and epsilon3, respectively), whereas for fractional values (we have used the values .33,.66 and .33 for epsilon1, epsilon2 and epsilon3, respectively) it does not?

[Walter Roberson]

Look at the denominator of flux1 . It is 0 under several circumstances:

1. eta = 0
2. eta = pi/2
3. eta = pi
4. w = pi/2
5. w = pi
6. w = 3*pi/2
7. w = 2*pi

You can show that the function crosses between -inf and +inf multiple times. These are unremovable singularities, and leave the integral undefined (because it becomes the sum of positive and negative infinities.)

[Hirak]

Thank you, I can understand now. I went to the kernel of my calcuations, which is

syms x y z w eta 
a1=25; a2=25; a3=25; epsilon1=.33; epsilon2=.33; epsilon3=0.66;
ellipsoid=[(a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1),(a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.*sign(sin(w)).* abs(sin(w)).^epsilon2),(a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3)];
F = [(a1.*abs(cos(eta)).*abs(cos(w))),(a2.*abs(cos(eta)).*abs(sin(w))),(a3.*abs(sin(eta)))];
ndS = cross(diff(ellipsoid,w),diff(ellipsoid,eta));
Fpar = subs(F,[(a1.*abs(cos(eta)).*abs(cos(w))),(a2.*abs(cos(eta)).*abs(sin(w))),(a3.*abs(sin(eta)))],ellipsoid);
realdot = @(u,v) u*transpose(v);
normvec=realdot(Fpar,ndS);
flux = vpaintegral(vpaintegral(normvec,eta,-pi/2,pi/2),w,-pi,pi);

Problem is, when (eps1,eps2,eps3) is set to (1,1,1) or (1,1,0.2) the solution converges. Now, I am trying with the configuration (.33,.33,.66). the code is running for the past ten hours, but could not result a solution.

Is it really running? What makes it so slow?

[Walter Roberson]

The denominator is 0 if:

1. eta = -pi/2
2. eta = 0
3. eta = pi/2
4. w = -pi
5. w = -pi/2
6. w = 0
7. w = pi/2
8. w = pi

So again you have unremovable singularities leading to sums of negative and positive infinities.

[Hirak]

WHAT TO DO? WHAT TO DO? I cannot replace my functions.

[Walter Roberson]

You have a few choices:

1. carefully control your integrations to not cross any of the singularities. As the singularties are not removable, then how close you get to the singularities will be important, as you will be adding and subtracting values approaching infinity.
2. Check your equations again multiple times, paying particular attention to denominators and how they arise and why they can be zero. Perhaps you will find a reason why those situations do not really arise and can be avoided for your situation
3. Check your equations again multiple times, paying particular attention to denominators and how they arise and why they can be zero. Perhaps you will find a reason why your situation is not solvable. For example, you have abs() and diff(abs) is not defined which tells you that your problem is not well posed.
4. You could give up.