双数值积分误差

Question

在下面的脚本中，在用于计算双重积分的 for 循环中，我不断收到错误消息，我不确定为什么：

Error using  * 
Inner matrix dimensions must agree.

Error in @(t,p)-sin(t)*(G(:,1)*cos(p)+H(:,1)*sin(p))


Error in @(t,p)B{k}(t,p).*A{k}(t,p).*(V(t)-Veq).*sin(t)


Error in integral2Calc>integral2t/tensor (line 228)
        Z = FUN(X,Y);  NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

这个脚本是由另一个脚本组成的，它只需要处理一个变量 t，所以我假设我在将它调整为两个变量函数时做错了。

谢谢！

%%Calculation of dH/dt for mode m=1 for the entire sphere, NH and SH
clear all
%%Radius of photosphere
r = 6.957*(10^5); %In km
R = 1/r; %This will come in handy later

%%Call in spherical harmonic coefficients, change the 535 figure as more
%%data is added to the spreadsheets
G(:,1) = xlsread('G Coefficients.xls', 'D3:D535');
G(:,2) = xlsread('G Coefficients.xls', 'F3:F535');
G(:,3) = xlsread('G Coefficients.xls', 'I3:I535');
G(:,4) = xlsread('G Coefficients.xls', 'M3:M535');
G(:,5) = xlsread('G Coefficients.xls', 'R3:R535');
G(:,6) = xlsread('G Coefficients.xls', 'X3:X535');
G(:,7) = xlsread('G Coefficients.xls', 'AE3:AE535');
G(:,8) = xlsread('G Coefficients.xls', 'AM3:AM535');
G(:,9) = xlsread('G Coefficients.xls', 'AV3:AV535');

H(:,1) = xlsread('H Coefficients.xls', 'D3:D535');
H(:,2) = xlsread('H Coefficients.xls', 'F3:F535');
H(:,3) = xlsread('H Coefficients.xls', 'I3:I535');
H(:,4) = xlsread('H Coefficients.xls', 'M3:M535');
H(:,5) = xlsread('H Coefficients.xls', 'R3:R535');
H(:,6) = xlsread('H Coefficients.xls', 'X3:X535');
H(:,7) = xlsread('H Coefficients.xls', 'AE3:AE535');
H(:,8) = xlsread('H Coefficients.xls', 'AM3:AM535');
H(:,9) = xlsread('H Coefficients.xls', 'AV3:AV535');

%%Set function v which always remains the same
nhztoradperday = 2*pi*86400*(10^(-9));
a = 460.7*nhztoradperday;
b = -62.69*nhztoradperday;
c = -67.13*nhztoradperday;

B{1} = @(t,p) -sin(t)*(G(:,1)*cos(p) + H(:,1)*sin(p));

B{2} = @(t,p) -3*sin(t)*cos(t)*(G(:,2)*cos(p) + H(:,2)*sin(p));

B{3} = @(t,p) -1.5*(5*(cos(t)^2)-1)*sin(t)*(G(:,3)*cos(p) + H(:,3)*sin(p));

B{4} = @(t,p) -2.5*(7*(cos(t)^3)-3*cos(t))*sin(t)*(G(:,4)*cos(p) + H(:,4)*sin(p));

B{5} = @(t,p) -1.875*sin(t)*(21*(cos(t)^4)-14*(cos(t)^2)+1)*(G(:,5)*cos(p) + H(:,5)*sin(p));

B{6} = @(t,p) -2.625*cos(t)*sin(t)*(33*(cos(t)^4)-30*(cos(t)^2)+5)*(G(:,6)*cos(p) + H(:,6)*sin(p));

B{7} = @(t,p) -0.4375*sin(t)*(429*(cos(t)^6)-495*(cos(t)^4)+135*(cos(t)^2)-5)*(G(:,7)*cos(p) + H(:,7)*sin(p));

B{8} = @(t,p) -0.5625*cos(t)*sin(t)*(715*(cos(t)^6)-1001*(cos(t)^4)+385*(cos(t)^2)-35)*(G(:,8)*cos(p) + H(:,8)*sin(p));



