python中任意网格偏导数的有限差分估计

Question

我目前正在用有限差分测试python上的一些代码来求解微分方程，但我正在处理的一些问题和方程如下：

1. 太大了，无法继续手动编写它，而不会出现很大的错误风险。
2. 通常包含非均匀网格，这使得通常的等距网格有限差分不适用。

因此，为了避免错误和避免穷尽的手工选择方程，我决定建立一个函数，估计两个任意网格的偏导数。

守则的作用是：

1. 获取因变量(要被区分的变量)
2. 获得自变量(基本上是离散网格点)
3. 检查要区分的变量的维数。
4. 检查应该执行什么维，例如，对于带有条目I，j，k的3D变量，偏导数可以相对于i、j或k。

然后利用内点正向有限差分法和边界单边差分法，对变量的每一项返回偏导数。

import matplotlib.pyplot as plt
import numpy as np
import sympy as sp
from scipy.optimize import fsolve
from scipy import optimize
import pandas as pd
from numba import njit

@njit
def partial(var,var2,dim):
    y = var
    x = var2
    vardim = var.ndim
    invalid_input = 0
    
    #Forward
    
    if vardim == 1 and dim == 1:
        dydx = np.zeros(len(y))
        
        for i in range(len(dydx)):
            if i == len(dydx) - 1:
                dydx[i] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i] 
                        -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1] 
                        +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2] )
            else:    
                dydx[i] = (y[i+1] - y[i])/(x[i+1] - x[i])
            
    elif vardim == 2 and dim == 1:
        dydx = np.zeros((len(y[:,0]), len(y[0,:])))
        
        for i in range(len(dydx[:,0])):
            for j in range(len(dydx[0,:])):
                if i == len(dydx[:,0]) - 1:
                    dydx[i,j] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i,j] 
                              -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1,j] 
                              +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2,j] )
                else:         
                    dydx[i,j] = (y[i+1,j] - y[i,j])/(x[i+1] - x[i])
    
    elif vardim == 2 and dim == 2:
        dydx = np.zeros((len(y[:,0]), len(y[0,:])))
        
        for i in range(len(dydx[:,0])):
            for j in range(len(dydx[0,:])):
                if j == len(dydx[0,:]) - 1:
                    dydx[i,j] =( (-1/(x[j-2]-x[j]) - 1/(x[j-1]-x[j]))        * y[i,j] 
                              -(x[j-2]-x[j])/((x[j-1]-x[j-2])*(x[j-1]-x[j])) * y[i,j-1] 
                              +(x[j-1]-x[j])/((x[j-1]-x[j-2])*(x[j-2]-x[j])) * y[i,j-2] )                 
                else:    
                    dydx[i,j] = (y[i,j+1] - y[i,j])/(x[j+1] - x[j])
                
    elif vardim == 3 and dim == 1:
        dydx = np.zeros((len(y[:,0,0]), len(y[0,:,0]), len(y[0,0,:])))
        
        for i in range(len(dydx[:,0,0])):
            for j in range(len(dydx[0,:,0])):
                for k in range(len(dydx[0,0,:])):
                    if i == len(dydx[:,0,0]) - 1:
                        dydx[i,j,k] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i,j,k] 
                                    -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1,j,k] 
                                    +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2,j,k] )                 
                    else:
                        dydx[i,j,k] = (y[i+1,j,k] - y[i,j,k])/(x[i+1] - x[i])
                        
    elif vardim == 3 and dim == 2:
        dydx = np.zeros((len(y[:,0,0]), len(y[0,:,0]), len(y[0,0,:])))
        
        for i in range(len(dydx[:,0,0])):
            for j in range(len(dydx[0,:,0])):
                for k in range(len(dydx[0,0,:])):
                    if j == len(dydx[0,:,0])-1:
                        dydx[i,j,k] =( (-1/(x[j-2]-x[j]) - 1/(x[j-1]-x[j]))        * y[i,j,k] 
                                    -(x[j-2]-x[j])/((x[j-1]-x[j-2])*(x[j-1]-x[j])) * y[i,j-1,k] 
                                    +(x[j-1]-x[j])/((x[j-1]-x[j-2])*(x[j-2]-x[j])) * y[i,j-2,k] ) 
                    else:
                        dydx[i,j,k] = (y[i,j+1,k] - y[i,j,k])/(x[j+1] - x[j])
                        
    elif vardim == 3 and dim == 3:
        dydx = np.zeros((len(y[:,0,0]), len(y[0,:,0]), len(y[0,0,:])))
        
        for i in range(len(dydx[:,0,0])):
            for j in range(len(dydx[0,:,0])):
                for k in range(len(dydx[0,0,:])):
                    if k == len(dydx[0,0,:]) - 1:
                        dydx[i,j,k] =( (-1/(x[k-2]-x[k]) - 1/(x[k-1]-x[k]))        * y[i,j,k] 
                                    -(x[k-2]-x[k])/((x[k-1]-x[k-2])*(x[k-1]-x[k])) * y[i,j,k-1] 
                                    +(x[k-1]-x[k])/((x[k-1]-x[k-2])*(x[k-2]-x[k])) * y[i,j,k-2] )    
                    else:
                        dydx[i,j,k] = (y[i,j,k+1] - y[i,j,k])/(x[k+1] - x[k])
    else:
        invalid_input = 1
                        
    if invalid_input == 1:
        raise ValueError('Invalid Input')                
    else:
        return dydx

