Численное решение задачи

Программа, комментарии

Численное решение задачи проводится следующим образом: вводятся параметры задачи (размеры сетки, шаг по времени), заводятся два (в данном случае двумерных) массива для хранения верхнего и нижнего слоев схемы, задается начальное условие; сам расчет производят в цикле, где значение на верхнем слое вычисляется непосредственно через значения на нижнем, т. е. по явной схеме,.

;

Само вычисление производится по формуле:

u_next[i][j] = u[i][j] + ht*(eps* ((u[i — 1][j] - 2 * u[i][j] + u[i + 1][j]) /(hx*hx)+ (u[i][j -1] - 2 * u[i][j] + u[i][j + 1]) / (hy*hy)) — 2 * u[i][j] * (u[i][j] - u[i — 1][j]) / hx — 2 * u[i][j] * (u[i][j]- u[i][j — 1]) / hy);

после чего происходит коррекция области. В нашем случае мы производим сужение области вычислений, записывая нулевые значения в узлы, оказывающиеся вне нее.

int iN = N;

int jM = M;

int i0 = 1;

int j0 = 1;

for (int k = 1; k <= K; k++) //цикл по времени.

{.

//цикл пробегающий двумерный массив.

for (int i = 0; i <= N; i++).

{.

for (int j = 0; j <= M; j++).

{.

if ((i >= i0 && i = j0 && j < jM)).

//если принадлежит области расчета.

{.

/* формула (явная схема) */.

}.

else.

{.

u_next[i][j] = 0;

//Вне текущий области рассчета всегда 0.

}.

}.

}.

temp = u;

u = u_next;

u_next = temp;

/*Область уменьшается с каждым шагом по времени, уменьшаем область*/.

i0++;

j0++;

iN—;

jM—;

}.

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

Таким образом, размер сетки выбирается в зависимости от того насколько долго нужно считать. Вывод производится в файл, по которому в дальнейшем можно построить график решения в данный момент времени.

if (k%kout == 0).

{.

write_file_gnuplot_small (u, file, path_name, N, M, hx, hy, k, X, realX);

}.

где kout частота вывода в файл.

Дополнительно, для усмотрения движения угловой точки на фронте вводим массив для записи положения вершины в каждую итерацию. Сама запись производится непосредственно после вычисления.

vertex[k] = find_vertex (u_next, N, hx, vertex_flag, jv);

Для поиска вершины вводится функция find_vertex (в ней vertex_flag порог, выше которого точка считается достигнутой, соответственно максимальная из этих точек находится на фронте; jv хранит в себе вершину на предыдущей итерации).

float find_vertex (float**u, int N, float hx, float vertex_flag, int* jv).

{.

int i= 0;

for (i = *jv; i vertex_flag; i++);

*jv = i;

return i*hx;

}.

Результаты

Приложим вычисление при углах на промежутке.

График в моменты + график движения вершины.

Приведем также графики движения на одном рисунке.