怎么在matlab中求解复合分段函数的二重积分？

Question

∫0 1.8∫0.28 1(1-brt(r,t)/art(r,t))2*pi*rdrdt，(1.8,0),(0.28,1)分别为积分上下限
且function a=art(r,t)
if(Ert(r,t)<=2.47233)
a=0;
else if(2.47233<Ert(r,t)&Ert(r,t)<=5.06625)
a=1.01325*0.022*((10^6*Ert(r,t))/1.01325e5-24.4)^2;
else if(5.06625<Ert(r,t)&Ert(r,t)<=12.159)
a=1.01325*0.05*((10^6*Ert(r,t))/1.01325e5-24.4)^1.75;
else if(12.159<Ert(r,t)&Ert(r,t)<=25.33125)
a=17.7*Ert(r,t)-1.01325*72.5;
else
a=1116.6015*exp(-271e-6/(Ert(r,t)/1.01325e4));
end
end
end
end
function brt(r,t)
b=11.26*exp(-4.56/(20.009*(exp(-t/80.5369)-exp(-t/0.7742))/(r*log(825/2.8))))/(20.009*(exp(-t/80.5369)-exp(-t/0.7742))/(r*log(825/2.8)))
art(r,t)，brt(r,t)分别为两自定义函数，求高手指点怎么求这个二重积分。
这是原题说明

Answer 1

建议楼主把原题贴出来（而不是把问题的matlab code贴出来）。追问

这是这个问题的说明，希望能帮解答一下，非常感谢！

追答

所给表达式的二重数值积分为NaN:

function phi = erchong
% 二重积分调用格式：Q = dblquad(integrnd, pi, 2*pi, 0, pi)
phi = -2*pi*dblquad(@calcetalpha,0,1.8,0.28,1.025);

function etalpha = calcetalpha(r,t)
% 定义二重积分中的被积函数.
E = calcE(r,t);
myeta = calceta(E);
myalpha = calcalpha(E);
etalpha = r.*(1-myeta./myalpha);
end
function E = calcE(r,t)
% p = 出现在 E 表达式中的参数;
p = [20.069 80.5369 0.7742 82.5 0.28];
E = p(1)*( exp(-t/p(2))-exp(-t/p(3)) )./r/log(p(4)/p(5));
end
function myeta = calceta(E)
% p = 出现在 eta 表达式中的参数;
p = [11.26 4.56];
myeta = p(1)*exp(-p(2)./E)./E;
end
function myalpha = calcalpha(E)
% p = 区间分段点;
p = [2.47233 5.06625 12.159 25.33125];
% q = 出现在 alpha 表达式中的参数;
q = [1.01325 2.2e-7 24.4 5e-7 17.7 72.5 1116.6015 2.71];
% u = alpha 中子式单独计算，化简 alpha 的表达式;
u = 10*E/q(1) - q(3);
% 分段函数 alpha 的 matlab 表达:
myalpha = 0 .* (Ep(1) & Ep(2) & Ep(3) & Ep(4));
end
end

本回答被提问者采纳