matlab求一个复杂双重积分的问题

悠若谷

n（x）是一个递减函数，积分下限改成1e-6
这里采用了recowoods前辈的quadl跟arrayfun结合的两次积分方法

1. 
   
2. k=2*pi/0.6328e-6;
   
3. R=quadl(@(z) exp(-z.^2)./((1/2)^2+z.^2),-0.5,0.5);
   
4. nx1=1.520167+0.0235811*1/(pi*2)*exp(-(1/2)^2)*R;
   
5. jifen=quadl(@(x) k*sqrt((1.520167+0.0235811*x/(pi*2).*exp(-(x/2).^2).*arrayfun(@(x) quadl(@(s) exp(-s.^2)./(s.^2+(x/2).^2),-0.5,0.5),x))^2-nx1^2),1e-6,1)
   

复制代码

??? Error using ==> mpower
Matrix must be square.
Error in ==> @(x)((1.520167+0.0235811*x/(pi*2).*exp(-(x/2).^2).*arrayfun(@(x)quadl(@(s)exp(-s.^2)./(s.^2+(x/2).^2),-0.5,0.5),x))^2-nx1^2)

Error in ==> quadl at 70
y = feval(f,x,varargin{:}); y = y(:).';
不知道为什么出错了
希望前辈们能帮忙看看

悠若谷

发现是点乘弄成乘了

1. R=quadl(@(z) exp(-z.^2)./((1/2)^2+z.^2),-0.5,0.5);
   
2. nx1=1.520167+0.0235811*1/(pi*2)*exp(-(1/2)^2)*R;
   
3. k=2*pi/0.6328e-6;
   
4. >> jifen=quadl(@(y) k*sqrt((1.520167+0.023581*y/(pi*2).*exp(-(y/2).^2).*arrayfun(@(y) quadl(@(s) exp(-s.^2)./(s.^2+(y/2).^2),-0.5,0.5),y)).^2-nx1^2),1e-6,1)
   
5. 
   
6. jifen =
   
7. 
   
8.      1.322407131543267e+006 +8.100124224041057e-003i

复制代码

这次对了 但是不知道为什么 积分结果成了复数
从理论上分析 结果不应该是复数

TBE_Legend

