Octave中一元自定义函数二重数值积分实现问题求助

我在Octave中尝试对一个一元自定义函数进行二重数值积分时遇到困难，具体情况如下：

自定义函数与单次积分表现 #

自定义函数Sigma1代码如下：

function result = Sigma1(ec)
    ec2 = -0.001;
    fcd = 20e6;
    if(ec >= 0)
        result = 0;
        return;
    endif
    if(ec > ec2)
        result = -fcd.*(1-(1-ec./ec2).^2);
        return;
    endif
    result = -fcd;
endfunction

单次积分可正常运行，执行integral(@Sigma1, 0, -0.003)得到结果6.0000e+04。

错误尝试1：嵌套调用integral #

嵌套调用integral(integral(@Sigma1, 0, -0.003), 0, -0.003)时触发错误：

error: feval: first argument must be a string, inline function, or a function handle
error: called from
    integral at line 123 column 36

错误尝试2：使用integral2 #

执行integral2(@Sigma1, 0, -0.003, 0, -0.003)报错“Sigma1: function called with too many inputs”；于是修改为接受两个参数的Sigma2：

function result = Sigma2(ec, foo)
    ec2 = -0.001;
    fcd = 20e6;
    if(ec >= 0)
        result = 0;
        return;
    endif
    if(ec > ec2)
        result = -fcd.*(1-(1-ec./ec2).^2);
        return;
    endif
    result = -fcd;
endfunction

执行integral2(@Sigma2, 0, -0.003, 0, -0.003)得到结果-180，不符合预期。

问题原因与正确解法 #

1. 嵌套integral调用的错误本质 #

integral(@Sigma1, 0, -0.003)返回的是数值结果（6e4），但integral的第一个参数必须是函数句柄、字符串或内联函数，把数值传给它自然会触发参数类型错误。

2. integral2使用错误的原因 #

integral2要求被积分函数是真正的二元函数——需要接受两个独立变量并输出对应点的函数值。你写的Sigma2只是把第二个参数当占位符，实际只用到第一个参数，相当于计算：
$$\int_{0}^{{-0.003}\int_{0}}{-0.003} \Sigma1(ec) , dfoo , dec$$
内层对foo积分时，$\Sigma1(ec)$是常数，积分区间长度为$-0.003$，所以结果是$6e4 \times (-0.003) = -180$，这和你想要的二重积分逻辑不符。

正确实现方式 #

根据实际需求选择对应方法：

• 如果是两次独立积分的乘积（即单次积分结果的平方）：

  single_res = integral(@Sigma1, 0, -0.003);
  double_res = single_res ^ 2;
  

  结果为$(6e4)^2 = 3.6e9$。

• 如果是累次积分（内层积分上限依赖外层变量）：
  比如计算$\int_{0}^{-0.003} \left( \int_{0}^{y} \Sigma1(x) , dx \right) dy$，需要构造嵌套函数句柄：

  double_res = integral(@(y) integral(@Sigma1, 0, y), 0, -0.003);
  

内容的提问来源于stack exchange，提问作者stenio