请教大牛们一个关于quad2d的问题

town1998

本人最近在做多重积分，从网上查到了rocwoods大牛从二重到多重的一整套解决方案，可惜有一个问题，网上查到的被积函数都是用匿名函数表示的，而我的问题无法用匿名函数表达，只能用嵌套函数这类的，但是对于一元函数，写成嵌套函数我试成功了，对于二元函数我就不知道怎么弄了。

具体来说，我的被积函数是一个方波函数，比如一重时，有一个函数y(x)，只有x处于0到1之间时函数为1，其它x值则有y=0，可以如下计算
f=quad(@(x)y(x),-10,10);
    function result=y(x)
        result=zeros(size(x));
        for k=1:length(x)
        if x(k)<=1 && x(k)>=0
                          result(k)=1;
            else
                result(k)=0;
      
        end
        end
    end

但是两重时，就不知道怎么回事算不出来了，比如我的函数是z(x,y)，只有当x和y均处于0~1之间时才有z=1，否则z=0，程序是否是直接从一重推广过来，像下面这样写？
f=quad2d(@(x,y)z(x,y),-10,10,-10,10);
function result=z(x,y)
result=zeros(size(x),size(y));
for k=1:length(x)
        for l=1:length(y)
                     if x(k)>=0 && x(k)<=1 && y(l)>=0 && y(l)<=1
                         result(k,l)=1;
                     else
                         result(k,l)=0;
                     end
                 end
             end
          end

这样运行老是显示Integrand output size does not match the input size.，有人知道怎么办吗？

winner245

你应该仔细看看 quad2d 的帮助文档，你便不难发现你的被积函数 z(x,y) 定义有问题

q = quad2d(fun,a,b,c,d) approximates the integral of fun(x,y)

All input functions must be vectorized. The function Z=fun(X,Y) must accept 2-D matrices X and Y of the same size and return a matrix Z of corresponding values.

可见 quad2d 要求 Z=fun(X,Y) 中 X、Y、Z 是同尺寸的矩阵，而按照你目前的定义，你是将X、Y看作了向量，并且，你的输出 Z 的尺寸和X、Y不同，这是造成错误的根本原因。按照这个要求，你的代码可以修改为：

f=quad2d(@(x,y)z(x,y),-10,10,-10,10)
function result=z(x,y)
    result=zeros(size(x));
    for k=1:size(x,1)
        for l=1:size(x,2)
            if x(k,l)>=0 && x(k,l)<=1 && y(k,l)>=0 && y(k,l)<=1
                result(k,l)=1;
            else
                result(k,l)=0;
            end
        end
    end
end

或者用线性索引，修改为：

f=quad2d(@(x,y)z(x,y),-10,10,-10,10)
function result=z(x,y)
    result=zeros(size(x));
    for k=1:numel(x)
        if x(k)>=0 && x(k)<=1 && y(k)>=0 && y(k)<=1
            result(k)=1;
        else
            result(k)=0;
        end
    end
end

最后，我想说的是，你一开始说 “而我的问题无法用匿名函数表达，只能用嵌套函数这类的”，这显然是不成立的。如果可以用嵌套函数表达的话，一般也是可以用匿名函数表达，至少，你这里的被积函数可以轻松用匿名函数实现，而且你的函数定义没有必要用循环实现，应该充分利用向量化计算的优点：

f = @(x,y) (0<=x & x<=1 & 0<=y & y<=1)*1;
quad2d(f,-10,10,-10,10)

town1998

winner245 发表于 2014-6-1 17:54
你应该仔细看看 quad2d 的帮助文档，你便不难发现你的被积函数 z(x,y) 定义有问题

q = quad2d(fun,a,b,c,d ...

非常感谢，说的非常清楚明白，不但解决了我的问题，还从原理上进行了解释，还优化了我的程序，谢谢。

mxlzhenzhu

winner245 发表于 2014-6-1 17:54
你应该仔细看看 quad2d 的帮助文档，你便不难发现你的被积函数 z(x,y) 定义有问题

q = quad2d(fun,a,b,c,d ...

我遇到了一个非常扯淡的问题，这是在2013版本上发现的，同样的输入，在2010上计算结果是正确的；想知道到底哪里出问题了？


2013版本下面用下面计算结果不一样：
syms x
y=@(x)(0.045./(0.045-x)).^1.3-1;
eval(int(y,0,0.04))
计算结果为-0.1;在2013下面数值积分调用quad(y,0,0.04)得到是0.1；

这个结果很让人郁闷；

如果用2010版本计算结果都是0.1，这个就正确了；

所以搞不清int符合积分在2013下面怎么有问题了？在2013环境下面测试简单的函数，定积分得到的结果也是正确的，偏偏这个表达式有问题，不晓得是2013语法有更新？


【最近为啥上传不了图了？】

winner245

mxlzhenzhu 发表于 2014-6-25 17:07
我遇到了一个非常扯淡的问题，这是在2013版本上发现的，同样的输入，在2010上计算结果是正确的；想知道到 ...

你的提问方式应该再具体些，否则，别人不可能完全理解你的问题

winner245

mxlzhenzhu 发表于 2014-6-25 17:07
我遇到了一个非常扯淡的问题，这是在2013版本上发现的，同样的输入，在2010上计算结果是正确的；想知道到 ...

正确结果确实应该是 0.1 （与手算吻合）
这应该是 MuPad 的一个bug，其实MuPad在计算相应的不定积分的时候就已经出错了。

mxlzhenzhu

winner245 发表于 2014-6-26 21:55
正确结果确实应该是 0.1 （与手算吻合）
这应该是 MuPad 的一个bug，其实MuPad在计算相应的不定积分的时候 ...

非常感谢，我说嘛，还是老版本的好，升级要谨慎啊；

winner245

mxlzhenzhu 发表于 2014-6-25 17:07
我遇到了一个非常扯淡的问题，这是在2013版本上发现的，同样的输入，在2010上计算结果是正确的；想知道到 ...

这个 bug 主要是由于MuPad在计算 -1.3 次方时 出了问题，似乎 MuPad 对于 (-1/x)^1.3 和 (-x)^(-1.3) 的处理完全不同。具体可能跟符号引擎内部的规则有关  

1. syms x
   
2. int((-1/x)^1.3,x)
   
3. int((-x)^(-1.3),x)

复制代码

ans =
-(10*(-1/x)^(3/10))/7

ans =
10/(3*(-x)^(3/10))


定积分计算的差异更最直观：

1. double(int((-1/x)^1.3,x,1,2))
   
2. double(int((-x)^(-1.3),x,1,2))

复制代码

ans =
   0.1577 + 0.2170i

ans =
  -0.3679 + 0.5063i