实验十四:水塔水流量估计模型
练习一
1.海水温度随着深度的变化而变化,海面温度较高,随着深度的增加,海水温度越来越低.通过验观测得一组海水温度t与深度h的数据如下:
|   h/m  |   0  |   1.5  |   2.5  |   4.6  |   8.2  |   12.5  |   16.5  |   26.5  | 
|   t/℃  |   23.5  |   22.9  |   20.1  |   19.1  |   15.4  |   11.5  |   9.5  |   8.2  | 
要求:
(1)分别用多种数据插值方法找出温度t与深度h之间的近似函数关系;
(2)找出温度变化最快的深度位置,通过查询相关资料,了解这个特殊位置的实际应用价值.
(1)
clc;clear;
format long
h=[0,1.5,2.5,4.6,8.2,12.5,16.5,26.5];
t=[23.5,22.9,20.1,19.1,15.4,11.5,9.5,8.2];
x=0:0.1:26.5;
%拉格朗日插值
y=lglrcz(h,t,x);
figure(1)
plot(h,t,'.',x,y);
x0=ones(length(h),1);xx=[x0];
for i=1:length(t)-1
    x0=x0.*h';
    xx=[xx,x0];
end
p=inv(xx)*t'%多项式系数(从低到高排列)
%三次样条插值
yy=interp1(h,t,x,'spline');
figure(2)
plot(x,yy,h,t,'.');
pp=spline(h,t)
function y=lglrcz(x0,y0,x)
n=length(x0);
m=length(x);
for i=1:m
    z=x(i);
    s=0.0;
    for k=1:n
        p=1.0;
        for j=1:n
            if j~=k
                p=p*(z-x0(j))/(x0(k)-x0(j));
            end
        end
        s=p*y0(k)+s;
    end
    y(i)=s;
end
end 


p =
23.500000000000000
5.519770260409278
-6.784280947302284
2.409485615066580
-0.388966596330604
0.030972950891118
-0.001174710590966
0.000016740427015
pp =
包含以下字段的 struct:
form: 'pp'
breaks: [0 1.500000000000000 … ]
coefs: [7×4 double]
pieces: 7
order: 4
dim: 1
>> pp.coefs%三次样条插值的各区间的系数如下
ans =
0.666350830669392 -3.625403322677568 3.538815615010219 23.500000000000000
0.666350830669392 -0.626824584665302 -2.839526246004088 22.899999999999999
-0.286563754712924 1.372227907342874 -2.094122923326514 20.100000000000001
0.050422214507364 -0.433123747348549 -0.122004187338434 19.100000000000001
-0.005736766041512 0.111436169330978 -1.280079468201690 15.400000000000000
-0.000611295331494 0.037431887395469 -0.639946824277967 11.500000000000000
-0.000611295331494 0.030096343417537 -0.369833901025942 9.500000000000000
(2)通过观察图像得知,变化最快的点应该在[0,8.2]内。
由于拉格朗日插值的龙格现象,我们不妨把x缩小范围,在小范围内求解。
x1=0:0.1:16.5;%由于龙格现象,我们把x缩小范围进行求解
%拉格朗日插值
y=lglrcz(h,t,x1);
figure(1)
plot(h,t,'.',x1,y);
x0=ones(length(h),1);xx=[x0];
for i=1:length(t)-1
    x0=x0.*h';
    xx=[xx,x0];
end
p=inv(xx)*t'%多项式系数(从低到高排列)
function y=lglrcz(x0,y0,x)
n=length(x0);
m=length(x);
for i=1:m
    z=x(i);
    s=0.0;
    for k=1:n
        p=1.0;
        for j=1:n
            if j~=k
                p=p*(z-x0(j))/(x0(k)-x0(j));
            end
        end
        s=p*y0(k)+s;
    end
    y(i)=s;
end
end
 

pd=polyder(p);
xl=abs(polyval(p,x1));
x1(find(xl==max(xl)))
 
ans =
16.5000
感觉龙格现象的影响,拉格朗日插值不太适合此题。
对于三次样条插值,我们在前四个区间内求解。
syms t
f1=0.6664*(t)^3-3.6254*(t)^2 +3.5388*(t) + 23.5;
f2=0.6664*(t-1.5)^3 -0.6268*(t-1.5)^2  -2.8395*(t-1.5) + 22.9;
f3=-0.2866*(t-2.5)^3+1.3722*(t-2.5)^2 -2.0941*(t-2.5)+ 20.1;
f4=0.0504*(t-4.6)^3 -0.4331*(t-4.6)^2  -0.1220*(t-4.6)+ 19.1;
f11=matlabFunction(diff(f1));
f22=matlabFunction(diff(f2));
f33=matlabFunction(diff(f3));
f44=matlabFunction(diff(f4));
t1=0:0.1:1.5;t2=1.5:0.1:2.5;t3=2.5:0.1:4.6;t4=4.6:0.1:8.2;
f111=f11(t0)
f222=f22(t0);
f333=f33(t0);
f444=f44(t0);
max(abs(f111))
max(abs(f222))
max(abs(f333))
max(abs(f444))
%显然最大值在第一个区间内
t1(find(f111==3.5388))
 
