双重数值积分误差 [英] Double Numerical Integration Error
问题描述
在下面的脚本中，在用于计算双重积分的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
不幸的是，您尝试尝试执行的操作可能无法按原样工作.考虑到我知道问题的历史，因此我知道您正在尝试集成数组值函数. 这在1d中有效，但恐怕它在2d中无效.
如果您在注释中已经建议使用integral2的帮助，则会看到以下内容:
所有输入函数必须接受数组作为输入并按元素进行操作.函数
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.
In the script below, in the for loop for calculating the double integrations, I keep receiving an error, and I'm unsure why:
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);
This script was formed from another which only had to deal with one variable t, so I'm assuming that I've done something wrong in adapting it to two variable functions.
Thanks!
%%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
Unfortunately what you are trying to attempt will probably not work as-is. Considering that I know the history of the question, I know you're trying to integrate an array-valued function. This worked in 1d, but I'm afraid it won't work in 2d.
If you look at the help of integral2 as it has already been suggested in comments, you'll see this:
All input functions must accept arrays as input and operate elementwise. The function
Z = FUN(X,Y)must accept arraysXandYof the same size and return an array of corresponding values.
This means that the output from the fun entering integral2 can't have more dimensions than the input; in other words integral2 will only integrate scalar functions for you.
After a cursory glance over the options, I don't think the built-in 2d integrators support this. You have two options, and each are inefficient, so you should try both and see which is better for your application.
Option one you already know: loop over each index of your array-valued function, and integrate these scalars using interp2. This will be slow because you need a loop over the elements of your array, and interp2d has to be called for each.
Option two is to perform the double integral as two single integrals. I mean formally doing
integrated_result = integral(@(t)integral(@(p) fun(t,p),p1,p2,'arrayvalued',true),...
t1,t2,'arrayvalued',true);
to integrate in p from p1 to p2 and from t1 to t2. This will be slow because for each value of the outer variables you need to call integral.
这篇关于双重数值积分误差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
