10楼: Originally posted by 艳雪加油 at 2017-01-13 20:09:15
实在太开心啦！运行结果没有警告更没有error！谢谢大神的程序，谢谢！！
...

ht位于第2区间：

        ht=1.64484

function question_21
clear;clc
%2017-1-12
a=6370;
Ns=320;
G=-39;
C1=0.1258;
C9=0.1404;
N1=281;
N9=104.6713;
hs=0.15;
Re=50;
sita=0.0262;
ns=1.0003;
f1=quadl(@fun1,hs,hs+1);
f2=quadl(@fun2,hs+1,9);
%ht的值位于[hs,hs+1]中
[ht(1),fval,exitflag(1),output] =fsolve(@(ht) quadl(@fun1,hs,ht)-Re,hs);
%ht的值位于[hs+1,9]中
[ht(2),fval,exitflag(2),output] =fsolve(@(ht) f1+quadl(@fun2,hs+1,ht)-Re,5);
%ht的值位于[9,60]中
[ht(3),fval,exitflag(3),output] =fsolve(@(ht) f1+f2+quadl(@fun3,9,ht)-Re,30);
%输出正确的结果
lb=[hs,hs+1,9];
ub=[hs+1,9,60];
for i=1:3
    if fval>0 || (ht(i)>=lb(i)&&ht(i)<=ub(i))
        fprintf('ht位于第%d区间：\n\n\tht=%4.5f\n\n',i,ht(i))
    end
end
function y=fun1(h)
n=1+(Ns+G*(h-hs))*1e-6;
y=n.^2.*(a+h)./(sqrt(n.^2.*(a+h).^2-ns.^2*(a+hs).^2*cos(sita).^2));
end
function y=fun2(h)
n=1+(N1*exp(-C1*(h-hs-1)))*1e-6;
y=n.^2.*(a+h)./(sqrt(n.^2.*(a+h).^2-ns.^2*(a+hs).^2*cos(sita).^2));
end
function y=fun3(h)
n=1+(N9*exp(-C9*(h-9)))*1e-6;
y=n.^2.*(a+h)./(sqrt(n.^2.*(a+h).^2-ns.^2*(a+hs).^2*cos(sita).^2));
end
end

2楼: Originally posted by pdl9527 at 2017-01-13 05:45:27
首先ht可能有三种情况，分别求这三种情况下的ht值，如果方程有解且求得的ht也满足假设的话则输出结果，发现当ht在第二区间的时候能够得到正确的结果。

ht位于第2区间：

        ht=1.64484

当然使用fsolve的时候初 ...

3楼: Originally posted by 艳雪加油 at 2017-01-13 11:31:32
首先，先谢谢！！通过运行呢，我发现了点问题:我感觉程序的最后好像多了个end，而且function y=fun1(h)，function y=fun2(h),function y=fun3(h)在这个长长的程序中是不被识别的，可以另建三个m文件，这样就可以调 ...