Problem with symbolic integration

2012 12 15
1 个回答
回答已采纳
2 次查看（30 天）

0 个投票

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!

4 个评论
显示 2更早的评论 隐藏 2更早的评论

Walter Roberson
Walter Roberson 2012-12-15
Please show enough code to reproduce the problem.
ZhiH
ZhiH 2012-12-15
All the code is showed there.thank you
Walter Roberson
Walter Roberson 2012-12-15
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
ZhiH 2012-12-15
编辑：ZhiH 2012-12-15
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

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2012-12-15

0 个投票

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.

3 个评论
显示 1更早的评论 隐藏 1更早的评论

ZhiH
ZhiH 2012-12-15
编辑：ZhiH 2012-12-15
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
Walter Roberson 2012-12-15
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
ZhiH 2012-12-15
Thank you.I'll check this code again.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by