matlab求双重积分的问题

关于这个问题我花了不少时间研究，但很遗憾，后面提供的程序仍然没能彻底解决你的问题。就目前的情况看，我认为自己对你程序想要描述的问题理解上存在偏差的可能性应该不大（但不绝对保证哦），所编的程序从思路上说也没有大问题。

我画了一个示意图，把各个变量的几何含义尽量标在上面，但是，只是到了Q1和Q2的含义还好理解，后面的y1和y2我想不出是什么意思了。

想和你核实两件事：

（1）方便的话，请把要解决问题的原始描述贴出来，以便检验你的程序是否准确反映了原问题，以及我是否误解了你要表达的意思；

（2）a和b是否有取值范围限制？从你目前的描述看，a和b好像没有特别的约束，我按照无约束优化问题来求解，迭代一些步骤后，会出现奇异值导致无法最终完成优化。里面求c用到的函数我试了fzero和fsolve，数值积分试过了多种函数都不行。

参考程序：

function zd

opt=optimset('Display','iter');
X = fminunc(@objfun, [8000 9000],opt)


function Y = objfun(ab)
persistent c
syms X
y=9.616*10^(-16)*X^5 - 5.964*10^(-11)*X^4 + 1.485*10^(-6)*X^3 ...
    - 0.01843*X^2 + 113.0*X - 2.669*10^5;
a = ab(1);
b = ab(2);
ya=subs(y,'X',a);
yb=subs(y,'X',b);
de=sqrt((a-9500)^2+(ya-4000)^2);
dr=sqrt((b-9500)^2+(yb-4000)^2);
dy_dX = diff(y,X);
ds = sqrt(dy_dX^2+1);
J = @(a,b) quadgk( @(x)subs(ds,X,x), a, b);
eq = @(c) abs(J(a,c)-J(c,b)) - abs(de-dr);
opt=optimset('Display','iter');
if isempty(c), c = (a+b)/2; end
c = fsolve(eq,c,opt)
I1 = @(t) J(t,a);
Q1 = @(t) abs(I1(t))+de;
y1 = quadgk(@(t)arrayfun(Q1,t),7650,c);
I2 = @(t) J(b,t);
Q2 = @(t) abs(I2(t))+dr;
y2 = quadgk(@(t)arrayfun(Q1,t),b,14550);
Y=y1+y2;