НАУЧНАЯ БИБЛИОТЕКА - РЕФЕРАТЫ - Решение задачи Дирихле для уравнения Лапласа методом сеток
Решение задачи Дирихле для уравнения Лапласа методом сеток
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 . Для данного решения выполняется принцип
максимума.
|