더블 수치 통합은 아래의 스크립트에서 오류

Question

, 나는 확실 해요 하나의 변수 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 

Answer 1

당신은 아마-그대로 작동하지 않습니다 시도하려고 불행하게도 무엇을. 질문의 역사를 알고 있다고 생각하면 배열 값 함수를 통합하려고합니다. This worked in 1d,하지만 2d에서 작동하지 않을까 걱정됩니다.

이미 의견 제안 된대로 integral2의 도움을 보면

,이 나타납니다 :

모든 입력 기능이 입력으로 배열을 받아 elementwise을 운영해야합니다. Z = FUN(X,Y) 함수는 동일한 크기의 배열 X 및 Y을 받아 들여야하며 해당 값의 배열을 반환해야합니다.

이것은 fun로부터의 출력이 입력 integral2 이상의 치수를 가질 수 있음을 입력하는 수단; 즉 integral2은 스칼라 함수 만 통합합니다.

옵션을 통해 피상적 인 눈 후, 나는 생각하지 않는다 내장 된 2D 통합이를 지원합니다. 두 가지 옵션이 있으며 각각은 비효율적이므로 두 가지 방법을 모두 시도해보고 어떤 것이 더 나은지 확인해야합니다.

옵션 당신이 이미 알고 하나를 배열 값 함수의 각 인덱스 돌이, 그리고 interp2를 사용하여 이러한 스칼라를 통합 할 수 있습니다. 배열 요소에 대해 루프가 필요하기 때문에 느려질 것이며 interp2d을 각각 호출해야합니다.

옵션 2는 두 개의 단일 적분으로 이중 적분을 수행하는 것입니다. 나는 공식적으로 p1에서 p2과 t1에서 t2에 p에 통합 할 수

integrated_result = integral(@(t)integral(@(p) fun(t,p),p1,p2,'arrayvalued',true),... 
          t1,t2,'arrayvalued',true); 

을하는 의미. 외부 변수의 각 값에 대해 integral에 전화해야하므로 속도가 느려집니다.