Problem with symbolic integration

Hello there!

I do have a problem with symbolic integration. I'm calculating a volume,but I can't even get a result. The original code looks like this:

 function f = fmbloss(phi)
n1 = 1.49;
n2 = 1.485;
ns1 = n1*n1;
ns2 = n2*n2;
lambda = 0.58e-6;
ra1 = 980e-6;
kw = 2*pi/lambda;
rad = 5e-3;
r0 = rad;
syms theta r
theta1 = acos(r0*cos(theta)/(rad+ra1));
theta2 = acos((rad+ra1)*cos(theta1)/(rad-ra1));
phip = 2*(theta1-theta2);
invarl = n1*cos(theta1);
invarls = invarl.*invarl;
p = invarl/n2;
ps = p.*p;
ln1 = sqrt(ns1-invarls);
ln2 = sqrt(invarls-ns2);
btf = 4*ln1.*ln2/(ns1-ns2);
bt = btf.*exp(-2*kw*ra1*invarl.*(log(p+sqrt(ps-1))-sqrt(ps-1)./p));
gamma = bt./phip
thetac = acos(n2/n1);
f1 = exp(-gamma*phi);
int1 = int(f1,theta,-thetac,thetac);
int2 = int(int1,r,rad-ra1,rad+ra1);
f = eval(int2);

and input fplot(@fmbloss,[0,pi/2]),grid on in the command window,

but

    Warning: Explicit integral could not be found. 
Warning: Explicit integral could not be found. 
Undefined function 'isfinite' for input arguments of type 'sym'.
Error in fplot (line 113)
  J = find(isfinite(y));

It is not working at all.I hope someone can give me a good advise, how to deal with that. Thank you!

[Walter Roberson]

Please show enough code to reproduce the problem.

[ZhiH]

All the code is showed there.thank you

[Walter Roberson]

Okay, please put a breakpoint in at the int, and show me the values of thetac and f1

Also: you should not eval() a symbolic formula. If you are expecting a definite numeric value from it, use double(int2)

[ZhiH]

the result is:

    >> fplot(@fmbloss,[0,pi/2]),grid on 
gamma = 
-(9007199254740992*exp(-(2174072679298543475*cos(theta)*(log((74500*cos(theta))/88803 + ((5550250000*cos(theta)^2)/7885972809 - 1)^(1/2)) - (88803*((5550250000*cos(theta)^2)/7885972809 - 1)^(1/2))/(74500*cos(theta))))/82188494176256)*(22201/10000 - (555025*cos(theta)^2)/357604)^(1/2)*((555025*cos(theta)^2)/357604 - 88209/40000)^(1/2))/(33495522228567*(2*acos((250*cos(theta))/201) - 2*acos((250*cos(theta))/299)))

 f1 = 
1 
31  int1 = int(f1,theta,-thetac,thetac);
33  int2 = int(int1,r,rad-ra1,rad+ra1);
36  f = int2;
End of function fmbloss.
K>>

and thetac=0.0819

so, Is there anything I can do now? thank you

Answer 1

[Walter Roberson]

Sorry, if there is an analytical solution for that first integral, then it is not at all easy to find. The integral can be done numerically, though.

Notice, by the way, that your gamma does not involve r (but does involve theta), so given a specific numeric phi, int1 would evaluate to a specific numeric value, so the integral of that with respect to r would just be (2*ra1) times int1.

[ZhiH]

thank you. And this is the new code:

 function f = fmbloss(phi)
n1 = 1.49;
n2 = 1.485;
ns1 = n1*n1;
ns2 = n2*n2;
lambda = 0.58e-6;
ra1 = 980e-6;
kw = 2*pi/lambda;
rad = 5e-3;
r0 = rad;
syms theta r
theta1 = acos(r0*cos(theta)/(rad+ra1));
theta2 = acos((rad+ra1)*cos(theta1)/(rad-ra1));
phip = 2*(theta1-theta2);
invarl = n1*cos(theta1);
invarls = invarl.*invarl;
p = invarl/n2;
ps = p.*p;
ln1 = sqrt(ns1-invarls);
ln2 = sqrt(invarls-ns2);
btf = 4*ln1.*ln2/(ns1-ns2);
bt = btf.*exp(-2*kw*ra1*invarl.*(log(p+sqrt(ps-1))-sqrt(ps-1)./p));
gamma = bt./phip;
thetac = acos(n2/n1);
f1 = exp(-gamma*phi);
int1 = int(f1,theta,-thetac,thetac);
int2 = 2*ra1*int1;
f = int2;

Error message:

 >> fplot(@fmbloss,[0,pi/2]),grid on
Warning: Explicit integral could not be found.
Undefined function 'isfinite' for input arguments of type 'sym'.
Error in fplot (line 113)
  J = find(isfinite(y));

How do I modify this code? Thank you.

[Walter Roberson]

There is no analytic integral for that expression.

Also, I am able to demonstrate that with those values, f1 is necessarily complex. The only values that affect it being necessarily complex are rad, ra1, ns1, and ns2. If I keep rad and ns1 and ns2 constant, then ra1 needs to be about 58 times smaller to prevent the problem. If I keep rad and ra1 and ns1 constant, then the maximum ns2 to prevent the problem is 1.245819398.

Are you expecting to be doing integration of complex values? If you are not then please recheck your constants and equations in detail.

[ZhiH]

Thank you.I'll check this code again.