학교수업, CS/수치해석
편미분 방정식(PDE, 라플라스, 포아송)
빨대도둑
2022. 11. 30. 16:03
편미분 방정식은 상미분 방정식에서처럼 차수나 선형성에 대한 분류도 똑같이 적용될 수 있다. 하지만 ODE와는 달리 PDE는 다음 세 가지로 크게 분류함에 따라 해법도 달라지고 있다. 타원형, 쌍곡선형, 포물선형으로 분류한다.
PDE에서 차분 방정식의 표현
어떤 종류의 PDE든지 미분향을 차분 방정식의 표현으로 고쳐서 생각할 수 있다. 주어진 영역을 사각형의 격자 모양으로 나누어서 차분 방정식들을 도출한 후, 연립 방정식의 해법으로 그 해를 구한다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
void LAP_matrix(int out_mat[][21],int n, int m, int method) {
int nl1, np1, msize, msln, i;
nl1 = n - 1;
np1 = n + 1;
msize = n * m;
msln = msize - n;
/*1. five point methd*/
if (method == 5) {
for (i = 0; i < msize; i++)
out_mat[i][i] = -4;
for (i = 1; i < msize; i++) {
out_mat[i - 1][i] = 1;
out_mat[i][i-1] = 1;
}
for (i = nl1; i < msln; i += n) {
out_mat[i][i + 1] = 0;
out_mat[i + 1][i] = 0;
}
for (i = np1 - 1; i < msize; i++) {
out_mat[i][i - n] = 1;
out_mat[i - n][i] = 1;
}
}
/*2. nine point method*/
else {
for (i = 0; i < msize; i++)
out_mat[i][i] = -20;
for (i = 1; i < msize; i++) {
out_mat[i - 1][i] = 4;
out_mat[i][i - 1] = 4;
}
for (nl1; i < msln; i += n) {
out_mat[i][i + 1] = 0;
out_mat[i + 1][i] = 0;
}
for (i = np1 - 1; i < msize; i++) {
out_mat[i][i-n] = 4;
out_mat[i - n][i] = 4;
out_mat[i][i - n + 1] = 1;
out_mat[i - n + 1][i] = 1;
out_mat[i+1][i - n ] = 1;
out_mat[i-n][i + 1] = 1;
}
for (i = np1 - 1; i < msize; i += n) {
out_mat[i - 1][i + 1] = 0;
out_mat[i][i + n - 1] = 0;
out_mat[i +n -1][i] = 0;
out_mat[i +n][i - 1] = 0;
}
}
}
int main(void) {
int out_mat[21][21], n, m, method, i, j;
n = 7, m = 3;
for (i = 0; i < 21; i++)
for (j = 0; j < 21; j++)
out_mat[i][j] = 0;
while (1) {
printf("determine whether you select 5-point or 9-point \n");
printf("for 5-point method type 5, for 9-point method type 9: ");
scanf_s("%d", &method);
if ((method == 5) || (method == 9))break;
else printf("retype \n");
}
LAP_matrix(out_mat, n, m, method);
printf("\t\t the output for %d-point \n", method);
for (i = 0; i < 21; i++) {
for (j = 0; j < 21; j++)
printf("%3d", out_mat[i][j]);
printf("\n");
}
}
라플라스 방정식
좀더 정확한 해를 구하기 위해서 h를 작게 하면 좋다. 그런데 h가 작아짐에 따라 증가하는 방정식의 개수와 이로 인해 이식들을 푸는데 필요한 엄청난 양의 메모리와 연산량가 소모된다는 문제점이 생긴다.
이러한 특성을 갖는 행렬식을 풀기 위해서 보통 가우스-자이달 반복법을 사용한다. 그리고 이러한 방법이 라플라스 방정식에 응용될 때 이를 이를 리브만의 방법이라고 부른다.
수렴의 속도는 사실상 초기값의 추정을 어떻게 하느냐에 따라 달려있으므로 이에 대한 많은 연구가 있었다. 내부점들의 값을 선형보간을 이용해 초기값을 정하는 것도 좋은 방법이다. 그러나 이 리브만의 방법이 소거법에서의 메모리 문제를 해결했다면, 그 부산물로서 특히 많은 점의 해를 구할 때 수렴 속도에 대한 문제가 또 생기가 된다. 이에 대한 대책으로 SOR(Successive OverRelaxation)이라는 방법을 사용할 수 있다. 이것은 Southwell의 안정화 방식인데 보다 빠르게 수렴할 수 있게 한다.
#include <stdio.h>
#include <math.h>
#include <conio.h>
#define SQ(X) ((X)*(X))
#define TOL 0.001
float funct(float x, float y) {
float out;
out = 0.0 * x + 0.0 * y;
return out;
}
int main(void) {
int row, column, i, j, n;
float h, u[100][100], temp, w, sum, average, x, y, residue, chgmax;
h = 1.25, row = 16, column = 8;
for (i = 0; i <= column; i++) {
temp = 0.0;
u[i][0] = (float)temp;
}
for (i = 0; i <= column; i++) {
temp = 100.0;
u[i][row] = (float)temp;
}
for (j = 1; j < row; j++) {
temp = 0.0;
u[0][j] = (float)temp;
}
for (j = 1; j < row; j++) {
temp = 0.0;
u[column][j] = (float)temp;
}
/*find initial matrix*/
sum = 0.0;
for (i = 0; i <= column; i++)
sum += u[i][0] + u[i][row];
for (j = 1; j < row; j++)
sum += u[0][j] + u[column][j];
if (sum != 0.0)
average = sum / (2.0 * (float)row + 2.0 * (float)column);
x = y = 0.0;
for (i = 1; i < column; i++)
for (j = 1; j < row; j++)
u[i][j] = average + h * h * funct(x, y);
/*S.O.R method for the various w, limit : 100 itration, tolerance : 0.001*/
printf("this program is finding the relation between w and number of iterations");
for (w = 1.0; w < 1.8; w += 0.1) {
for (n = 0; n < 100; n++) {
chgmax = 0.0;
for (i = 1; i < column; i++) {
y = (float)i * h;
for (j = 1; j < row; j++) {
x = (float)j * h;
residue = w * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1]
- 4.0 * u[i][j] + h * h * funct(x, y)) / 4.0;
if (chgmax < fabs(residue))
chgmax = fabs(residue);
u[i][j] += residue;
}
}
if (chgmax < TOL)break;
}
printf("\n w=%f, count = %d", w, n + 1);
printf("\t press the key");
for (i = 1; i < column; i++)
for (j = 1; j < row; j++)
u[i][j] = average + h * h * funct(x, y);
}
}
포아송 방정식
라플라스의 방정식 풀이법과 거의 같다.
<수식>
#include <stdio.h>
#include <math.h>
#include <conio.h>
#define M_PI 3.14159265358979323846
#define SQ(X) ((X)*(X))
#define TOL 0.001
float funct(float x, float y) {
float out;
out = 0.0 * x + 0.0 * y + 2.;
return out;
}
float w_opt(int row, int column) {
float out, temp;
temp = cos(M_PI / (float)row) + cos(M_PI / (float)column);
out = 4.0 / (2.0 + sqrt(4.0 - SQ(temp)));
return out;
}
int main(void) {
int row, column, i, j, n;
float h, u[100][100], temp, w, sum, average, x, y, residue, chgmax;
printf("\n\n enter the step size h=");
scanf_s("%f", &h);
printf("\n enter the number of meshes in the x-direction : ");
scanf_s("%d", &row);
printf("\n enter the number of neshes in the y-direction : ");
scanf_s("%d", &column);
/*for boundary condition*/
printf("\n");
printf("suppose the matrix notation\n");
for (i = 0; i <= column; i++) {
printf("enter the boundary condition in left side, u[%d][0]=", i);
scanf_s("%f", &temp);
u[i][0] = (float)temp;
}
printf("\n");
printf("suppose the matrix notation\n");
for (i = 0; i <= column; i++) {
printf("enter the boundary condition in right side, u[%d][%d]=", i, row);
scanf_s("%f", &temp);
u[i][row] = (float)temp;
}
printf("\n");
printf("suppose the matrix notation\n");
for (j = 1; j < row; j++) {
printf("enter the boundary condition in upper side, u[0][%d]=", j);
scanf_s("%f", &temp);
u[0][j] = (float)temp;
}
printf("\n");
printf("suppose the matrix notation\n");
for (j = 1; j < row; j++) {
printf("enter the boundary condition in lower side, u[%d][%d]=", column, j);
scanf_s("%f", &temp);
u[column][j] = (float)temp;
}
/*find w_opt*/
w = w_opt(row, column);
/*finding initial matrix*/
sum = 0.0;
for (i = 0; i <= column; i++)
sum += u[i][0] + u[i][row];
for (j = 1; j < row; j++)
sum += u[0][j] + u[column][j];
if (sum != 0.0)
average = sum / (2.0 * (float)row + 2.0 * (float)column);
x = y = 0.0;
for (i = 1; i < column; i++)
for (j = 1; j < row; j++)
u[i][j] = average + h * h * funct(x, y);
/*S.O.R method for the various w, limit : 100 itration, tolerance : 0.001*/
for (n = 0; n < 100; n++) {
chgmax = 0.0;
for (i = 1; i < column; i++) {
y = (float)i * h;
for (j = 1; j < row; j++) {
x = (float)j * h;
residue = w * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1]
- 4.0 * u[i][j] + h * h * funct(x, y)) / 4.0;
if (chgmax < fabs(residue))
chgmax = fabs(residue);
u[i][j] += residue;
}
}
if (chgmax < TOL)break;
}
/*output the solution matrix*/
printf("\n--- the solution matrix with overrelaxation facotr w=%f---", w);
for (i = 0; i < column; i++) {
for (j = 0; j <= row; j++)
printf("%f\t", u[i][j]);
printf("\n");
}
}