SPDE for Markovian Model¶
In [1]:
import gstlearn as gl
import gstlearn.plot as gp
import numpy as np
import matplotlib.pyplot as plt
In [2]:
Markov = True
In [3]:
dx = [.5,.5]
nx = [400,400]
db = gl.DbGrid.create(nx, dx)
In [4]:
mesh = gl.MeshETurbo(nx,dx)
In [5]:
scale = 1
sill = 5.
param = 2.
ranges = [1.8 * scale, 1. * scale]
angles = [30.,0]
if Markov :
model = gl.Model.createFromParam(gl.ECov.MARKOV, ranges=ranges, angles = angles,
sill=sill, flagRange=False)
cova = model.getCova(0)
cova.setMarkovCoeffsBySquaredPolynomials([1,-1.5],[1,-1],0.001)
#coeffs = np.array(cova.getMarkovCoeffs())
#for i in range(len(coeffs)):
# coeffs[i]*= scale**(2*i)
#cova.setMarkovCoeffs(coeffs)
else :
model = gl.Model.createFromParam(gl.ECov.MATERN, param=param,
ranges=ranges, sill=sill,
angles =angles, flagRange=False)
model
Out[5]:
Model characteristics ===================== Space dimension = 2 Number of variable(s) = 1 Number of basic structure(s) = 1 Number of drift function(s) = 0 Number of drift equation(s) = 0 Covariance Part --------------- Markov (Third Parameter = 1) - Sill = 5.000 - Ranges = 12.471 6.928 - Theo. Ranges = 1.800 1.000 - Angles = 30.000 0.000 - Rotation Matrix [, 0] [, 1] [ 0,] 0.866 -0.500 [ 1,] 0.500 0.866 Total Sill = 5.000 Known Mean(s) 0.000
In [6]:
Q = gl.PrecisionOp(mesh,model.getCova(0))
In [7]:
indx = int(nx[0]/2)
indy = int(nx[1]/2)
ind = indy + nx[1] * indx
covy = np.array(Q.evalCov(ind)).reshape(nx)
In [8]:
xx = np.arange(len(covy)/2)*dx[0]
if not Markov :
covM = [model.eval(gl.SpacePoint([0,0]),gl.SpacePoint([i,0])) for i in xx]
plt.plot(xx,covM,"--",label = "Covariance")
plt.plot(xx,covy[indx,indy:(indy+nx[1])],label = "discretized")#/covy[indx,indx])
covy[indx,indx]
ax = plt.legend()
plt.show()
i = plt.imshow(covy.reshape(nx),origin = "lower")
In [9]:
dbg = gl.DbGrid.create(nx,dx)
x = Q.simulateOne()
dbg["simu"]=x
gp.plot(dbg)
In [10]:
N = 2**8
cova = model.getCova(0)
mH1 = 3.5 * np.max(cova.getRanges())
mH2 = 3.5 * np.max(cova.getRanges())
result = np.array(cova.evalCovFFT([mH1,mH2],N).getValues())
res = np.array(result).reshape((N,N))
i = plt.imshow(res,extent = [-mH1,mH1,-mH2,mH2],origin = "lower")
c = plt.colorbar()
In [11]:
X1 = np.linspace(-mH1,mH1,N)
X2 = np.linspace(-mH2,mH2,N)
plt.plot(X1 , res[int(N/2),:], label="FFT")
uu = covy[indx,indy:(indy+100)]
plt.plot(np.arange(len(uu))*dx[0],uu,"--",label = "Discretized covariance")
ax = plt.legend()
In [12]:
round(res[int(N/2),int(N/2)]/uu[0],3)
Out[12]:
1.009
In [ ]: