GeophyAI

Python与地球物理数据处理

0%

pyseis--有限差分系数计算

发表于 分类于 Valine：
有限差分法中系数计算，基于python的numpy包实现，包括了常规网格二阶导数和交错网格一阶导数的任意偶数阶精度差分系数求解。

正演模拟中二阶导数差分系数计算

整数网格点二阶导数的任意偶数阶精度有限差分系数计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
from numpy.linalg import solve
import numpy as np

def Normal_Grid_Coefficients(N, diff_order = 2):
    """
        Coefficients of arbitrary even order Taylor
        expansion when the discrete values are at 
        integral grid point. Such as for normal grid
        second-order acoustic forward modeling.
        
        Parameters:
        ----------
        N : Half-order of Taylor expansion or the
            length of unilateral operator.
        diff_order : Not used in this method. Just
            for notice.
        
        Returns:
        ----------
        out : 1-D ndarray.
        The length of output is (2*N+1) with the 
        first coefficient is C0.
    """
    
    holder = []
    for i in range(N):
        matrix = []
        for j in range(N):
            matrix.append(np.power((i+1)**j, 2))
        holder.append(matrix)
    holder = np.array(holder).T
    constant = np.zeros(N)
    constant[0] = 1.0
    x = solve(holder, constant)
    
    t1 = 0
    for i in range(N):
        x[i] = x[i]/np.power(i+1, 2)
        t1+=x[i]
    
    coes = np.zeros(N+1)
    coes[0] = -2*t1
    coes[1:len(coes)] = x
    return coes

半网格点一阶导数的任意偶数阶精度有限差分系数计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def Staggered_Grid_Coefficients(N, diff_order = 1):
    """
        Coefficients of arbitrary even order Taylor
        expansion when the discrete values are at 
        half grid point. Such as for staggered grid
        first-order velocity-stress elastic forward
        modeling.
        
        Parameters:
        ----------
        N : Half-order of Taylor expansion or the
            length of unilateral operator.
        diff_order : Differential order of discretization,
            not used in this method. Just for notice.
        
        Returns:
        ----------
        out : 1-D ndarray.
        The length of output is (2*N).
    """
    
    holder = []
    for i in range(N):
        matrix = []
        for j in range(N):
            matrix.append(np.power(2*j+1, 2 * (i + 1) - 1))
        holder.append(matrix)
    holder = np.array(holder)
    constant = np.zeros(N)

    constant[0] = 1.0#*f
    x = solve(holder, constant)
    
    coes = np.zeros(N)
    coes[0:len(coes)] = x
    return coes