数值积分误差

Question

当我在下面运行下面的脚本时，我得到了错误Too many input arguments (我已经发现，在计算第二个for循环中的intNH(5,2)时，这是特别发生的)。基本上发生的事情是，fun2 for i=5等于0 (因为相关的G元素是0)，而MATLAB不能对等于零的函数进行数值积分(或者这就是我对此的解释)。

有人对我如何克服这个问题有什么建议吗？我尝试过一些if语句，如果函数等于0，那么intNH = 0，否则按常规进行数值积分，但是当我尝试时会发生更多的错误！

%%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

%%Call in spherical harmonic coefficients, change the 535 figure as more
%%data is added to the spreadsheets
G(:,1) = xlsread('G Coefficients.xls', 'C3:C535');
G(:,2) = xlsread('G Coefficients.xls', 'E3:E535');
G(:,3) = xlsread('G Coefficients.xls', 'H3:H535');
G(:,4) = xlsread('G Coefficients.xls', 'L3:L535');
G(:,5) = xlsread('G Coefficients.xls', 'Q3:Q535');
G(:,6) = xlsread('G Coefficients.xls', 'W3:W535');
G(:,7) = xlsread('G Coefficients.xls', 'AD3:AD535');
G(:,8) = xlsread('G Coefficients.xls', 'AL3:AL535');
G(:,9) = xlsread('G Coefficients.xls', 'AU3:AU535');

%%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;

syms t %t short for theta here
V = a + (b*cos(t)^2) + (c*cos(t)^4);

zero = @(t) 0*t;

%%Set B and A functions within separate matrices

for i = 1:length(G)

    B(i,1) = G(i,1)*cos(t);

    B(i,2) = 0.5*G(i,2)*(3*(cos(t)^2)-1);

    B(i,3) = 0.5*G(i,3)*(5*(cos(t)^3)-3*cos(t));

    B(i,4) = 1/8*G(i,4)*(35*(cos(t)^4)-30*(cos(t)^2)+3);

    B(i,5) = 1/8*G(i,5)*(63*(cos(t)^5)-70*(cos(t)^3)+15*cos(t));

    B(i,6) = 1/16*G(i,6)*(231*(cos(t)^6)-315*(cos(t)^4)+105*(cos(t)^2)-5);

    B(i,7) = 1/16*G(i,7)*(429*(cos(t)^7)-693*(cos(t)^5)+315*(cos(t)^3)-35*cos(t));

    B(i,8) = 1/128*G(i,8)*(6435*(cos(t)^8)-12012*(cos(t)^6)+6930*(cos(t)^4)-1260*(cos(t)^2)+35);

    B(i,9) = 1/128*G(i,9)*(12155*(cos(t)^9)-25740*(cos(t)^7)+18018*(cos(t)^5)-4620*(cos(t)^3)+315*cos(t));


    A(i,1) = 0.5*R*G(i,1)*sin(t);

    A(i,2) = 0.5*R*G(i,2)*csc(t)*(cos(t)-(cos(t)^3));

    A(i,3) = 0.5*R*G(i,3)*csc(t)*(1.5*(cos(t)^2)-1.25*(cos(t)^4)-0.25);

    A(i,4) = 1/8*R*G(i,4)*csc(t)*(-7*(cos(t)^5)+10*(cos(t)^3)-3*cos(t));

    A(i,5) = 1/8*R*G(i,5)*csc(t)*(-21/2*(cos(t)^6)+35/2*(cos(t)^4)-15/2*(cos(t)^2)+0.5);

    A(i,6) = 1/16*R*G(i,6)*csc(t)*(-33*(cos(t)^7)+63*(cos(t)^5)-35*(cos(t)^3)+5*cos(t));

    A(i,7) = 1/16*R*G(i,7)*csc(t)*(-429/8*(cos(t)^8)+231/2*(cos(t)^6)-315/4*(cos(t)^4)+35/2*(cos(t)^2)-5/8);

    A(i,8) = 1/128*R*G(i,8)*csc(t)*(-715*(cos(t)^9)+1716*(cos(t)^7)-1386*(cos(t)^5)+420*(cos(t)^3)-35*cos(t));

    A(i,9) = 1/128*R*G(i,9)*csc(t)*(-2431/2*(cos(t)^10)+6435/2*(cos(t)^8)-3003*(cos(t)^6)+1155*(cos(t)^4)-315/2*(cos(t)^2)+7/2);

