Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Центр помощи > [Численные методы] Конечные разности

Автор: Kostik88 4.12.2008, 20:31

Здравствуйте!

сталкнулся с такой проблемой:

Есть конечно разностное уравнение вида  (h[i+1][j]^3)*P[i+1][j]  -  (  h[i+1][j] ^3  +  h[i-1][j] ^3  +  h[i][j+1] ^3 +  h[i][j-1] ^3)*P[i][j]  +  (h[i-1][j] ^3)*P[i-1][j]  +  (h[i][j+1] ^3)*P[i][j+1] + (h[i][j-1]^3)*P[i][j-1] = A*b*(h[i+1][j]  + h[i][j] );

где h - известные величины, A  и  b  - постоянные величины. Требуется найти P в каждом узле расчётной области.

К примеру, есть прямоугольная область, задаю количество узлов по i скажем 5 (т.е. 0<=i<5) и по j столько же. Получается матрица области - 5 на 5.
Плюс ко всему даны граничные условия:
                                                                       P[1][j] - P[1][j-1] = 0

                                                                       P[i][1] + P[i+1][1] = 2F

                                                                       P[i+1][4] + P[i][4] = 2F

                                                                       P[4][j] + P[4][j+1] = 2F

, где F - постоянная (любое значение, к приперу 771)

Так вот вопрос в том как найти все значения величины P в каждом узле. Я знаю что должна получится матрица системы уравнений размером 25х25, но как её получить + как учитывать такие граничные условия не могу сообразить?

Автор: Sannis 5.12.2008, 02:35

Цитата

задаю количество узлов по i скажем 5 (т.е. 0<=i<5)

Видимо всё-таки 1<=i<=5 smile

См. граничные условия.

Теоретическая часть: http://csa.ru/~stan/multigrid/

Практически нужно из уравнения выразить P[i][j] через остальные переменные. Далее итеративно проходить по всем внутренним точкам области и рассчитывать P[i][j] на следующем шаге.

Единственное, чего нам не хватает - удобных граничных условий:

P[1][j] - P[1][j-1] = 0 => P[1][j] = P[1][j-1] = C и не зависит от j => P[1][1] = C

P[i][1] + P[i+1][1] = 2F => P[i+1][1] = 2F - P[i][1] => P[2][1] = 2F-C, P[3][1] = C, P[4][1] = 2F-C, P[5][1] = C.

P[i+1][4] + P[i][4] = 2F - аналогично предыдущему.

P[4][j] + P[4][j+1] = 2F - а вот отсюда получаем, что если точек нечётное число, то мы не сможем определить отсюда C. а если четное, то C = F smile

Неоднозначность нехорошая, я бы переспросил условие smile

Автор: Kostik88 5.12.2008, 15:24

Цитата

Видимо всё-таки 1<=i<=5 См. граничные условия.

Нет так не пойдёт, ведь все граничные условия записаны в неявном виде.

Можно ГУ записать в таком виде:
                                                                      P[0][j] - P[0][j-1] = 0  (условие симметрии)

                                                                       P[i][0]  = F

                                                                        P[i][4] = F

                                                                       P[4][j] = F

Но вот с первым граничным условием не знаю что делать.   smile

записать его чтоли так P[0][j] = P[0][j-1]?

И ещё вопрос:
Sannis, я посмотрел сайт который вы мне кинули... там есть пример алгоритма для уравнения лапласа.... так вот я попробовал закодить его на С++, но решения не получается Оо

Может где то я напутал?

Код


#include<stdio.h>
#include<iostream>
#include<stdio.h>
#include<stdlib.h>
#include<cstring>
#include<math.h>
using namespace std;
void main(  )
{
     double a[10][10];
    

    
    
    
 double D;
double p;



for(int i=0;i<10;i++)
{for(int j =0;j<10;j++)
  {
   
      a[i][j] = 0;

}
}


for(int i=0;i<10;i++)
{
 for(int j=0;j<10;j++)
 {
  a[i][0] = 1;
  a[9][j]=1;
  a[i][9]=2;
  a[0][j]=2;
 
 }
}

for(int i=0;i<10;i++)
{ cout<<"\n";
 for(int j=0;j<10;j++)
 {
  cout<<a[i][j]<<" ";
 
 }
}




double eps = 0.01;

int num_iter = 0;
int max_iter = 1000;

D = 0;

while(D > eps)
{
    if(num_iter > max_iter) { break;}

    num_iter+=1;
 

 for(int i =1;i<=8;i++)
 { for(int j=1;j<=8;j++)
     { p = (a[i-1][j] + a[i+1][j] + a[i][j-1] + a[i][j+1])/4;  
    
 if(fabs(p - a[i][j]) > D)
 {
  D = fabs(p - a[i][j]);
 }
 else
  {
   a[i][j] = p; 
  }
  

     }
 
 
 }


}



cout<<"\n\n";
for(int i=0;i<10;i++)
{ cout<<"\n";
 for(int j=0;j<10;j++)
 {
  cout<<a[i][j]<<" ";
 
 }
}

}

А вот собственно сайт: http://csa.ru/~stan/multigrid/lesson1.html

Автор: Kostik88 5.12.2008, 22:03

Я в коде заменил вместо D=0 поставил D=1.... вроде чтото считает, но невязку fabs(p - a[i][j]) выдаёт какуюто странную.... она больше единицы Оо и условие while(D>eps) не выполняется, т.е. программа пробигает все итерации пока не дойдёт до максимальной и следовательно выходит из цикла... я немогу понять в чём проблема. Помоготи резобраться пожалуйста.

Автор: Kostik88 5.12.2008, 23:56