发现是点乘弄成乘了R=quadl(@(z) exp(-z.^2)./((1/2)^2+z.^2),-0.5,0.5);
nx1=1.520167+0.0235811*1/(pi*2)*exp(-(1/2)^2)*R;
k=2*pi/0.6328e-6;
>> jifen=quadl(@(y) k*sqrt((1.520167+0.023581*y/(pi*2).*exp(- ...
悠若谷 发表于 2009-9-24 15:03

你的被积函数在某点奇异，值无穷大。

悠若谷

前辈是怎么发现的？
这个点不知道有没有办法去除

悠若谷

我主要是要编一个function来解这么一个方程
xc就是newsolve函数的输入变量x，输出变量是要求的N
其他都已知

1. 
   
2. function N=newsolve(x)
   
3. ns=1.520167;
   
4. nf=1.548167;
   
5. k=2*pi/0.6328e-6;
   
6. nc=1;
   
7. D=2e-6;
   
8. N1=1.5436744；
   
9. dN=0.0235072;
   
10. N=fzero(@(N) quadl(@(y) k*sqrt((ns+dN*y/(pi*D).*exp(-(y/D).^2).*arrayfun(@(y) quadl(@(s) exp(-s.^2)./(s.^2+(y/D).^2),-0.5,0.5),y)).^2-N.^2),1e-10,x)-pi/4-atan(sqrt((N.^2-nc^2)/(N1^2-N.^2))+(k^2*N1*(-1.855800789988048e4))/((k^2*N1^2-k^2*N.^2)^1.5*2)),[1.52017,1.543]);
    

复制代码

运行结果报错newsolve(1.537433e-6)
??? Error using ==> fzero at 260
Function values at interval endpoints must be finite and real.
Error in ==> newdraw2 at 16
N=fzero(@(N) quadl(@(y) k*sqrt((ns+dN*y/(pi*D).*exp(-(y/D).^2).*arrayfun(@(y) quadl(@(s)
exp(-s.^2)./(s.^2+(y/D).^2),-w/D,w/D),y)).^2-N.^2),1e-10,x)-pi/4-atan(sqrt((N.^2-nc^2)/(N1^2-N.^2))+(k^2*N1*(-
这个问题琢磨了好几天 都没解决
想了很久还是决定来simwe问问
哪位前辈有空帮忙看看吧！

TBE_Legend

不清楚你的问题。 一点小小的提醒：

f(x)=1/x^2, F(x)是f(x)的不定积分。 你本来要求：F(2)-F(-2)

但你却用f(x)对x在[-2,2]上的定积分来球，是求不出来的，因为在x=0处发散。

rocwoods

发现是点乘弄成乘了R=quadl(@(z) exp(-z.^2)./((1/2)^2+z.^2),-0.5,0.5);
nx1=1.520167+0.0235811*1/(pi*2)*exp(-(1/2)^2)*R;
k=2*pi/0.6328e-6;
>> jifen=quadl(@(y) k*sqrt((1.520167+0.023581*y/(pi*2).*exp(- ...
悠若谷 发表于 2009-9-24 15:03

检查了半天发现是红色部分你丢了个1，下次小心点啊。呵呵
第一次看到有人用我帖子里的方法解决这样的复杂积分，还是比较高兴的，帖子没白写:lol
关于第二个问题，暂时没有时间考虑，等有时间了，帮你考虑一下。楼主可以试着自己考虑下，看是否有类似拼写方面的错误。有什么问题贴上来再来讨论。

悠若谷

7# rocwoods
嗯~点头，真不好意思，犯了个小错误~
其实不知道这个方法这个思维最后能不能解决我这个问题
却还是想要继续试下去
另外，还想问问 fsolve能不能将x（1）限制在实数范围内迭代？以防止它自动变成了虚数

rocwoods

fsolve本来就是求实数最优解的,出现虚数很可能是初始值选取不合理.你这个问题用fzero应该可以解出来啊,除非方程有问题或者无解

waynebuaa

答案可是  0.554162

TBE_Legend

7# rocwoods  
嗯~点头，真不好意思，犯了个小错误~
其实不知道这个方法这个思维最后能不能解决我这个问题
却还是想要继续试下去
另外，还想问问 fsolve能不能将x（1）限制在实数范围内迭代？以防止它自动变成了 ...
悠若谷 发表于 2009-9-25 20:31

虚数应该是由于开方引起的。

matlab的嵌套应该是可以解决你的问题的,你最好能先用个计算量小点的试试。

能的话，先画个图出来看看，不要着急求根。 附件是我用m@画的一个泥的这类问题的图，实际上比你的还要复杂些。

TBE_Legend

答案可是  0.554162
waynebuaa 发表于 2009-9-26 11:19

哦？ wayne兄怎么做的？

rocwoods

5楼方程我算出来的结果是:1.524606031629952
楼主程序没有问题,问题出在fzero求解区间:[1.52017,1.543]上了,
假设

1. ff = (@(N) quadl(@(y) k*sqrt((ns+dN*y/(pi*D).*exp(-(y/D).^2).*arrayfun(@(y) quadl(@(s) exp(-s.^2)./(s.^2+(y/D).^2),-0.5,0.5),y)).^2-N.^2),1e-10,x)-pi/4-atan(sqrt((N.^2-nc^2)/(N1^2-N.^2))+(k^2*N1*(-1.855800789988048e4))/((k^2*N1^2-k^2*N.^2)^1.5*2))

复制代码

ff在N=1.543上值是虚数,所以必须重新界定fzero的求解区间,方法就是用ff函数试验1.543附近的值,直到区间右侧的值是负数,或者实部是负数,虚部很小为止.我选的区间是[1.52017,1.5246095]
求解出上面的结果.

悠若谷

13# rocwoods
太好了，我画图看的那个范围大概就是1.52460左右。
您赶紧讲一下，是怎么发现问题的，
我好久都不明白是哪出的问题，这个太迫切了~

悠若谷

仔细看了下，现在明白了~
我看方程时候觉得N只要在ns跟N1之间取值就不会出现虚数了
一直以为是内部算法的问题呢
然后我其实积分上限这个数是不知道的 它跟N之间的关系还有另一个方程 所以我用到了fsolve
初值稍微大一些fsolve解出来就是复数
我直接把解设成初值 x0=[1.537433;1.5246060229];
options=optimset('display','iter');
[x,fval]=fsolve('newmy',x0,options)
Optimization terminated: relative function value changing by less
than max(options.TolFun^2,eps) and sum-of-squares of function
values is less than sqrt(options.TolFun).
x =
   1.538085169196832
   1.524606089165542
我不明白它怎么会不显示successfully
下面这个是随便选的初值 x（1）是积分上限值，x（2）是上面求的N
x0=[1.537;1.524];
options=optimset('display','iter');
[x,fval]=fsolve('newmy',x0,options)
Optimization terminated: relative function value changing by less
than max(options.TolFun^2,eps) and sum-of-squares of function
values is less than sqrt(options.TolFun).
x =
  1.538085720162842 - 0.000000014693514i
  1.524606086062150 + 0.000000000084311i

fval =
  1.0e-007 *
-0.029188718109197 - 0.214689046896666i
                  0 + 0.000000000012215i
虚数部分很小，好像可以直接忽略掉 然后把实部当做解 前辈们应该也碰到过这种情况吧 不知道你们是怎么看待这个问题的？

悠若谷

11# TBE_Legend
谢谢斑竹~这个画图的确挺管用 就是值还不够精确呢

maplelab

16# 悠若谷

maplelab

悠若谷

18# maplelab
里边的digits是什么呢？
不知道为什么这个值不一样虚数部分就不一样哦

gentlefu

18# maplelab  
里边的digits是什么呢？
不知道为什么这个值不一样虚数部分就不一样哦
悠若谷 发表于 2009-9-28 21:10

The Digits environment variable controls the number of digits that Maple uses  when calculating with software floating-point numbers