Fipy 中的狄拉克增量源项

Question

我想知道如何将 Dirac delta 函数表示为 Fipy 中的源项。我想解决以下等式 

我试过下面的代码

from fipy import *

nx = 50

ny = 1

dx = dy = 0.025  # grid spacing

L = dx * nx

mesh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)

phi = CellVariable(name="solution variable", mesh=mesh, value=0.)

Gamma=1

delta=1  # I want knowing how to make this right. 

eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+delta

valueTopLeft = 0

valueBottomRight = 1

X, Y = mesh.faceCenters

facesTopLeft = ((mesh.facesLeft & (Y > L / 2)) | (mesh.facesTop & (X < L / 2)))

facesBottomRight = ((mesh.facesRight & (Y < L / 2)) |

                    (mesh.facesBottom & (X > L / 2)))

phi.constrain(valueTopLeft, facesTopLeft)

phi.constrain(valueBottomRight, facesBottomRight)

timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)

steps = 100

results=[]

for step in range(steps):

    eqG.solve(var=phi, dt=timeStepDuration)

    results.append(phi.value)

该代码正在运行，但我想要确切的 Dirac delta 函数。我查找了 numerix 模块，但找不到这样的功能。Sx1 和 Sy1 是常数。我正在使用 python 2.7

Answer 1

像使用扩散接口方法一样平滑 Dirac delta 函数可能是一个好主意（请参见此处的方程 11、12 和 13 ）。所以，这是一种选择def delta_func(x, epsilon):&nbsp; &nbsp; return ((x < epsilon) & (x > -epsilon)) * \&nbsp; &nbsp; &nbsp; &nbsp; (1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilon2 * epsilon是狄拉克三角函数的宽度，选择为几个网格间距宽。您也可以只使用1 / dx并选择离狄拉克三角函数位置最近的网格点。但是，我认为这变得更加依赖网格。这是一维的工作代码。from fipy import *nx = 50dx = dy = 0.025&nbsp; # grid spacingL = dx * nxmesh = Grid1D(dx=dx, nx=nx)phi = CellVariable(name="solution variable", mesh=mesh, value=0.)Gamma=1def delta_func(x, epsilon):&nbsp; &nbsp; return ((x < epsilon) & (x > -epsilon)) * \&nbsp; &nbsp; &nbsp; &nbsp; (1 + numerix.cos(numerix.pi * x / epsilon)) / 2 / epsilonx0 = L / 2.eqG = TransientTerm() == DiffusionTerm(coeff=Gamma)+ delta_func(mesh.x - x0, 2 * dx)valueTopLeft = 0valueBottomRight = 1timeStepDuration = 10 * 0.9 * dx ** 2 / (2 * 0.8)steps = 100viewer = Viewer(phi)for step in range(steps):&nbsp; &nbsp; res = eqG.solve(var=phi, dt=timeStepDuration)&nbsp; &nbsp; print(step)&nbsp; &nbsp; viewer.plot()input('stopped')在这里，epsilon = 2 * dx，任意选择，并且 delta 函数以 为中心L / 2。2D 只需要乘以函数。