숫자 적분 오류

Question

아래 스크립트를 실행하면 Too many input arguments 오류가 발생합니다.이 문제는 두 번째 for 루프에서 intNH(5,2)을 계산할 때 발생하는 것으로 나타났습니다. 기본적으로 발생하는 것은 i=5에 대한 fun2이 0 (해당 G 요소가 0이기 때문에)이고 MATLAB은 숫자가 0 인 함수를 숫자로 통합 할 수 없습니다 (또는이 방법을 해석 한 방법입니다).

사람이 이것을 극복 할 수있는 방법에 대한 지침이 있습니까? 나는 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 

등을 할 수 있습니다.