我正在将一些代码从MATLAB转换为python,并且正在努力获得一个函数,该函数以数组参数(作为参数)来使用Scipy进行集成。在
我已经将代码简化为一个在Scipy中产生相同错误的基本示例,而等效的MATLAB代码的功能与预期相同。在
我试图将长度为m的行向量参数参数参数和长度为N的列向量参数参数参数传递给被集成在另一个单独的积分参数上的函数,期望我的集成输出将具有MxN的形状。在
以下python代码会产生此错误:File "C:\Anaconda3\lib\site-packages\scipy\integrate\quadrature.py", line 196, in quadrature
if err < tol or err < rtol*abs(val):
ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()
从
^{pr2}$
等效的MATLAB代码(输出MxN形状的结果)是c = [1, 2];
r = transpose([2, 1]);
out = integral(@(p)intgd(p, r, c), -pi/2, pi/2, 'ArrayValued', true);
function intgd = intgd(p, r, c)
c_bcr = repmat(c, length(r), 1);
r_bcc = repmat(r, 1, length(c));
A = ones(length(r), length(c));
s = sin(r_bcc) - sin(p);
A(s ~= 0) = sin(c_bcr(s ~= 0).*s(s ~= 0))./(c_bcr(s ~= 0).*s(s ~= 0));
intgd = 1./fun(p, c).^2.*(A.^2);
end
function d = fun(p, c)
p1 = zeros(1, length(c));
mask = sqrt(pi./(2*c)) < 1;
p1(mask) = acos(sqrt(pi./(2*c(mask))));
d = zeros(1, length(c));
mask = abs(p) <= p1;
d(mask) = 1./(pi./(2*c(mask).^2) + cos(p));
mask = and(abs(p) > p1, abs(p) <= pi/2);
d(mask) = 1./(pi./(2*c(mask).^2) + ((cos(p1(mask)) - cos(p))./2));
end
上面的MATLAB输出out是[6.58727018139280, 0.963083280848789;
6.78600314283299, 1.05994693990888]
我不确定scipy.integrate.quadrature是如何处理通过它的对象的维度的,但想法是它应该产生相同的MxN输出。在
我知道numpy有它自己的内置广播,通常可以避免这里所示的显式广播的需要,但我不确定它如何处理m=N的数组,如本例中所示,因此我将其保持为显式。任何关于这个次要问题的建议也将是受欢迎的。在