Решение задачи Дирихле для уравнения Лапласа методом сеток

Решение задачи Дирихле для уравнения Лапласа методом сеток

1. ПОСТАНОВКА ЗАДАЧИ

Решить численно задачу Дирихле для уравнения Лапласа :

[pic]

(x,y) ?D ; u|Г=xy2=f(x,y) ;

область D ограничена линиями: x=2 , x=4 , y=x , y=x+4 ;

(x0, y0 ) = (3, 5) .

Следует решить задачу сначала с шагом по x и по y : h=0.2, потом с

шагом h=0.1 . Точность решения СЛАУ ?=0.01 .

2. ОПИСАНИЕ МЕТОДА РЕШЕНИЯ ПОСТАВЛЕННОЙ ЗАДАЧИ

Поставленная задача решается численно с помощью программы,

реализующей метод сеток , разработанный для численного решения задачи

Дирихле для уравнений эллептического типа.

Программа написана на языке C++ , в среде Borland C++ версии 3.1. Ниже

описан алгоритм работы этой программы.

1. На первом шаге область D дискретизируется. Она заменяется на

область Dh путем разбиения области D параллельными прямыми по следующему

правилу: yi=y0 ± ih, xj=x0 ± ih , i,j=0,1,2…. РР

Разбиение производится до тех пор, пока текущая прямая не будет лежать

целиком вне области D. Получается множество точек (xi,yj).

2. За область Dh принимают те точки множества (xi,yj) , которые

попали внутрь области D, а также дополняют это множество граничными

точками.

3.Во всех точках области Dh вычисляются значения функции f(xi,yj)

.

4. За область Dh* принимаются все внутренние точки области Dh, т.е.

удовлетворяющие требованию:

(xi,yj) ? Dh* , если (xi+1,yj) ? Dh , (xi-1,yj) ? Dh ,

(xi,yj+1) ? Dh , (xi,yj-1) ? Dh .

5. Во всех точках области Dh* вычисляется функция F(N)*[i,j] (

индекс N обозначает номер итерации, на которой производится вычисление):

F(N)*[i,j]=(f(xi+1,yj) + f(xi-1,yj) + f(xi,yj+1)

+ f(xi,yj-1))/4

6. Теперь если max | F(N+1)*[i,j] - F(N)*[i,j]|< ? ,взятый по всем

точкам области Dh* ,то задача решена;

если нет , то выполнять шаг 5 ( пересчитывать функцию

F(N)*[i,j] через значения F(N-1)*[i,j]) до тех пор, пока не выполнится

указанное условие.

3.ТЕКСТ ПРОГРАММЫ

#include

#include

#include

#include

#include

int i,j,k; // Variables

float h,x,y,tmp,E1;

struct point {

float xx;

float yy;

int BelongsToDh_;

int BelongsToDh;

float F;

float F_;

} p0,arrayP[13][33];

float arrayX[13];

float arrayY[33];

float diff[500];

void CreateNet(void); // Procedure Prototypes

int IsLineFit(float Param);

void CrMtrD(void);

void RegArrayX();

void RegArrayY();

void CreateDh_();

int IsFit(point Par);

void FillF();

void CreateDh();

int IsInner(int i,int j);

void FillF_();

void CountDif();

void MakeFile();

void main(void) //MAIN

{

clrscr();

p0.xx = 3;

p0.yy = 5;

h = 0.2;

p0.BelongsToDh_=1;

p0.BelongsToDh=1;

CreateNet();

RegArrayX();

RegArrayY();

CrMtrD();

CreateDh_();

FillF();

CreateDh();

FillF_();

CountDif();

while (E1>=0.005) {

for(i=0;iarrayX[i+1]) {double tmp=arrayX[i];

arrayX[i]=arrayX[i+1];

arrayX[i+1]=tmp;}}

LastUnreg=LastUnreg-1; }

for (i=0;iarrayY[i+1]) { tmp=arrayY[i];

arrayY[i]=arrayY[i+1];

arrayY[i+1]=tmp;}}

LastUnreg=LastUnreg-1; }

for (i=0;i=1.99) && (Par.yy>=Par.xx)

&& (Par.yyE1) E1=diff[k];}

}

void MakeFile()

{

ofstream f;

FILE *f1=fopen("surf.dat","w1");

fclose(f1);

f.open("surf.dat",ios::out,0);

for(i=0;i<13;i++)

for(j=0;j<33;j++) { if (arrayP[i][j].BelongsToDh==1) {

f<

" "<

f.close() ;

}

3. ГРАФИКИ РЕШЕНИЙ

РИС.1 шаг h=0.2

[pic]

РИС.2 шаг h=0.1

5.ВЫВОД

Функция f(x,y) является неотрицательной в области D. Полученное решение

лежит целиком над плоскостью XOY . Для данного решения выполняется принцип

максимума.