VPF::[Численные методы] Конечные разности

Kostik88

Дата 4.12.2008, 20:31 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 76
Регистрация: 9.5.2008

Репутация: 2
Всего: 2

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

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

Есть конечно разностное уравнение вида  (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 (ссылка)

(нет голосов)

Загрузка ...

Опытный

Профиль
Группа: Участник
Сообщений: 375
Регистрация: 14.8.2007
Где: Москва

Репутация: 1
Всего: 13

Цитата

задаю количество узлов по 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

--------------------

Я Критик Джим, раньше писал моды для IPB 1.3.x-2.3.x, а теперь продвигаю Node.js в массы.

Kostik88

Дата 5.12.2008, 15:24 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 76
Регистрация: 9.5.2008

Репутация: 2
Всего: 2

Цитата

Видимо всё-таки 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:04

Kostik88

Дата 5.12.2008, 22:03 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 76
Регистрация: 9.5.2008

Репутация: 2
Всего: 2

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

Kostik88

Дата 5.12.2008, 23:56 (ссылка)

(нет голосов)

Загрузка ...

Шустрый

Профиль
Группа: Участник
Сообщений: 76
Регистрация: 9.5.2008

Репутация: 2
Всего: 2

1 Пользователей читают эту тему (1 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Центр помощи \| Следующая тема »