Добрый день! Решаю численно задачу (дифференциальное уравнение в частных производных 4 порядка - бигармоническое уравнение).
Для решения использую следующую разностную схему:
Пытаюсь заполнить матрицу системы размером
, где каждая строка -- это соответствующее уравнение, а столбцы -- это коэффициенты перед переменными..
Т.е. каждая строка матрицы -- уравнение с переменными
и т.д.
Возникла проблема именно программно, как в
эту матрицу заполнить циклически? Т.е. у меня каждый элемент матрицы
,
-номер уравнения,
-- номер переменной. С номером уравнения еще ладно, а вот номер переменной
-- как получить этот номер. Сижу уже 2 день, не могу никак сдвинуться. Т.е. вопрос в следующем -- как системе разностных уравнений сопоставить матрицу коэффициентов?
Попробую объяснить простыми словами: вот у меня есть самое верхнее уравнение и для каждых
. Если расписывать например при фиксированных индексах, получится линейное уравнение относительно
. Понятно, что некоторые коэффициенты будут равны нулю. Я пытаюсь грубо говоря что сделать -- каждая строка матрицы
-- это каждое уравнение системы. А каждый столбец -- это коэффициент при соответствующей переменной. Так вот я не понимаю, как найти связь между позицией переменной
и номером столбца, где она должна стоять. Я например, пробую на самом простом случае найти эту закономерность, т.е. если взять
Ну вот например, при
самое верхнее уравнение будет только одно. Переменная
-- первая по счету в уравнении, потом
-- вторая и т.д.
Внизу прилагаю код того, как хотя бы я пытался заполнить матрицу коэффициентов системы для первого уравнения. Я просто не могу выявить закономерность между переменной
и номером столбца, в которой она стоит.
//Начальные данные для задачи
TIP T, h_x, h_z, M, N, a, b;
cout << "Введите a = ";
cin >> a;
cout << "Введите b = ";
cin >> b;
cout << "Введите шаг по x-координате h_x =";
cin >> h_x;
cout << "Введите шаг по x-координате h_z =";
cin >> h_z;
int m = a / h_x ; // количество точек сетки по пространству
int n = b / h_z ; // количество точек сетки по времени
cout << "m = " << n << endl;
cout << "n = " << m << endl;
Matrix A((n+1)*(m+1), vector <TIP>((n + 1) * (m + 1))); // матрица системы
for (int i = 0; i <= m; i++)
{
for (int j = 0; j <= n; j++)
{
A[i][j] = 0.0;
}
}
Menu(); //вызов меню пользователя для выбора решаемой задачи
int selector;
cin >> selector;
vec X(n+1); // сетка по X
vec Z(m+1); // сетка по Z
vec B((m+1) * (n+1));
for (int i = 0; i <= n * m; i++)
{
B[i] = 0;
}
for (int i = 0; i <= m; i++)
{
X[i] = i * h_x;
}
for (int j = 0; j <= n; j++)
{
Z[j] = j * h_z;
}
cout << "Сетка по X:" << endl;
out_of_right_vector(X);
cout << endl;
cout << "Сетка по Z:" << endl;
out_of_right_vector(Z);
cout << endl;
cout << "Формирование матрицы системы - 1 этап" << endl;
for (int k = 0; k < (m - 3) * (n - 3) ; k++)
{
for (int i = 2; i <= (m - 2); i++)
{
for (int j = 2; j <= (n - 2); j++)
{
A[k][j + i ] = 20 / pow(h_x, 4);
A[k][j + i + 1 ] = -8 / pow(h_x, 4);
A[k][j + 1 + i ] = -8 / pow(h_x, 4);
A[k][j - 1 + i] = -8 / pow(h_x, 4);
A[k][i + 1 + j + 1 ] = 2 / pow(h_x, 4);
A[k][j + 1 + i -1 ] = 2 / pow(h_x, 4);
A[k][j - 1 + i - 1 ] = 2 / pow(h_x, 4);
A[k][j - 1 + i + 1 ] = 2 / pow(h_x, 4);
A[k][j + i + 2 ] = 1 / pow(h_x, 4);
A[k][j + 2 + i ] = 1 / pow(h_x, 4);
A[k][j + i - 2 ] = 1 / pow(h_x, 4);
A[k][j - 2 + i ] = 1 / pow(h_x, 4);
}
}
}