라플라스방정식의 전산적 풀이

이런 형태의 편미분 방정식을 라플라스 방정식이라 한다.

이걸 전산적으로 푸는 법을 알아보자.

 

1차원 라플라스방정식

은 대충 이렇다

이걸 양변을 두 번 적분하면 y=ax+b라는 직선이 나온다.

이 2계 상미분방정식은 경계조건(양끝점)을 알면 a와 b를 구할 수 있다.

또한 이 두 점 사이의 중간값도 아주 쉽게 구할 수 있을 것이다.

 

양 끝점의 평균값 말이다.

 

이걸로 알 수 있는 것은 라플라스방정식은 양 끝점을 알면 모든 값을 알 수 있다는 것이다.

 

2차원 라플라스방정식을 보자

이 방정식의 해는 고무막을 아주 조금 위로 당긴 것과 같은 해를 갖는다.(더 많이 당긴 건 막의 진동이라고 베셀 방정식쓰는 거 있다.)

 

그럼 여기서 특정 위치에서의 값은 어떻게 구하냐?

요렇게 구한다.

미적분학 모르는 사람들을 위해 설명하자면 아주 작은 반경의 원을 그려서 그 원 위의 점들이 가지는 값들의 평균값이 구하고자 하는 점의 값이라는 뜻이다.

 

그럼 코드를 보자

#include <iostream>
#include <vector>
#include <tuple>
#include <cassert>

template<typename Data>
class Matrix
{
    int myRows;
    int myColumns;
    std::vector<Data> myData;

    class Proxy
    {
        Matrix *myOwner;
        int myRow;
    public:
        Proxy(Matrix *owner, int row) : myOwner(owner), myRow(row)
        {
        }

        Data &operator[](int column)
        {
            return myOwner->get(myRow, column);
        }
    };

public:
    Matrix(int rows, int columns) : myRows(rows), myColumns(columns), myData(rows * columns)
    {
    }

    Data &get(int row, int column)
    {
        assert(0 <= row && row < myRows && 0 <= column && column < myColumns);
        return myData[row * myColumns + column];
    }

    int getMax_x()
    {
        return myRows;
    }

    int getMax_y()
    {
        return myRows;
    }

    Proxy operator[](int row)
    {
        return Proxy(this, row);
    }
};

struct init
{
    double v;
    double prev_v;
    bool memo;
};

void Laplace2D_solve(Matrix<init> &V, double c)
{
    while (1)
    {
        int check = 0;

        for (int i = 1; i < V.getMax_x() - 1; ++i)
        {
            for (int j = 1; j < V.getMax_y() - 1; ++j)
            {
                if (V[i][j].memo)
                    continue;

                V[i][j].v = (V[i - 1][j].v + V[i + 1][j].v + V[i][j - 1].v + V[i][j + 1].v) / 4.0;
            }
        }

        for (int i = 1; i < V.getMax_x() - 1; ++i)
        {
            for (int j = 1; j < V.getMax_y() - 1; ++j)
            {
                if (V[i][j].memo)
                    continue;

                if(std::abs(V[i][j].v - V[i][j].prev_v) > c)
                    check = 1;
            }
        }

        if (check == 0)
            break;

        for (int i = 1; i < V.getMax_x() - 1; ++i)
        {
            for (int j = 1; j < V.getMax_y() - 1; ++j)
            {
                V[i][j].prev_v = V[i][j].v;
            }
        }

    }
}

int main()
{
    Matrix<init> mat(12,12);

    for (int i = 0; i < 12; ++i)
    {
        mat[0][i].v = 1;
        mat[0][i].memo = 1;
    }

    Laplace2D_solve(mat, 1e-8);

    for (int i = 1; i < 11; ++i)
    {
        for (int j = 1; j < 11; ++j)
        {
            printf("%.6f    ", mat[i][j].v);
        }
        std::cout<<std::endl;
    }

    return 0;
}

코드는 아주 쉽다. 3차원으로도 확장해 볼 수 있을 것이다.

이 시뮬레이션의 경계조건은 한 변이 모두 1의 값을 갖는 것을 시뮬레이션 한 것이다.

 

위의 코드에서 Laplace2D_solver부분을 보자

코드 설명을 하자면:

1. 외곽 부분은 사용하지 않는다(초기조건을 넣거나 0으로 초기화한다.)

2. memo: 만약 지정한 값이 초기조건이라면, 값을 바꾸지 않는다.

3. 값을 바꿀 때

이렇게 주변 네 개 점의 평균값을 취한다.

4. 모든 값이 이전의 계산값과 일정 이하의 오차가 나지 않는다면 루프를 탈출하지 않는다.

아-주 간단하다.

 

대충 이런 그림을 그릴 수 있다.

그림은 GNUPlot등으로 그릴 수 있으나 생략하겠다.

 

라플라스방정식은 전자기학, 중력 등에 쓰인다.

나중에 자세히 설명하겠다.