数值积分误差

Question

当我运行下面的脚本时，我得到错误 Too many input arguments （我发现这在第二个 for 中计算 intNH(5,2) 时特别发生环形）。基本上发生的事情是 i=5 的 fun2 等于 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

等等。

数值积分误差

Numerical integration error

matlab

numerical-integration