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

8楼: Originally posted by 艳雪加油 at 2017-01-13 17:03:41
程序编的太好了，非常感谢！！只不过还有个小小的问题，要求最后是求角的。就是求出的ht要用到求角的公式里（一个公式），大神帮人帮到底，能不能把最后求角公式的matlab程序加到原先的程序里，能运行出来真是太谢 ...

function question_21_new
clear;clc
%2017-1-13
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;
[ht,fval,exitflag,output] =fsolve(@(ht) quadl(@fun1,hs,ht)-Re,30);
if exitflag>0
    angle=ns*(a+hs)*cos(sita)*quadl(@fun2,hs,ht);
    fprintf('\n\n\t\tht = %4.5f\n\n\n\t\tAngle = %4.5f\n\n\n',ht,angle);
end
    function y=fun1(h)
        n=(h>=hs & h<=hs+1).*(1+(Ns+G*(h-hs))*1e-6)  +  (h>hs+1&h<=9).*(1+(N1*exp(-C1*(h-hs-1)))*1e-6)  +  (h>9&h<60).*(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
    function y=fun2(h)
        n=(h>=hs & h<=hs+1).*(1+(Ns+G*(h-hs))*1e-6)  +  (h>hs+1&h<=9).*(1+(N1*exp(-C1*(h-hs-1)))*1e-6)  +  (h>9&h<60).*(1+(N9*exp(-C9*(h-9)))*1e-6);
        y=1./((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的时候初 ...

    if exitflag(i)>0 && (ht(i)>=lb(i)&&ht(i)<=ub(i))

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

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文件，这样就可以调 ...

4楼: Originally posted by 艳雪加油 at 2017-01-13 11:36:31
就是ht运行的在三个区间都有，怎么确定在哪个区间呢？
...

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文件，这样就可以调 ...

7楼: Originally posted by pdl9527 at 2017-01-13 16:12:36
不好意思，还真有个地方错了，其它地方不要改，新建一个m文件，把最开始的代码复制进去，并且把第27行改为：

    if exitflag(i)>0 && (ht(i)>=lb(i)&&ht(i)<=ub(i))
...

9楼: Originally posted by pdl9527 at 2017-01-13 18:17:05
简化了一下程序，以前的代码太冗长了，可读性太差：

function question_21_new
clear;clc
%2017-1-13
a=6370;
Ns=320;
G=-39;
C1=0.1258;
C9=0.1404;
N1=281;
N9=104.6713;
hs=0.15;
Re=50;
sita=0 ...