数值积分误差 [英] Numerical integration error
问题描述
当我在下面运行以下脚本时，出现错误Too many input arguments(我发现这是在第二个for循环中计算intNH(5,2)时特别发生的).基本上发生的是fun2对于i=5等于0(因为相关的G元素是0)，并且MATLAB无法在数值上积分等于零的函数(或者这就是我对此的解释)
When I run the following script below, I get the error Too many input arguments (I have figured out that this specifically happens when calculating intNH(5,2) in the second for loop). Basically what is happening is that fun2 for i=5 equals 0 (because the relevant G element is 0) and MATLAB cannot numerically integrate a function equal to zero (or that's how I've interpreted this).
有人对我如何克服这个问题有任何指示吗?我已经尝试了一些if语句，说如果函数等于0然后是intNH = 0，则其他数值都正常进行积分，但是当我尝试执行此操作时，还会发生更多错误！
Does anyone have any pointers about how I can overcome this? I have tried some if statements saying that if the function is equal to 0 then intNH = 0, else numerically integrate as normal, but more errors keep happening when I try this!
%%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')
推荐答案
您的代码有很多问题.您使用syms和手动设置函数值数组所做的所有魔术都是完全不必要的，并且很难分辨出确切的问题出在哪里(如果将这段代码放入函数中，至少我们会知道这一行错误的来源...).
There are many problems with your code. All this magic you're doing with syms and setting function-valued arrays manually is entirely unnecessary, and it's hard to tell where the exact problem is coming from (if you put this code into a function, at least we'd know the line where the error comes from...).
您应该只使用匿名函数，就像在代码中称为zero的函数一样.您所拥有的只是分析定义的函数和数字.但是，您只能将函数(句柄)存储在cell对象中，而不能存储在数组中.因此，应该为此指定fun{k}和朋友.但是，如果您重写代码以使用数字值(这实际上也是MATLAB的最佳选择)，那么您将成为一个更幸福的人.
You should just use anonymous functions, just like the one called zero in your code. All you have are analytically defined functions and numbers. You can only store function (handles) in cell objects, though, not in arrays. So fun{k} and friends should be assigned for this. But if you rewrite your code to work with numerical values (which is actually optimal for MATLAB too), you'll end up a happier person.
例如，您的第一个循环与此等效:
For instance, your first loop is equivalent to this:
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)的数组.
These are all cells containing array-valued function handles. Each of these functions gets a single value t and outputs an array of length 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
以此类推.
这篇关于数值积分误差的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
