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.BESSEL_K, 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)
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")
No description has been provided for this image
No description has been provided for this image
In [9]:
dbg = gl.DbGrid.create(nx,dx)
x = Q.simulate()[0]
dbg["simu"]=x
gp.plot(dbg)
No description has been provided for this image
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()
No description has been provided for this image
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()
No description has been provided for this image
In [12]:
round(res[int(N/2),int(N/2)]/uu[0],3)
Out[12]:
1.009