Решение задачи Дирихле для уравнения Лапласа методом сеток
1.
ПОСТАНОВКА ЗАДАЧИ
Решить численно задачу Дирихле для уравнения Лапласа :
(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 . Точность решения СЛАУ e=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]|< e ,взятый по всем точкам области Dh* ,то задача решена;
если нет , то
выполнять шаг 5 ( пересчитывать функцию F(N)*[i,j] через значения F(N-1)*[i,j])
до тех пор, пока не выполнится указанное условие.
3.ТЕКСТ ПРОГРАММЫ
#include
<stdio.h>
#include <fstream.h>
#include <conio.h>
#include <iostream.h>
#include <math.h>
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;i<13;i++)
FillF_();
CountDif();
}
cout<<(0-arrayP[7][17].F_);
MakeFile();
getchar();
} //MAIN END
int IsLineFit(float par,char Axis) // does the line
belong to the defined area
{
switch(Axis) {
case 'y': if ((par>8.0) || (par<2.0))
return 1;
else return 0;
case 'x': if (par<1.9) return 1;
else if (par>4.0) return 1;
else return 0;
}
}
void CreateNet(void) // Creation of Net (area
D )
{
x = p0.xx;
i=0;
arrayX[i]=x;
while (!IsLineFit(x,'x'))
{
x += h;
i++;
arrayX[i] = x;
}
x = p0.xx-h;
i++;
arrayX[i]=x;
while (!IsLineFit(x,'x'))
{
x -= h;
i++;
arrayX[i] = x;
}
for (i=0;i<13;i++) { printf("%g
",arrayX[i]); }
printf("\n");
y = p0.yy;
i = 0;
arrayY[i]=y;
while (!IsLineFit(y,'y'))
{
y += h;
i++;
arrayY[i] = y;
}
y = p0.yy - h;
i++;
arrayY[i]=y;
while (!IsLineFit(y,'y'))
{
y -= h;
i++;
arrayY[i] = y;
}
for(i=0;i<33;i++) { printf("%g
",arrayY[i]);}
printf("\n");
} // end CreateNet
void RegArrayX() // Regulation of arrays X & Y
{
int LastUnreg = 13 ;
while (LastUnreg != 0) {
for(i=0;i<LastUnreg-1;i++) {
if (arrayX[i]>arrayX[i+1]) {double
tmp=arrayX[i];
arrayX[i]=arrayX[i+1];
arrayX[i+1]=tmp;}}
LastUnreg=LastUnreg-1; }
for (i=0;i<13;i++) { printf("%g
",arrayX[i]);
}
void RegArrayY()
{
int LastUnreg = 33 ;
while (LastUnreg != 0) {
for(i=0;i<LastUnreg-1;i++) {
if (arrayY[i]>arrayY[i+1]) {
tmp=arrayY[i];
arrayY[i]=arrayY[i+1];
arrayY[i+1]=tmp;}}
LastUnreg=LastUnreg-1; }
for (i=0;i<33;i++) { printf("%g
",arrayY[i]); }
printf("\n");} // End of
Regulation
void CrMtrD(void) //Create general Matrix
{
for(i=0;i<13;i++)
for(j=0;j<33;j++)
{arrayP[i][j].BelongsToDh_=0;
arrayP[i][j].BelongsToDh=0;}
for(i=0;i<13;i++)
for(j=0;j<33;j++) {
arrayP[i][j].xx=arrayX[i];
arrayP[i][j].yy=arrayY[j];
}
// printf("%g
%g",arrayP[12][0].xx,arrayP[12][0].yy);
// printf("\n");
}
int IsFit(point Par) //does point belong to area D?
{
if ((Par.xx<=4) &&
(Par.xx>=1.99) && (Par.yy>=Par.xx)
&& (Par.yy<=Par.xx+4)) return
1;
else return 0;
}
void CreateDh_(void) //Create area Dh_
{
for(i=0;i<13;i++)
for(j=0;j<33;j++)
if (IsFit(arrayP[i][j]))
arrayP[i][j].BelongsToDh_=1;
cout << arrayP[1][1].BelongsToDh_<<
"\n";
cout << arrayP[1][1].xx <<
" " << arrayP[1][1].yy<<"\n";
}
void FillF(void) // calc function F(x,y) at area Dh_
{
for(i=0;i<13;i++)
for(j=0;j<33;j++)
if (arrayP[i][j].BelongsToDh_==1)
arrayP[i][j].F=arrayP[i][j].xx*pow(arrayP[i][j].yy,2);
else arrayP[i][j].F=0;
}
int IsInner(int i,int j) //Is point inner?
{
if ((arrayP[i-1][j].BelongsToDh_==1)
&&
(arrayP[i+1][j].BelongsToDh_==1)
&&
(arrayP[i][j+1].BelongsToDh_==1)
&&
(arrayP[i][j-1].BelongsToDh_==1)) return
1;
else return 0;
}
void CreateDh(void) //Create area Dh
{
for(i=0;i<13;i++)
for(j=0;j<33;j++)
if ((arrayP[i][j].BelongsToDh_==1)
&&
IsInner(i,j))
arrayP[i][j].BelongsToDh=1;
}
void FillF_() //calc new appr. values of F
{
for(i=0;i<13;i++)
for(j=0;j<33;j++) {
if (arrayP[i][j].BelongsToDh==1)
arrayP[i][j].F_=(arrayP[i-1][j].F+arrayP[i+1][j].F+
arrayP[i][j-1].F+arrayP[i][j+1].F)/4;
}
void CountDif() // find maximal difference abs(F-F_)
{
k=0;
for(i=0;i<13;i++)
for(j=0;j<33;j++)
{ if (arrayP[i][j].BelongsToDh==1) {
diff[k]=fabs(arrayP[i][j].F_-arrayP[i][j].F);
k++;}}
E1=diff[0];
for (k=1;k<500;k++) {
if (diff[k]>E1) 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<<arrayP[i][j].xx<<"
"<<arrayP[i][j].yy<<
"
"<<arrayP[i][j].F_<<"\n";}}
f.close() ;
}
3.
ГРАФИКИ РЕШЕНИЙ
РИС.1 шаг h=0.2
РИС.2 шаг h=0.1
5.ВЫВОД
Функция f(x,y) является
неотрицательной в области D. Полученное решение лежит целиком над плоскостью XOY .
Для данного решения выполняется принцип максимума.