实验十四：水塔水流量估计模型
练习一
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
本文由作者自创，由于时间原因，难免出现些许错误，还请大家多多指正。创作不易，请大家多多支持。
扫一扫
2115
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
