请教含贝塞尔函数的积分怎么做啊，我的程序哪里有错

djh713

f=300;
c=340;
l=c/f;
k=2*pi/l;
z=0;
r0=0.01*l;
syms r kr pp pe y1 y2;
r=0:0.001*l:0.01*l;
y1=besseli(1,r0*kr);
y2=besseli(0,r*kr);
f=r0*y1*y2;
pp=int(f,kr,0,k);
pe=int(f,kr,k,+inf);
Pp=double(pp);
Pe=double(pe);
R=double(r);
plot(R,Pp);
plot(R,Pe);

一段积分程序，含有零阶和一阶的贝塞尔函数，程序采用符号积分，是想做出变量Pe与Pp随R变化的曲线。哪位高手可帮一下。
这个程序运行后总是提示Warning: Explicit integral could not be found. 看Pp和Pe的值也没有具体值，而且贝塞尔函数好象计算不出来。如下，但实际是应为对应于R的一组值。

pp =

[ besseli(0,1/50*pi)-1,  int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/15000*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/7500*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/5000*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/3750*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/3000*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/2500*kr),kr=0..30/17*pi), int(17/1500*besseli(1,17/1500*kr)*besseli(0,119/15000*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,17/1875*kr),kr=0..30/17*pi),   int(17/1500*besseli(1,17/1500*kr)*besseli(0,51/5000*kr),kr=0..30/17*pi),  1/2*besseli(0,1/50*pi)^2-1/2]

eight

原帖由 djh713 于 2007-11-29 12:31 发表

                               
登录/注册后可看大图

f=300;
c=340;
l=c/f;
k=2*pi/l;
z=0;
r0=0.01*l;
syms r kr pp pe y1 y2;
r=0:0.001*l:0.01*l;
y1=besseli(1,r0*kr);
y2=besseli(0,r*kr);
f=r0*y1*y2;
pp=int(f,kr,0,k);
pe=int(f,kr,k,+inf);
...

请搜索版面，这个警告信息已经有很多很多讨论

djh713

f=300;
c=340;
l=c/f;
k=2*pi/l;
z=0;
r0=0.01*l;
r=linspace(0,0.01*l);
for i=1:length(r)
    rr=r(i);
    strp=[num2str(r0),'.*besseli(1,1.*',num2str(r0),'.*kr)*besseli(0,1.*',num2str(rr),'.*kr.)'];
    funp=inline(strp,'kr');
    pp(i)=quadl(funp,0,k);
end
plot(r,pp)

看了一些以前的帖，我把程序改了，但还是有错，实在找不出来了，帮帮忙吧

eight

原帖由 djh713 于 2007-11-29 16:23 发表

                               
登录/注册后可看大图

f=300;
c=340;
l=c/f;
k=2*pi/l;
z=0;
r0=0.01*l;
r=linspace(0,0.01*l);
for i=1:length(r)
    rr=r(i);
    strp=[num2str(r0),'.*besseli(1,1.*',num2str(r0),'.*kr)*besseli(0,1.*',num2str(rr), ...

请不要让别人通过运行你的代码才能知道出现什么问题了。发帖前，多看看置顶帖

djh713

eight 同志，你今天怎么就和我置气呀，我发什么帖都有问题了，看我今天不顺呀，你有这时间帮我解决一下多好呀。别人发也没怎样吧，我哪里招你了吗:@@


??? Error using ==> inlineeval
Error in inline expression ==> 0.0113*besseli(1,0.0113*kr.)*besseli(0,1.*0.*kr.)


??? Error: "identifier" expected, ")" found.

Error in ==> D:\ENGINEER\MATLAB\toolbox\matlab\funfun\@inline\feval.m
On line 34  ==>         INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);

Error in ==> D:\ENGINEER\MATLAB\toolbox\matlab\funfun\quadl.m
On line 60  ==> y = feval(f,x,varargin{:}); y = y(:).';

eight

自己检查一下strp即可：
strp=[num2str(r0),'.*besseli(1,1.*',num2str(r0),'.*kr)*besseli(0,1.*',num2str(rr),'.*kr.)'];

仔细看就能看出来了

ps：没有人要跟你过不去，对每一个会员，我都这样纠正他的提问，明白吗

[ 本帖最后由 eight 于 2007-11-29 16:43 编辑 ]

djh713

eight，我都要无语了，我要检查出来就不问了，到底错在哪你就不能说一下吗？不会就是你列的这个吧，要是，那还是错
:@Q

eight

原帖由 djh713 于 2007-11-29 16:47 发表

                               
登录/注册后可看大图

eight，我都要无语了，我要检查出来就不问了，到底错在哪你就不能说一下吗？不会就是你列的这个吧，要是，那还是错
:@Q

把strp最后的 点 去掉试试吧

djh713

我彻底无语了，看来这里的人都是“天才”吧，说着一些莫名其妙的话，eight同志，你到底会不会，要会就做出来试试，要不会就别说了，哪来的点，去掉什么试试，不是让我把错误提示写出来吗，那这倒底是什么错，:@Q

最后一次提问，太伤心了。在这里从来没解决过问题，到头来只有靠自己。

xjzuo

1. 问问题应该先将公式一并贴出，以便他人先理解你的问题，或确认你没有写错；
2. 问问题时请尊重给你建议的人，别人是无偿帮助你，没有义务帮你解决问题。
------------------------------------------------------------------------------------------------------------
你的问题本来很简单，只是你没有理解我的代码，而且套用不对，所以出错。
在你的代码中＂strp=＂之后加一个命令＂strcat＂即可，记得[]改为()。
还有，1.*之类的可去掉（虽然不影响计算），最关键的点运算在两个besseli（向量）之间，其余的点号均可去除。
-------------------------------------------------------------------------------------------------------
希望以后能稍微思考一下示例贴的代码，其意义是让大家提高编程能力，而不只是一个用来套用的工具。

xjzuo

结果图片也附上，你自己参考一下吧：

djh713

我做出来了，不过还是要谢谢xjzuo的帮助,还想问一下，为什么我用同样的方法计算另一个公式就可得出结果，也没加strcat.
程序前面都一样，仅积分这个变为
strp=['exp(i.*2.*pi.*',num2str(c),'.*',num2str(rr),'.^2./',num2str(l),'.*kr.^2./2).*besselj(0,2.*pi.*',num2str(f),'.*',num2str(rr),'./',num2str(c),'.*kr).*kr'];

花如月

遇事要先找自己的原因，eight的水平和贡献大家人所共知。

不要盲目怀疑，也希望对帮你解决问题的版友有起码的耐心和尊重

回答本身便意味着付出和帮助，他原本可以什么都不说的，什么都不做

有不少这样的帖子沉了下去，不是因为没人懂。

而是提问本身的问题实在太多、、、、、

[ 本帖最后由 花如月 于 2007-11-30 09:56 编辑 ]

djh713

谢谢大家的批评教育，我接受。
看来我提问的水平确实有问题，所以屡屡下沉，最终只能自己解决。昨天我发牢骚的原因也是太急了，搞了一天也没出来，即然没人回答又提问不合格干脆删帖算了。如果能解决讲一下不就更好了，大家都满意，何必一而再再而三的教育呢。
顺便说一下，我按xjzuo提的做了，还是一样的问题，我也搞不懂为啥了，也不想再问了。但我用另一个公式还是我原来的方法就可行。我晕了。自己啃吧，我就不信搞不定，也不是第一次了。