c++ - 如何从点源模拟二维球面波？

Question

我正在尝试通过使用欧拉积分对波动方程进行数值积分来模拟波(直到我解决了问题，然后我将切换到 runge-kutta)。我正在使用一组 float 作为网格。然后我通过在一个点上改变网格的值来制造干扰。现在，波不是从这一点向所有方向辐射，而是只在一个方向上传播，朝向左上角，即朝向减小的 x 和 y。那么，我的问题是如何让波辐射出去？

这是我的代码

void Wave::dudx(float *input,float *output) //calculate du/dx
{
    for(int y=0;y<this->height;y++)
    {
        for(int x=0;x<this->width;x++)
        {
            output[x+y*this->width]=(this->getPoint((x+1)%this->width,y)-this->getPoint(x,y)); //getPoint returns the value of the grid at (x,y)
        }
    }
}

void Wave::dudy(float *input,float *output) //calculate du/dy
{
    for(int x=0;x<this->width;x++)
    {
        for(int y=0;y<this->height;y++)
        {
            output[x+y*this->width]=(this->getPoint(x,(y+1)%this->height)-this->getPoint(x,y));
        }
    }
}

void Wave::simulate(float dt)
{
    float c=6.0f;

    //calculate the spatial derivatives
    this->dudx(this->points,this->buffer);
    this->dudx(this->buffer,this->d2udx2);

    this->dudy(this->points,this->buffer);
    this->dudy(this->buffer,this->d2udy2);

    for(int y=0;y<this->height;y++)
    {
        for(int x=0;x<this->width;x++)
        {
            this->points[x+y*this->width]+=c*c*(this->d2udx2[x+y*this->width]+this->d2udy2[x+y*this->width])*dt*dt; //I know that I can calculate c*c and dt*dt once, but I want to make it clear what I'm doing.
        }
    }
}

Answer 1

只是为了其他人来这里解决同样的问题。通常将拉普拉斯算子转换为规则网格上的有限差分表达式的方法是:

∆u(x,y) -> idx2*[u(x+1,y) + u(x-1,y) - 2*u(x,y)] +
           idy2*[u(x,y+1) + u(x,y-1) - 2*u(x,y)]

其中 idx2 和 idy2 分别是维度 x 和 y 中网格间距的平方反比。在两个维度的网格间距相同的情况下，这简化为:

∆u(x,y) -> igs2*[u(x+1,y) + u(x-1,y) + u(x,y+1) + u(x,y-1) - 4*u(x,y)]

可以通过将乘法系数隐藏在其他系数中来移除乘法系数，例如c，通过改变它们的测量单位:

∆u(x,y) -> u(x+1,y) + u(x-1,y) + u(x,y+1) + u(x,y-1) - 4*u(x,y)

顺便说一下，因为球体是 3D 物体，所以不可能有 2D 球面波。二维波称为圆波。