[求助,计算d]利用三层风速廓线计算零平面位移距离d

暴龙

大家好！
我有三层（2米、4米和10米）风速数据（见本贴最后附上的一部分实测三层风速数据 WS.txt），需要根据Roland B. Stull 的《An Introduction to Boundary Layer Meteorology》(Page 382)9.7.3b式,求零平面位移（zero plane displacement height,d）的数值解。

下垫面是半退化草地，植被高度h只有约5cm，按照一般的估算d=2/3h,因此d应该只有 约0.03m，我算出来的结果明显不对,不知道哪里理解错或者算错了？

以下是我的matlab程序，利用牛顿迭代法求解:

%++++++++++++++++++++++++++++
a=importdata('F:\WS.txt');

u1=a(:,1);%2m风速
u2=a(:,2);%4m风速
u3=a(:,3);%10m风速

%测量高度
z1=2;
z2=4;
z3=10;

d_Matrix=[];
for i=1:500

    %牛顿迭代法求解
    syms d;
    f=(  u1(i)-u2(i)  )./(  u1(i)-u3(i)  )-  (log(z1-d)-log(z2-d))./ (  log(z1-d)-log(z3-d));% 函数：Stull书本里9.3.7b式
    df=diff(f,d,1);%导数
   
    tol=0.001;%截断误差
    nmax=1000;%最大迭代次数
    x0=1; % 迭代初始猜值


    f = inline(f);
    df = inline(df);
    x(1) = x0 - (f(x0)/df(x0));
    ex(1) = abs(x(1)-x0);
    k = 2;
    while (ex(k-1) >= tol) && (k <= nmax)
        x(k) = x(k-1) - (f(x(k-1))/df(x(k-1)));
        ex(k) = abs(x(k)-x(k-1));
        k = k+1;
        df(x(k-1));
    end
    x=x(end);
    d_Matrix=[d_Matrix ;x];
end

d_Matrix=real(d_Matrix);%求得各时次的零平面位移距离d

%去掉异常值
d_Matrix(d_Matrix<0)=nan;
d_Matrix(d_Matrix>2)=nan;
%++++++++++++++++++++++++++++


请各位大神帮忙看看哪里出问题了，谢谢！


%——————————————————————————————————————
WS.txt：
（初次发帖，不知道怎么设置下载附件不需要扣5个金币，于是把一部分风速数据和我算的零平面位移距离d粘贴在下面：）

2米风速（m/s）    4米风速（m/s）    10米风速(m/s)    零平面位移距离d(m)
3.934    4.625    5.08    1.660561262
5.749    6.798    7.246    1.912416638
4.825    5.689    6.145    1.822537656
4.348    5.128    5.497    1.875262885
3.496    4.146    4.476    1.842958056
2.653    3.097    3.189    NaN
2.212    2.596    2.749    1.932136933
2.04    2.349    2.381    NaN
1.611    1.885    1.772    NaN
1.59    1.906    2.024    1.947190606
1.336    1.511    1.363    NaN
0.945    0.954    0.709    NaN
1.154    1.289    1.267    NaN
1.198    1.364    1.413    NaN
0.993    1.096    0.809    NaN
1.363    1.538    1.327    NaN
1.426    1.559    1.326    NaN
1.842    2.074    1.974    NaN
2.006    2.26    2.252    NaN
3.303    3.774    3.952    1.944632736
3.774    4.204    4.369    1.941238591
4.39    4.909    4.995    NaN
4.298    4.804    5.022    1.909634063
4.31    4.82    3.945    NaN
4.587    5.168    5.529    1.711627008
3.859    4.361    4.542    1.954132924
4.497    5.103    5.288    NaN
5.685    6.437    6.691    1.96510959
6.401    7.257    7.463    NaN