我已经在使用njit来加速这个函数，它确实比没有它快得多，但即便如此，这段代码仍然比手工写出方程慢得多，并且增加了解决问题的瓶颈。由于我所处理的系统将包含数千个变量，所以时间上的差异变得非常显著。

有什么方法可以提高上述代码的效率吗？

谢谢。

编辑:时间测试示例

compile_function = partial(np.linspace(0,5,6),np.linspace(0,5,6),1)

Test_variable = np.linspace(0,10000000,10000001)
Test_grid = np.linspace(0,10000000,10000001)
t0 = time.time()
partial(Test_variable,Test_grid,1)
elapsed = time.time() - t0
print(elapsed)
0.04687905311584473

Answer 1

正如@J_H所说。在其当前状态下，这是不正确的，而是作为对代码和标准的回顾。

将代码拆分为简单函数

与其将您的方法都放在同一个partial函数中的块中，不如将代码划分为子函数：

def calc_diff_1d_dx(x, y)
    dydx = np.zeros(len(y))
        
    for i in range(len(dydx)):
        if i == len(dydx) - 1:
           dydx[i] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i] 
                    -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1] 
                    +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2] )
        else:    
            dydx[i] = (y[i+1] - y[i])/(x[i+1] - x[i])

然后，您的大select语句变成：

    if vardim == 1 and dim == 1:
        calc_diff_1d_dx(x, y)           
    elif vardim == 2 and dim == 1:
        calc_diff_2d_dx(x, y)
    elif vardim == 2 and dim == 2:
        calc_diff_2d_dy(x, y)
    ...

这意味着我们可以专注于优化/计时，每一个都是独立的。

测试

然而，在您开始优化之前，您需要进行良好的测试，特别是回归测试，以确保您的更改不会以可疑的方式影响您的答案。许多优化可以(而且将会)改变你的答案，但是你必须意识到在你的环境中什么是可以接受的变化。1e-10的绝对差异会影响您的结果吗？1%怎么样？等。

如果没有测试和完全工作的代码，就不可能明智地进行优化。

测试还为我们提供了与时间进行比较的稳定状态的好处。

明显的效率增益

现在有一件显而易见的事情，就是在循环中使用ifs来检查是否正在进行最后一次迭代。尽早停止循环，并作为最后一步对if中的部分进行评估。

       for i in range(len(dydx)):
            if i == len(dydx) - 1:
                dydx[i] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i] 
                        -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1] 
                        +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2] )
            else:    
                dydx[i] = (y[i+1] - y[i])/(x[i+1] - x[i])

化作

       for i in range(len(dydx)-1):
           dydx[i] = (y[i+1] - y[i])/(x[i+1] - x[i])

       i = len(dydx) - 1
       dydx[i] =( (-1/(x[i-2]-x[i]) - 1/(x[i-1]-x[i]))        * y[i] 
                -(x[i-2]-x[i])/((x[i-1]-x[i-2])*(x[i-1]-x[i])) * y[i-1] 
                +(x[i-1]-x[i])/((x[i-1]-x[i-2])*(x[i-2]-x[i])) * y[i-2] )

1. 这使事情变得更清楚了
2. 您并不是在执行循环的每一步的比较和潜在的分支。我不能百分之百地确定Numba的JIT会如何优化它，但是如果它们不够聪明的话，这可能会扼杀许多编译语言的性能。

进一步优化

关于优化的一般好技术，我的第一个建议是找到一个库，它可以做您想做的事情。通常，您不希望自己实现高性能的数学算法，除非作为一种学习经验(除非它们确实不存在/不再受支持)。

在您的代码中，您正在使用numpy，但是在Python中循环它们(在优化中应该避免)，这浪费了numpy的潜力，因为numpy已经优化了C中的向量化例程来完成繁重的工作。如果您希望它是快速的，您应该尽可能多地利用这些，而不是依赖于JIT来获得性能。