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

    Одним из наиболее распространенных подходов численного решения дифференциальных уравнений является метод конечных разностей (метод сеток). Следуя этому подходу, область решения D представляется в виде дискретного (как правило, равномерного) набора (сетки) точек (узлов). Так, например, прямоугольная сетка в области D может быть задана в виде:

где величина N задает количество узлов по каждой из координат области D.

    Обозначим оцениваемую при подобном дискретном представлении аппроксимацию функции u(x, y) в точках (x_i, y_j) u_ij. Тогда, используя пятиточечный шаблон для вычисления значений производных, мы можем представить уравнение Пуассона в конечно-разностной форме:

    Данное уравнение может быть разрешено относительно u_ij:

    Разностное уравнение, записанное в подобной форме, позволяет определять значение u_ij по известным значениям функции u(x, y) в соседних узлах используемого шаблона. Данный результат служит основой для построения различных итерационных схем решения задачи Дирихле, в которых в начале вычислений формируется некоторое приближение для значений u_ij, а затем эти значения последовательно уточняются в соответствии с приведенным соотношением. Так, например, метод Гаусса-Зейделя для проведения итераций уточнения использует правило:

по которому очередное k-ое приближение значения u_ij вычисляется по последнему k-ому приближению значений u_i-1,j и u_i,j-1 и предпоследнему (k-1)-ому приближению значений ui+1,j и ui,j+1. Выполнение итераций обычно продолжается до тех пор, пока получаемые в результате итераций изменения значений u_ij не станут меньше некоторой заданной величины (требуемой точности вычислений). Последовательность решений, получаемых методом сеток, равномерно сходится к решению задачи Дирихле, а погрешность решения имеет порядок h².

    Рассмотренный алгоритм (метод Гаусса-Зейделя) на псевдокоде, приближенном к алгоритмическому языку С++, может быть представлен в виде:

do { 
  dmax = 0; // максимальное изменение значений u 
  for ( i=1; i<N+1; i++ ) 
    for ( j=1; j<N+1; j++ ) { 
      temp = u[i][j]; 
      u[i][j] = 0.25*(u[i-1][j]+u[i+1][j]+u[i][j-1]+u[i][j+1]-h*h*f[i][j]); 
      dm = fabs(temp-u[i][j]); 
      if ( dmax < dm ) dmax = dm; 
    } 
} while ( dmax > eps ); 

(значения u_ij при индексах i, j = 0, N+1 являются граничными, задаются при постановке задачи и не изменяются в ходе вычислений).

    Со следующего шага мы начнем рассматривать использование параллельных алгоримов.

Предыдущий шаг Содержание Следующий шаг