A{1} = @(t,p) 0.5*R*cos(t)*(G(:,1)*cos(p) + H(:,1)*sin(p));

A{2} = @(t,p) 0.5*R*cos(2*t)*(G(:,2)*cos(p) + H(:,2)*sin(p));

A{3} = @(t,p) 0.125*R*cos(t)*(15*(cos(t)^2)-11)*(G(:,3)*cos(p) + H(:,3)*sin(p));

A{4} = @(t,p) 0.125*R*(-3*cos(2*t)+7*(cos(t)^4-3*sin(t)^2*cos(t)^2))*(G(:,4)*cos(p) + H(:,4)*sin(p));

A{5} = @(t,p) 0.0625*R*(21*(cos(t)^5)-(cos(t)^3)*(14+84*(sin(t)^2))+cos(t)*(1+28*(sin(t)^2)))*(G(:,5)*cos(p) + H(:,5)*sin(p));

A{6} = @(t,p) 0.0625*R*(33*(cos(t)^6)-(cos(t)^4)*(165*(sin(t)^2)+30)+90*(sin(t)^2)*(cos(t)^2)+5*cos(2*t))*(G(:,6)*cos(p) + H(:,6)*sin(p));

A{7} = @(t,p) 1/128*R*(429*(cos(t)^7)-(cos(t)^5)*(495+2574*(sin(t)^2))+(cos(t)^3)*(135+1980*(sin(t)^2))-cos(t)*(5+270*(sin(t)^2)))*(G(:,7)*cos(p) + H(:,7)*sin(p));

A{8} = @(t,p) 1/128*R*(715*(cos(t)^8)-1001*(cos(t)^6)+385*(cos(t)^4)-35*(cos(t)^2)+(sin(t)^2)*(-5005*(cos(t)^6)+5005*(cos(t)^4)-1155*(cos(t)^2)+35))*(G(:,8)*cos(p) + H(:,8)*sin(p));

V = @(t) a + (b*cos(t)^2) + (c*cos(t)^4);
Veq = V(0);

intNH = zeros(length(G),9);
intSH = zeros(length(G),9);
intSun = zeros(length(G),9);
for k=1:8
   fun{k} = @(t,p) B{k}(t,p).*A{k}(t,p).*(V(t)-Veq).*sin(t);
   intNH(:,k) = integral2(fun{k},0,pi/2,0,2*pi);
   intSH(:,k) = integral2(fun{k},pi/2,pi,0,2*pi);
   intSun(:,k) = integral2(fun{k},0,pi,0,2*pi);
end

for i=1:length(G)
   NH(i) = sum(intNH(i,:));
   SH(i) = sum(intSH(i,:));
   Sun(i) = sum(intSun(i,:));
end

Answer 1

不幸的是，您尝试尝试的内容可能无法按原样工作。考虑到我知道这个问题的历史，我知道您正在尝试集成一个数组值函数。 This worked in 1d，但我担心它不会在 2d 中工作。

如果您查看 integral2 的帮助，因为它已经在 cmets 中建议，您会看到：

所有输入函数都必须接受数组作为输入并按元素进行操作。函数Z = FUN(X,Y) 必须接受大小相同的数组X 和Y，并返回对应值的数组。

这意味着从fun 输入integral2 的输出不能比输入具有更多的维度；换句话说，integral2 只会为您集成标量函数。

粗略浏览一下选项后，我认为内置的 2d 积分器不支持这一点。您有两种选择，每种都效率低下，因此您应该尝试两种选择，看看哪个更适合您的应用程序。

您已经知道的选项一：循环遍历数组值函数的每个索引，并使用interp2 整合这些标量。这会很慢，因为您需要对数组元素进行循环，并且必须为每个元素调用 interp2d。

选项二是将双积分作为两个单积分来执行。我的意思是正式做

integrated_result = integral(@(t)integral(@(p) fun(t,p),p1,p2,'arrayvalued',true),...
                             t1,t2,'arrayvalued',true);

从p1 到p2 和从t1 到t2 集成p。这会很慢，因为对于外部变量的每个值，您需要调用 integral。