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

Question

我想绘制不同时间的方程式。所以时间应该是恒定的，x 应该变化然后绘制 y？这个方程是瞬态扩散方程的解析解。

到目前为止我的代码：

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()

与输出

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

Plotting the solution of diffusion equation for multiple times using SciPy

python

scipy

differential-equations