-2

我正在使用 C++ 求解拉普拉斯方程。我正在求解压力场 P。在我的代码中，您可以在下面看到，P 是 x 和 y 的函数。我将其声明为长度为 nx 和 ny 的二维数组。

我得到的最终矩阵结果 P 是正确的。然而，当我尝试使用“printf”打印这个矩阵 P 时，即使数字是完美的，我也得到了这样的错误：

*** stack smashing detected ***
/bin/bash: line 1: 11738 Aborted

下面是我的程序：

#include <cmath>
#include <stdio.h>
const int nx = 5; // number of elements in x-direction
const int ny = 5; // number of elements in x-direction
const int niter = 100; //# of iterations 

int main(){
    double dx = 2/double((nx-1)); 
    double dy = 1/double((ny-1));
    int xmax = 2; int xmin = 0;
    int ymax = 1; int ymin = 0;
    double p[nx][ny];
    double pn[nx][ny];
    double x[nx];
    double y[ny];

    //populate p with zeros     
    for (int xi = 1; xi <=nx;xi++){
        for (int yi = 1;yi<=ny;yi++){
            p[xi][yi] = 0;          
        }
    }

    //populate x and y

    //X
    for (int xnum = 1; xnum <=nx; xnum++){
        x[1] = 0;
        x[xnum+1] = x[xnum] + dx;
    }

    //Y
    for (int ynum = 1; ynum <=ny; ynum++){
        y[1] = 0;
        y[ynum+1] = y[ynum] + dy;
    }

    //initial condition

    for (int yrange = 1; yrange<=ny;yrange++){
        p[nx][yrange] = y[yrange];
    }

    //SOLVING FOR P     

    for (int iter = 1; iter<=niter; iter++){
        //copy values
        for (int xiter= 1; xiter<=nx;xiter++){
            for (int yiter = 1; yiter<=ny;yiter++){
                pn[xiter][yiter]=p[xiter][yiter];
            }
        }
        //main loop
        for (int i = 2; i<=nx-1;i++){
            for (int j = 2; j<=ny-1;j++){
                p[i][j] = ((pow(dy,2)*(pn[i+1][j]+pn[i-1][j]))+(pow(dx,2)*(p[i][j+1]+pn[i][j-1])))/(2*(pow(dx,2)+pow(dy,2)));
            }
        }

        for (int xrange = 2; xrange<=nx-1;xrange++){
            p[xrange][1] = p[xrange][2];
            p[xrange][ny] = p[xrange][ny-1];
        }

    }

    //Testing matrix

        for (int x = 1; x<=nx;x++){
            for (int y =1; y<=ny;y++){
            printf("%1.3f\t",p[x][y]);
            }
        printf("\n");
        }

return 0;
}

一开始我以为是迭代次数太高的问题。换句话说，硝酸盐太高了，所以我决定降低它。不过，我得到了同样的错误。你能给我帮助吗？矩阵 P 的输出结果是正确的，所以我不认为这是求解 P 时的逻辑。非常感谢任何帮助，谢谢！

4

1 回答 1

2

你的循环一直到nx这是未定义的行为，因为你正在访问越界。索引从 0 开始，一直到i < nx.

这个：

for (int xi = 1; xi <= nx; xi++){
    for (int yi = 1; yi <= ny; yi++){
        p[xi][yi] = 0;          
    }
}

应该

for (int xi = 0; xi < nx; xi++){
    for (int yi = 0; yi < ny; yi++){
        p[xi][yi] = 0;          
    }
}

对于所有其他循环也是如此。您的程序可能看起来运行良好，并且可能在与 UB 原因不明显相关的某个地方崩溃。当你有未定义的行为时，任何事情都可能发生（包括什么都没有，所以实际上你可以很高兴你遇到了崩溃并有机会修复它）。

PS：如果你真的想要（不知道为什么，也许你是 FORTRAN 的忠实拥护者），你可以从 1 开始你的循环。但是，我不建议这样做，因为它会使你的循环难以阅读并且容易出错。例如，以下循环中的错误很容易被忽略：

    //copy values (!!! WRONG CODE !!! spot the error to win a FORTRAN fanshirt)
    for (int xiter= 1; xiter < nx+1; xiter++){
        for (int yiter = 1; yiter < ny+1; yiter++){
            pn[xiter-1][yiter-1] = p[xiter][yiter-1];
        }
    }

...顺便说一句，我无意以任何方式抨击 FORTRAN。事实上，我喜欢 FORTRAN，但它是一种不同的语言（实际上是我知道的唯一一种从 1 开始索引的语言）。

于 2017-05-16T07:42:56.377 回答