Question

怎样用Matlab求复杂符号表达式的数值积分

那位高手能帮我看看这个程序哪里出错了，最后一步积分积不出来，能帮我修改一下吗，非常感谢！
clc;
clear;
a=4.15.*(10).^-6;
l=1.6.*(10).^-6;
b=62.5.*(10).^-6;
n1=1.4504;
n2=1.4447;
n3=1.4786;
k=2.*pi./l;
V=k.*a.*((n1.^2-n2.^2).^(1./2));
cn=((n1.*k).^2-(2.27645.*V-0.306.*V.^2-0.992)./a.^2).^(1./2);
r=(cn.^2-n2.^2.*k.^2).^(1./2);
ka=(n1.^2.*k.^2-cn.^2).^(1./2);
ev=2;
R=8;
syms y;
X20=(R.*10.^-3./(2.*k.^2.*n2.^2)).^(2./3).*(cn.^2+y.^2-k.^2.*n2.^2);
X2a=(R.*10.^-3./(2.*k.^2.*n2.^2)).^(2./3).*(cn.^2+y.^2-k.^2.*n2.^2.*(1+2.*a./(R.*10.^-3)));
X2b=(R.*10.^-3./(2.*k.^2.*n2.^2)).^(2./3).*(cn.^2+y.^2-k.^2.*n2.^2.*(1+2.*b./(R.*10.^-3)));
X3b=(R.*10.^-3./(2.*k.^2.*n3.^2)).^(2./3).*(cn.^2+y.^2-k.^2.*n3.^2.*(1+2.*b./(R.*10.^-3)));
X2=(2.*k.^2.*n2.^2./(R.*10.^-3)).^(2./3).*(-X2b);
X3=(2.*k.^2.*n2.^2./(R.*10.^-3)).^(2./3).*(-X3b);
Oy=2./3.*(-X2b).^(3./2)+pi./4;
f=exp(-a.*(r.^2+y.^2).^(1/2))./(r.^2+y.^2).^(1/2).*1./pi.*(X20./3).^(1./2).*besseli(1./3,2./3.*X20.^(3./2))./(X2a./3).^(1./2)./(besselk(-1./3,2./3.*X2a.^(3./2))+besselk(1./3,2./3.*X2a.^(3./2))).*(X2).^(1/2).*(X3).^(1/2)./(X2.*(cos(Oy)).^2+X3.*(sin(Oy)).^2);
J=quadl(@(t)subs(f,y,t),0,inf)

Answer 1

你的函数里面，y 是干嘛的？不给值怎么数值积分

追问

y是定义的符号变量，范围是0到inf，最后的积分是想将y转化成数值t，t的范围也是0到inf，来求解积分。但是这个积分表达式积不出来。

追答

你这写法不好，你要定义f这个函数不要用符号定义，直接新建一个文件，比如这样
建一个f.m
function z = f(y)
....... % 这里写函数 f 的计算代码，有字数限制我贴不过来

然后在主界面积分的命令是
quadl(@(y)test(y), 0, inf)

现在积分结果有问题，你检查一下这个函数有没有写错，我这边看，你的ka 和 ev 根本没用到

追问

谢谢你的答复喔，不过积分结果还是有问题，函数没错，积分范围要是缩短到（0,2），还可以积出来。ka和ev程序后面的语句会用到，这里没贴上来。