Математика

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


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;i<13;i++)
                 for(j=0;j<33;j++) arrayP[i][j].F=arrayP[i][j].F_;
                 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;iarrayX[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]);
            } printf("\n");
            }
   void RegArrayY()
            {
             int LastUnreg = 33 ;
             while (LastUnreg != 0)       {
                 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<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;
                 else arrayP[i][j].F_=0; }
            }
   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<

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