Управляющие структуры в FreeFem++

Условные операторы if и if/else. Эта структура работает аналогично C++. В первом варианте if, если условие истинно, во выполняется оператор. Во втором варианте if/else, если условие истинно, во выполняется первый оператор, а если ложно, то второй.

real a = 1.5, b = 2.;
if ( a < b )
  cout << a << " < " << b << endl;

if ( a < b )
  cout << a << " < " << b << endl;
else
  cout << b << " < " << a << endl;

Если операторов несколько, то они помещаются в фигурные скобки { ... }. В следующем примере, проверяется целочисленное и вещественное деление (сравните с m = 3).

int m = 4;
real d1, d2;
if ( m / 2 == m / 2. ) {
  d1 = m / 2;
  d2 = m / 2.;
  cout << " m / 2 = m / 2., " << d1 << " = " << d2 << endl;
} else {
  d1 = m / 2;
  d2 = m / 2.;
  cout << " m / 2 != m / 2., " << d1 << " != " << d2 << endl;
}

Иногда бывают случаи когда во время отладки нужно отключать или включать выполнение части кода. Можно поместить этот код в конструкцию if( 0 ) { ... }. Для его включения достаточно поменять ее на if( 1 ) { ... }. Более универсально ввести логическую переменную, например, bool doSamething и присваивая ей соответствующее значение включать или отключать выполнение части кода if( doSamething ) { ... }.

Еще одним типом условия является арифметическое условие (тернарный оператор) - условный оператор используемый в арифметических выражениях. Его формат условие?выражение 1:выражение 2. Если условие истинно, то подставляется выражение 1, иначе - выражение 2. Пример определения максимального из двух чисел:

real a1 = 3.5, a2 = 2.5, amax;
amax = a1 > a2 ? a1 : a2;
cout << "max(" << a1 << ", " << a2 << ") = " << amax << endl;

Цикл, управляемый счетчиком (for). Эта структура работает аналогично C++. Структура for имеет три аргумента: инициализация счетчика, условие повторения и правило обновления счетчика. Переменную счетчика можно описать при инициализации (int i = 0), тогда область существования этой переменной находится внутри структуры цикла.

for (int i1 = 0; i1 < 3; i++)
  cout << " i1 = " << i1 << endl;

Если переменная не была описана до начала цикла, то по его окончанию ее не больше не существует, например, в следующем примере компилятор выдаст ошибку на последней строке

for (int i1 = 0; i1 < 3; i1++)
  cout << " i1 = " << i1 << endl;
cout << i1 << endl; // Ошибка. i1 не существует

Если переменная счетчика описана заранее и в области инициализации она не описывается повторно, то по окончанию цикла переменная имеет последнее значение счетчика не удовлетворяющее условию повтора (например, если правило обновления счетчика - увеличение на единицу i++, то по завершению цикла переменная будет иметь значение на единицу больше, чем последнее значение счетчика использованное внутри структуры цикла).

int i2;
for (i2 = 0; i2 < 3; i2++)
  cout << " i2 = " << i2 << endl;          // i2 = 0, 1, 2
cout << " End value i2 = "  << i2 << endl; // i2 = 3

Если переменная счетчика описана заранее и в области инициализации она описывается повторно, то по окончанию цикла переменная имеет то значение, что она имела до начала цикла.

int i3 = 10;
for (int i3 = 0; i3 < 3; i3++)
  cout << " i3 = " << i3 << endl;          // i3 = 0, 1, 2
cout << " End value i3 = "  << i3 << endl; // i3 = 10

Если тело цикла содержит несколько операторов, то они заключаются внутри фигурных скобок.

int s = 0;
for (int i4 = 0; i4 < 5; i4++) {
  s += i4;
  cout << " i4 = " << i4 << ", s = " << s << endl;
}

Цикл с предусловием (while). Эта структура работает аналогично C++. Структура while имеет один аргумент логического типа, цикл повторяется до тех пор, пока условие остается верным.

int j = -2;
while ( j < 2 )
  cout << " j = " << ++j << endl; // j = -1, 0, 1, 2

Структура цикла while часто используется в итерационных методах - поиск решения с заданной точностью. Рассмотрим пример численного решения нелинейного уравнения x²+x-2=0 методом Ньютона с начальным приближением x=0.5 (точное решение x=1). При инициализации переменной точности eps присваивается заведомо большое значение для старта цикла, внутри цикла точность пересчитывается для нового значения x пока не будет найден корень уравнения с заданной в условии цикла точностью.

real eps = 1.e+10;
real xx = 0.5, dx, f, jac;
while ( eps > 1.e-10 ) {
  f = xx * xx + xx - 2.;
  jac = 2.*xx + 1.;
  dx = - f / jac;
  xx += dx;
  eps = abs(xx * xx + xx - 2.);
  cout << " x = " << xx << " dx = " << dx << ", eps = " << eps << endl;
}

Результат расчета

 x = 1.125, dx = 0.625, eps = 0.390625
 x = 1.00481, dx = -0.120192, eps = 0.0144462
 x = 1.00001, dx = -0.00480001, eps = 2.30401e-05
 x = 1, dx = -7.68e-06, eps = 5.89822e-11

Для тренировки можно запустить тот же пример с начальным приближением x=-1.5, для того, чтобы найти второй корень уравнения (точное значение x=-2), а так же запустить для различных значений точности eps.