end
    %Turn vector V from a sym to a function and compute V at equator
    Vfunc = matlabFunction(V);
    Veq = Vfunc(0);

    for i=1:length(G)

        fun(1) = 0; fun(2) = 0; fun(3) = 0; fun(4) = 0; fun(5) = 0; fun(6) = 0; 
        fun(7) = 0; fun(8) = 0; fun(9) = 0;

        %Create functions for integration
        fun(1) = matlabFunction(B(i,1).*A(i,1).*(V-Veq).*sin(t));
        fun(2) = matlabFunction(B(i,2).*A(i,2).*(V-Veq).*sin(t));
        fun(3) = matlabFunction(B(i,3).*A(i,3).*(V-Veq).*sin(t));
        fun(4) = matlabFunction(B(i,4).*A(i,4).*(V-Veq).*sin(t));
        fun(5) = matlabFunction(B(i,5).*A(i,5).*(V-Veq).*sin(t));
        fun(6) = matlabFunction(B(i,6).*A(i,6).*(V-Veq).*sin(t));
        fun(7) = matlabFunction(B(i,7).*A(i,7).*(V-Veq).*sin(t));
        fun(8) = matlabFunction(B(i,8).*A(i,8).*(V-Veq).*sin(t));
        fun(9) = matlabFunction(B(i,9).*A(i,9).*(V-Veq).*sin(t));

        %Compute numerical integral for each order and store values

        %Northern Hemisphere
        intNH(i,1) = integral(fun1,0,pi/2);
        intNH(i,2) = integral(fun2,0,pi/2);
        intNH(i,3) = integral(fun3,0,pi/2);
        intNH(i,4) = integral(fun4,0,pi/2);
        intNH(i,5) = integral(fun5,0,pi/2);
        intNH(i,6) = integral(fun6,0,pi/2);
        intNH(i,7) = integral(fun7,0,pi/2);
        intNH(i,8) = integral(fun8,0,pi/2);
        intNH(i,9) = integral(fun9,0,pi/2);   

        %Southern Hemisphere
        intSH(i,1) = int(fun1,t,pi/2,pi);
        intSH(i,2) = int(fun2,t,pi/2,pi);
        intSH(i,3) = int(fun3,t,pi/2,pi);
        intSH(i,4) = int(fun4,t,pi/2,pi);
        intSH(i,5) = int(fun5,t,pi/2,pi);
        intSH(i,6) = int(fun6,t,pi/2,pi);
        intSH(i,7) = int(fun7,t,pi/2,pi);
        intSH(i,8) = int(fun8,t,pi/2,pi);
        intSH(i,9) = int(fun9,t,pi/2,pi);

        %Whole Sun
        intSun(i,1) = int(fun1,t,0,pi);
        intSun(i,2) = int(fun2,t,0,pi);
        intSun(i,3) = int(fun3,t,0,pi);
        intSun(i,4) = int(fun4,t,0,pi);
        intSun(i,5) = int(fun5,t,0,pi);
        intSun(i,6) = int(fun6,t,0,pi);
        intSun(i,7) = int(fun7,t,0,pi);
        intSun(i,8) = int(fun8,t,0,pi);
        intSun(i,9) = int(fun9,t,0,pi);

        %Sum over all orders to get a value for dH/dt
        NH(i) = sum(dintNH(i,:));
        SH(i) = sum(intSH(i,:));
        Sun(i) = sum(intSun(i,:));

    end

x = linspace(1643,2175,533);

figure(1)
plot(x,NH)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Northern Hemisphere magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')

figure(2)
plot(x,SH)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Southern Hemisphere magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')

figure(3)
plot(x,Sun)
xlabel('Carrington Rotation')
ylabel('Magnetic helicity transport')
title('Whole sun magnetic helicity transport via twisting motions on the boundary vs Carrington Rotations')

