Ошибка численного интегрирования
Когда я запускаю следующий скрипт ниже, я получаю сообщение об ошибке Too many input arguments (Я понял, что это конкретно происходит при расчете intNH(5,2) во-вторых for петля). В основном происходит то, что fun2 за 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')
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
и так далее.