Еще один пример использования цикла while - это создание бесконечного цикла while (1). Зачем это может понадобится? Например, цикл должен стартовать даже при несоблюдении условия выполнения. Или же необходимо реализовать структуру цикла с постусловием (в FreeFem++ отсутствует структура do/while). Для выхода из цикла используется команда break.

int k = 3;
while ( 1 ) {
  cout << " k = " << --k << endl; // k = 2, 1, 0
  if ( !k ) break;
}
cout << "End with k = " << k << endl;

Прерывание цикла (break). Прерывает текущий цикл (for или while), выполнение программы продолжается со следующего за циклом оператором.

for (int i = 0; i < 5; i++)
  if (i > 3)
    break;
  else
    cout << " i = " << i << endl;  // i = 0, 1, 2, 3

Продолжение цикла (continue). Игнорирует операторы после команды continue и переходит на начало цикла.

for (int i = 0; i < 5; i++)
  if (i < 3)
    continue;
  else
    cout << " i = " << i << endl;  // i = 3, 4

Неявный цикл (for). Данная структура появилась в FreeFem++ начиная с версии 3.43 и не имеет аналоги в C++. Цикл имеет три аргумента - переменная цикла, ссылка на элемент вектора или матрицы и сам вектор или матрица:

for [i, vi : v] { ... }

Для одномерного вектора v размерности n переменная цикла i пробегает значения от 0 до n-1, а ссылка указывает на элемент вектора v(i). В следующем примере массиву присваивается значение половины квадрата его индекса:

real [int] v(5);
for [i, vi : v] vi = 0.5 * i^2;
cout << " v = " << v << endl;

Для двумерного массива u размерности n×m переменные цикла пробегают значения i=0..n-1 и j=0..m-1, а ссылка указывает на элемент массива u(i,j).

real [int, int] u(3, 5);
for [i, j, uij : u] {
  uij = 1.5 + i - j;
  if( uij < 0. ) uij = 0.;
}
cout << " u = " << u << endl;

Обработка исключительных ситуаций (try/catch). В случае фатальной ошибки выполнение программы останавливается. Если в коде есть опасные участки, где может произойти ошибка, то они могут быть помещены в блок try. При возникновении ошибки выполнение переключается на блок catch, где можно ее исправить без аварийного прерывания программы.

real c;
try{
    c = 1./0.;
}
catch (...) // обработка всех исключительных ситуаций
{
    cout << "Catch an ExecError" << endl;
    c = 0.;
}
cout << " c = " << c << endl;

Приведем еще один пример обработки исключительной ситуации при появлении ошибки в численной схеме. Метод Ньютона для решения нелинейных уравнений обеспечивает хорошую (квадратичную сходимость), но если выбрано начальное приближение в точке экстремума, то метод не сходится, т.к. якобиан равен нулю. Эту ситуацию можно обработать внутри команды исключив ее аварийное завершение.

real eps = 1.e+10;
real xx = -0.5, dx, f, jac;
while ( eps > 1.e-10 ) {
  f = xx * xx + xx - 2.;
  jac = 2.*xx + 1.;
  try {
    dx = - f / jac;
    xx += dx;
    eps = abs(xx * xx + xx - 2.);
  }
  catch (...) {
    cout << " Error: extremum" << endl;
    dx = 0.;
    eps = 1.e10;
    xx += 0.1;
  }
  cout << " x = " << xx << ", dx = " << dx << ", eps = " << eps << endl;
}

Результат выполнения программы:

2.25/0 : d d d
  current line = 107
Exec error :  Div by 0
   -- number :1
 Error: extremum
 x = -0.4, dx = 0, eps = 1e+10
 x = 10.8, dx = 11.2, eps = 125.44
 x = 5.24956, dx = -5.55044, eps = 30.8074
 x = 2.57045, dx = -2.67911, eps = 7.17764
 x = 1.40162, dx = -1.16883, eps = 1.36616
 x = 1.04241, dx = -0.359209, eps = 0.129031
 x = 1.00058, dx = -0.0418276, eps = 0.00174955
 x = 1, dx = -0.000582957, eps = 3.39839e-07
 x = 1, dx = -1.1328e-07, eps = 1.33227e-14

real eps = 1.e+10;
real xx = -0.4999999, dx, f, jac;
int ic = 0;
while ( eps > 1.e-10 ) {
  ic++;
  f = xx * xx + xx - 2.;
  jac = 2.*xx + 1.;
  if ( abs(jac) < 0.25 ) {
    cout << " Updating x value from " << xx << " to ";
    xx += 0.1;
    cout << xx << endl;
    continue;
  }
  dx = - f / jac;
  xx += dx;
  eps = abs(xx * xx + xx - 2.);
  cout << " x = " << xx << ", dx = " << dx << ", eps = " << eps << ", jac = " << jac << endl;
}
cout << " Total steps: " << ic << endl;

 Updating x value from -0.5 to -0.4
 Updating x value from -0.4 to -0.3
 x = 5.225, dx = 5.525, eps = 30.5256, jac = 0.4
 x = 2.55901, dx = -2.66599, eps = 7.10751, jac = 11.45
 x = 1.39727, dx = -1.16174, eps = 1.34963, jac = 6.11801
 x = 1.04159, dx = -0.355677, eps = 0.126506, jac = 3.79454
 x = 1.00056, dx = -0.041031, eps = 0.00168355, jac = 3.08318
 x = 1, dx = -0.000560972, eps = 3.1469e-07, jac = 3.00112
 x = 1, dx = -1.04897e-07, eps = 1.06581e-14, jac = 3
 Total steps: 9