ans =
3.5388
ans =
3.0357
ans =
2.0941
ans =
1.3624
ans =
0
这两种方法得到的结果差别很大,估计是由于选择的插值方法不是很合适,个人疑问:这个数据所给的为何会插值得到奇怪的曲线,刚开始为何深度变大而温度升高呢?
个人建议还是改成拟合较好。
温度变化快的地方可以用来海水温差发电,且发电效率更高。
2.表14.8给出了在低潮时某一平面区域内若干点(x,y)处的水深z值(单位:ft).已知船的吃水深度为5ft.试画出海底的地貌图,并在平面矩形区域(80,196)x(-70,145)内标注哪些地方船要避免进人·
表14.8 水域坐标数据
|   x  |   129.0  |   140.5  |   103.5  |   88.0  |   185.5  |   195.0  |   105.5  |   157.5  |   107.5  |   77.0  |   81.0  |   162.0  |   162.0  |   117.5  | 
|   y  |   7.5  |   141.5  |   23.0  |   147.0  |   22.5  |   137.5  |   85.5  |   -6.5  |   -81.0  |   3.0  |   56.5  |   -66.5  |   84.0  |   -33.5  | 
|   z  |   4  |   8  |   6  |   8  |   6  |   8  |   8  |   9  |   9  |   8  |   8  |   9  |   4  |   9  | 
clc;clear;
x=[129.0,140.5,103.5,88.0,185.5,195.0,105.5,157.5,107.5,77.0,81.0,162.0,162.0,117.5];
y=[7.5,141.5,23.0,147.0,22.5,137.5,85.5,-6.5,-81.0,3.0,56.5,-66.5,84.0,-33.5];
z=[4,8,6,8,6,8,8,9,9,8,8,9,4,9];
xx=77:196;
yy=-81:145;
[xxx,yyy]=meshgrid(xx,yy);
zzz=griddata(x,y,z,xxx,yyy,'cubic');
figure(1)
mesh(xxx,yyy,zzz);
figure(2)
contourf(xxx,yyy,zzz);%等高线
figure(3)
rectangle('Position',[80,-70,116,215]);%画矩形,[左下横,左下纵,长,宽]
for xi=80:196
    for yi=-70:145
        if zzz(find(yy==yi),find(xx==xi))>5
            hold on
            plot(xi,yi,'r.');
        end
    end
end
 



3.估计煤矿的储量.表14.9给出了某露天煤矿在平面矩形区域(1100mx700m)内,纵横均匀的网格交点处测得的煤层厚度(单位:m).由于客观原因,有些点无法测量煤层厚度,用‘—’标出,其中每一网格均为100mx100m的小矩形,试根据这些数据,用不同的方法估算该矩形区域煤矿的储藏量(体积)。
表14.9 煤层厚度
|   A  |   B  |   C  |   D  |   E  |   F  |   G  |   H  |   I  |   J  |   K  | |
|   1  |   —  |   —  |   12.5  |   13.5  |   17.2  |   —  |   8.8  |   14.7  |   8.0  |   13.0  |   —  | 
|   2  |   —  |   —  |   —  |   15.6  |   18.2  |   13  |   6.4  |   8.9  |   9.2  |   11.7  |   —  | 
|   3  |   —  |   12  |   13.5  |   13.5  |   17.8  |   16.9  |   13.2  |   —  |   —  |   —  |   —  | 
|   4  |   7.5  |   12.6  |   14.9  |   18.7  |   17.7  |   17.5  |   14.7  |   13  |   —  |   —  |   6.5  | 
|   5  |   8.9  |   7.8  |   12.4  |   13.5  |   15.7  |   17.6  |   11.7  |   9.6  |   9.2  |   9.5  |   8.6  | 
|   6  |   —  |   —  |   —  |   13.7  |   13.6  |   16.5  |   12.5  |   8.7  |   9.7  |   —  |   —  | 
|   7  |   —  |   —  |   8.6  |   11.8  |   12.5  |   11.3  |   13.4  |   —  |   —  |   —  |   —  | 
首先,我对该题有一个疑问:图中只给出了7*11个数据(姑且将‘—’也称为数据),每个数据是网格交点处的煤层厚度,而题中又说该平面区域为1100mx700m,可是表中的数据怎么画出1100mx700m的区域呢(十一个点只有十个线段啊),那该怎么办呢?
姑且每行都拟合让其拥有7*12个数据,把平面区域改为1100mx600m,这样才能有效地估计煤矿的仓储量。
clc;clear;
format long
x=0:100:1100;
y=0:100:600;
z=xlsread("C:\Users\dell\Desktop\煤层厚度.xlsx");
zz=[];
for i=1:7
    l=length(find(isnan(z(i,:))==0));
    ll=find(isnan(z(i,:))==0);
    lll=z(i,:);
  p=polyfit(1:l,lll(ll),2);
  m=polyval(p,1:12);
  if length(find(m<0))>0
      p=polyfit(1:l,lll(ll),1);
      m=polyval(p,1:12);
  end      
  zz=[zz;m];
end
[x,y]=meshgrid(x,y);
x0=0:5:1100;
y0=0:5:600;
[xx0,yy0]=meshgrid(x0,y0);
zz0=interp2(x,y,zz,xx0,yy0,'cubic');
mesh(xx0,yy0,-zz0);
v=0;
for k=1:length(x0)-1
    for i=1:length(y0)-1
        h=(zz0(i,k)+zz0(i+1,k)+zz0(i,k+1)+zz0(i+1,k+1))/4;
        v=v+25*h;
    end
end
v
 

v =
8.597807745430814e+06
本文由作者自创,由于时间原因,难免出现些许错误,还请大家多多指正。创作不易,请大家多多支持。



















