一维地下水运动MATLAB模拟.doc

Problem 1The ground flow equation for 1D heterogeneous, isotropic porous medium with a constant aquifer thickness is given by:Where h(x,t) is the hydraulic head, T(x) is the transmissivity, and S is the storativity. The boundary conditions imposed are constant head hL (15m) and hR (5m) at the left and right ends of soil column, respectively. The initial condition is 0 or random number at each node. The length of aquifer is 1m (L=1m). The storativity is 1 (S = 1).Develop a MATLAB program which can handle a heterogeneous transmissivity field using the implicit method.1) Test the program for the case of homogenous parameters, and compare the results with the results generated from flow equation with homogenous transmissivity field (i.e., )2) Test the program for patch- heterogeneous T field, which comprises a two-zone T field with a fivefold difference in transmissivity (T1=5T2, T1=0.01m2s-1, T2=0.002m2s-1 ; or, T1=0.2T2, T1=0.002m2s-1, T2=0.01m2s-1) and an interface in the middle of the domain.Problem 1解答1 理论推导非均质一维地下水运动方程为： 如右图，各个点的水头分别为，水头和之间水力传导系数为。将方程(1)进行离散化：设,则 设 则 所以： 设, 则, , 2 编码实现根据第一部分的理论推导编写代码2.1 均质情况下一维地下水运动function plot_Data = oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt)% Finite difference method to solve 1-D flow equation% Writed By Dongdong Kong, 2014-01-07% Sun Yat-Sen University, Guangzhou, China% emal: % -% one dimension groundwater flow model script% which can handle a heterogeneous transmissivity field% -% because of T value, this function only suit for patch-heterogeneous T % field comprises a two-zone T field, or homogeneous field. you can Modify% input discreted T value to suit other heterogeneous situation.% -% hL % left constant hydraulic head (m)% hR % left constant hydraulic head (m)% L % the length of aquifer% S % storativity% T % transmissivity% method% forward or backward approximation% -x=0:dx:L; nx=length(x);w=T/S*dt/dx/dx;h=rand(nx,1);h(1)=hL; h(end)=hR; times = 100;Result = zeros(nx,times+1); if strcmp(method,forward) % forward approximation W=(1-2*w)*diag(ones(nx,1),0)+. w*diag(ones(nx-1,1),1)+. w*diag(ones(nx-1,1),-1); W(1,1)=1;W(1,2)=0; W(nx,nx)=1;W(nx,nx-1)=0; for i=1:times+1 Result(:,i)=h; h=W*h; endelseif strcmp(method,backward) % backward approximation w=T/S*dt/dx/dx; W=(1+2*w)*diag(ones(nx,1),0) - . w*diag(ones(nx-1,1),1) - . w*diag(ones(nx-1,1),-1); W(1,1)=1;W(1,2)=0; W(nx,nx)=1;W(nx,nx-1)=0; for i=1:times+1 Result(:,i)=h; h=Wh; endelse warning(error method input, method can be only set forward,backward)endfigure timeID=0 1 5 10 20 40 100+1;plot_Data=Result(:,timeID);plot(x,plot_Data,linewidth,1.5)legend(t= 0,t= 1,t= 5,t= 10,t= 20,t= 40,t= 100)2.2 非均质情况下一维地下水运动function plot_Data = oneDimGroudwaterFlowHet(hL, hR, L, S, T)% Finite difference method to solve 1-D flow equation% Writed By Dongdong Kong, 2014-01-07% Sun Yat-Sen University, Guangzhou, China% emal: % -% one dimension groundwater flow model script% which can handle a heterogeneous transmissivity field% -% because of T value, this function only suit for patch-heterogeneous T % field comprises a two-zone T field, or homogeneous field. you can Modify% input discreted T value to suit other heterogeneous situation.% -% hL % left constant hydraulic head (m)% hR % left constant hydraulic head (m)% L % the length of aquifer% S % storativity% T % transmissivity% -nx=50;x=linspace(0,L,nx*2+1);dx=x(2)-x(1);dt = 0.1;N = length(x);Ti=repmat(T(1),1,nx),repmat(T(2),1,nx); W=zeros(N);w=dt/dx/dx/S;for i=2:N-1 W(i, i+1) =-w*Ti(i); W(i, i) =1+w*(Ti(i)+Ti(i-1); W(i, i-1)=-w*Ti(i-1);endW(1,1)=1;W(end,end)=1;h=rand(N,1);h(1)=hL; h(end)=hR;times=10000; %simualtion timesResult=zeros(N,times+1); figure,hold onfor i=1:times+1 Result(:,i)=h; h=Wh;endtimeID=0 1 5 10 20 40 100 1000*10+1;plot_Data=Result(:,timeID);plot(x,plot_Data,linewidth,1.5)legend(t= 0,t= 1,t= 5,t= 10,t= 20,t= 40,t= 100,t=1000)2.3 主函数clear,clc% Writed By Dongdong Kong, 2014-12-30% one dimension groundwater flow model script% which can handle a heterogeneous transmissivity fieldhL = 15; % left constant hydraulic head (m)hR = 5; % left constant hydraulic head (m)L = 1; % the length of aquiferS = 1; % storativity% # question1.1dt=1;dx=0.2; %for delta_w 0.5T=0.01; method=forward;Hom_forward = oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt);dx=0.01;T=0.01; method=backward;Hom_backward = oneDimGroudwaterFlowHom(hL,hR,L,S,T,method,dx,dt);T=0.01 0.01;Hom = oneDimGroudwaterFlowHet(hL,hR,L,S,T);% # question1.2T=0.01 0.002;Het_situ1 = oneDimGroudwaterFlowHet(hL,hR,L,S,T);T=0.002 0.01;Het_situ2 = oneDimGroudwaterFlowHet(hL,hR,L,S,T);save(oneDimFlowPlotData.mat,Hom_forward,Hom_backward,Hom,Het_situ1,Het_situ2)3. 结果输出Question1.1 result：由图一可以看出非均质插分和均质情况下向前插分、向后插分，计算结果和收敛速度大致相同，考虑到向前插分中，因此向前差分时取，其他情况取。图一. 非均质与均质求解下的差异，从左到右分别为非均质T插分、均质情况下向后插分、均质情况下向前插分Question1.2 result：T1=0.01m2/s, T2=0.002m2/s 和 T1=0.2T2, T1=0.002m2/s, T2=0.01m2/s非均质T情况下插分结果如图二。图二. 非均质情况下两个情景下插分结果Problem 2Use rainfall from 8 stations to estimate the annual rainfall for 6 stations using the kriging ordinary methods.The locations of the 8 stations where rainfall is available are as follows:No.station namelatitudelogitude1boulder40.033105.2672castle39.383104.8673green9ne39.267104.7504lawson39.767105.6335longmont40.252105.1506manitou38.850104.9337parker39.533104.6508woodland39.100105.083The 15-year rainfall records for 8 the stations are available in the Excel file “Rain_8Stations_Known.xls”The locations of the 6 stations where rainfall need to be estimated are as follows:No.station namelatitudelogitude1byers39.750104.1332estes40.383105.5173golden39.700105.2174green9se39.100104.7335lakegeor38.917105.4836morrison39.650105.200The annual rainfall for 6 the stations are available in the Excel file “AnnualRain_6Stations_ToBeEstimated.xls”. This will be used to calculate the estimate error by using different interpolation methods1) For each of 6 stations (i.e., byers, estes, golden, green9se, lakegeor, and morrison) where the annual will be estimated, estimate the annual rainfall based on rainfall at 8 stations (i.e., boulder, castle, green9ne, lawson, longmont, manitou, parker and woodland) using the kriging method.2) Estimate the annual rainfall using the inverse-distance-square method and the arithmetic mean method, and compare the error from the kriging method.3) Calculate the error variance for kriging method at each of 6 stationsNote: use an exponential function to model covariance: , where a and b are parameters to be determined, and d is the distance between two points.Problem 2解答1. Kriging MATLAB 编码如下KrigineInterpolate.mclear,clc% Writed By Dongdong Kong, 2014-01-08% Sun Yat-Sen University, Guangzhou, China% email: % -dataDir = ./data/; stationInfo_known= xlsread(dataDir,LatitudeLogitude_8Stations_Known.xls);stationInfo_pre = xlsread(dataDir,LatitudeLogitude_6Stations_ToBeEstimated.xls);Rain_known = xlsread(dataDir,Rain_15year_8Stations_Known.xls);Rain_preReal = xlsread(dataDir,AnnualRain_6Stations_ToBeEstimated.xls);stationInfo_known = stationInfo_known(:,3,4);stationInfo_pre = stationInfo_pre(:,3,4);nYear = size(Rain_known,1); Np_know = 8; % station number of know dataNp_pre = 6; % station number to predictRainPredict = zeros(nYear, Np_pre);for k = 1:Np_pre stationInfo = stationInfo_pre(k,:); stationInfo_known; nStation = size(stationInfo,1); dispPoint=zeros(nStation); % calculate two station distance for i=1:nStation dispPoint(:,i) = distance(stationInfo, stationInfo(i,:); end % covariance mat_cor = cov(Rain_known); % simulate Cov function x = dispPoint(2:end,2:end); y = mat_cor; x_tri = tril(x); y_tri = tril(y); X = x_tri(:); ind = X=0; Y = y_tri(:); % delete diag zeros data X_nozeros = X(ind); Y_nozeros = Y(ind); Y_zeros = diag(y); % General model Exp1: % y = a*exp(b*x) fitResult, gof = expFit(X_nozeros, Y_nozeros); a = fitResult.a; b = fitResult.b; C0 = mean(Y_zeros) - b; Cov_fit = reshape(feval(fitResult,x),8,8); for i=1:length(Cov_fit) Cov_fit(i, i) = C0 + b; end C = Cov_fit,ones(8,1); ones(1,8),0; dh = dispPoint(2:end,1); D = feval(fitResult,dh);1; W = CD; %last element is mu w = W(1:end-1); RainPredict(:,k) = Rain_known*w;endRain_kriging = mean(RainPredict)expFit.mfunction fitresult, gof = expFit(X1, Y1)%MATLAB AUTO GENERATE CODE expFit FUNCTIONxData, yData = prepareCurveData( X1, Y1 );% Set up fittype and options.ft = fittype( exp1 );opts = fitoptions( ft );opts.Display = Off;opts.Lower = -Inf -Inf;opts.StartPoint = 871049.657333104 -6.29140824382876;opts.Upper = Inf Inf;% Fit model to data.fitresult, gof = fit( xData, yData, ft, opts );% -计算结果保留在第二问的表格里。2. 反距离权重插值、算术平均值和kriging插值的误差比较反距离权重插值采用R语言下的工具包gstat的idw函数进行，对于每个站点采用8个雨量站进行插值。IDW Interpolation R 代码如下：library(gstat)library(maptools)# Loading required package: sp# Checking rgeos availability: TRUErm(list=ls()setwd(F:/Users/kongdd/Documents/MATLAB/finalWork/krigineInterpolation) stationInfo_known=read.table(./data/LatitudeLogitude_8Stations_Known.txt,sep=t,head=T)stationInfo_pre = read.table(./data/LatitudeLogitude_6Stations_ToBeEstimated.txt,sep=t,head=T)Rain_known = read.table(./data/Rain_15year_8Stations_Known.txt,sep=t,head=T)Rain_preReal = read.table(./data/AnnualRain_6Stations_ToBeEstimated.txt,head=T)Rain_known - t(Rain_kn