Answer 1

您的代码有很多问题。使用syms和手动设置函数值数组所做的所有这些魔术都是完全没有必要的，而且很难判断问题的确切来源(如果将这些代码放入函数中，至少我们知道错误的来源是……)。

您应该使用匿名函数，就像代码中名为zero的函数一样。你所拥有的只是解析定义的函数和数字。但是，您只能将函数(句柄)存储在cell对象中，而不能存储在数组中。因此，fun{k}和朋友应该为此分配。但是，如果您重写代码来处理数值(这对于MATLAB来说也是最优的)，那么您最终将成为一个更快乐的人。

例如，您的第一个循环相当于以下内容：

B{1} = @(t)  G(:,1)*cos(t);

B{2} = @(t)  0.5*G(:,2)*(3*(cos(t)^2)-1);

B{3} = @(t)  0.5*G(:,3)*(5*(cos(t)^3)-3*cos(t));

B{4} = @(t)  1/8*G(:,4)*(35*(cos(t)^4)-30*(cos(t)^2)+3);

B{5} = @(t)  1/8*G(:,5)*(63*(cos(t)^5)-70*(cos(t)^3)+15*cos(t));

B{6} = @(t)  1/16*G(:,6)*(231*(cos(t)^6)-315*(cos(t)^4)+105*(cos(t)^2)-5);

B{7} = @(t)  1/16*G(:,7)*(429*(cos(t)^7)-693*(cos(t)^5)+315*(cos(t)^3)-35*cos(t));

B{8} = @(t)  1/128*G(:,8)*(6435*(cos(t)^8)-12012*(cos(t)^6)+6930*(cos(t)^4)-1260*(cos(t)^2)+35);

B{9} = @(t)  1/128*G(:,9)*(12155*(cos(t)^9)-25740*(cos(t)^7)+18018*(cos(t)^5)-4620*(cos(t)^3)+315*cos(t));


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

A{2} = @(t)  0.5*R*G(:,2)*csc(t)*(cos(t)-(cos(t)^3));

A{3} = @(t)  0.5*R*G(:,3)*csc(t)*(1.5*(cos(t)^2)-1.25*(cos(t)^4)-0.25);

A{4} = @(t)  1/8*R*G(:,4)*csc(t)*(-7*(cos(t)^5)+10*(cos(t)^3)-3*cos(t));

A{5} = @(t)  1/8*R*G(:,5)*csc(t)*(-21/2*(cos(t)^6)+35/2*(cos(t)^4)-15/2*(cos(t)^2)+0.5);

A{6} = @(t)  1/16*R*G(:,6)*csc(t)*(-33*(cos(t)^7)+63*(cos(t)^5)-35*(cos(t)^3)+5*cos(t));

A{7} = @(t)  1/16*R*G(:,7)*csc(t)*(-429/8*(cos(t)^8)+231/2*(cos(t)^6)-315/4*(cos(t)^4)+35/2*(cos(t)^2)-5/8);

A{8} = @(t)  1/128*R*G(:,8)*csc(t)*(-715*(cos(t)^9)+1716*(cos(t)^7)-1386*(cos(t)^5)+420*(cos(t)^3)-35*cos(t));

A{9} = @(t)  1/128*R*G(:,9)*csc(t)*(-2431/2*(cos(t)^10)+6435/2*(cos(t)^8)-3003*(cos(t)^6)+1155*(cos(t)^4)-315/2*(cos(t)^2)+7/2);

这些都是包含数组值函数句柄的单元格。每个函数都得到一个值t，并输出一个长度为length(G)的数组。

稍后，你可以做

V = @(t) a + (b*cos(t)^2) + (c*cos(t)^4);
Veq = V(0);
% todo: pre-allocate fun
intNH = zeros(length(G),9);
for k=1:9
   fun{k} = @(t) B{k}(t).*A{k}(t).*(V(t)-Veq).*sin(t);
   intNH(:,k) = integral(fun{k},0,pi/2,'arrayvalued',true);
end

诸若此类。