python - 使用 SciPy 多次绘制扩散方程的解

Question

我想在不同时间绘制这个方程。那么时间应该是恒定的，x应该变化，然后绘制y？该方程是 transient 扩散方程的解析解。

到目前为止我的代码:

import numpy as np
from scipy.sparse import diags
import scipy as sp
import scipy.sparse
from scipy import special
import matplotlib.pyplot as plt

def Canalytical(intervals, D=1):

    points = 1000

    x=np.linspace(0, 1, intervals+1)
    t=1
    c=np.ones([intervals+1])

    sm = 0
    pos = 0

    for xi in x:

        for i in range(points):

            sm+=sp.special.erfc((1-xi+2*i)/(2*np.sqrt(D*t))) + 
                sp.special.erfc((1+xi+2*i)/(2*np.sqrt(D*t)))

        c[pos] = sm
        pos += 1
        sm = 0

    return c, x

c, xi = Canalytical(intervals=1000) 
plt.plot(xi, c)
plt.show()

Answer 1

图中的方程是错误的。将 x = 0 代入其中，您会发现它不为零。第二个erfc函数前面的符号应该是-。

时间t 应作为参数传递给 Canalytical，因此该函数可用于 t 的多个值。

使用 1000 项求和是过多的，因为 erfc 在无穷大处衰减得非常快。 erfc(10) 大约是 2e-45，远远超出了机器精度，更不用说绘图的分辨率了。

此外，在使用 NumPy 评估函数时，请考虑使用向量化。整个数组 x 可以立即传递给函数，从而消除循环。剩下的就是:

import numpy as np
from scipy import special
import matplotlib.pyplot as plt
def Canalytical(intervals, t=1, D=1):
    points = 1000
    x = np.linspace(0, 1, intervals+1)
    c = np.zeros_like(x)
    for i in range(points):
        c += special.erfc((1-x+2*i)/(2*np.sqrt(D*t))) - special.erfc((1+x+2*i)/(2*np.sqrt(D*t)))
    return x, c
plt.plot(*Canalytical(intervals=1000, t=1)) 
plt.plot(*Canalytical(intervals=1000, t=0.1)) 
plt.plot(*Canalytical(intervals=1000, t=0.01)) 
plt.plot(*Canalytical(intervals=1000, t=0.001)) 
plt.show()

与输出