数值积分错误

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

等等。