Заполнение матрицы системы

Решаю численно задачу (дифференциальное уравнение в частных производных 4 порядка - бигармоническое уравнение). Такое уравнение возникло в курсовой работе.
Для решения использую следующую разностную схему:

Пытаюсь заполнить матрицу системы размером (M+1)^2, где каждая строка – это соответствующее уравнение, а столбцы – это коэффициенты перед переменными…

Т.е. каждая строка матрицы – уравнение с переменными y_{0}^{0}, y_1^0 , y_2^0 и т.д.

Возникла проблема именно программно, как в С++ эту матрицу заполнить циклически? Т.е. у меня каждый элемент матрицы A[i][j], i-номер уравнения, j-- номер переменной. С номером уравнения еще ладно, а вот номер переменной j – как получить этот номер. Сижу уже 2 день, не могу никак сдвинуться. Просто именно это меня и останавливает в дальнейшем продвижении в работе.

Внизу прилагаю код того, как хотя бы я пытался заполнить матрицу коэффициентов системы для первого уравнения. Я просто не могу выявить закономерность между переменной y[i][j] и номером столбца, в которой она стоит.

//Начальные данные для задачи


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);

		}
	}
}

А не программно как? Как эти переменные вычисляются?

Так может её и нет?

Попробую объяснить простыми словами: вот у меня есть самое верхнее уравнение и для каждых i,j=2…M-2! . Если расписывать например при фиксированных индексах, получится линейное уравнение относительноy_0^0, y_1^0… Понятно, что некоторые коэффициенты будут равны нулю. Я пытаюсь грубо говоря что сделать – каждая строка матрицы А-- это каждое уравнение системы. А каждый столбец – это коэффициент при соответствующей переменной. Так вот я не понимаю, как найти связь между позицией переменной y_i^j
и номером столбца, где она должна стоять. Я например, пробую на самом простом случае найти эту закономерность, т.е. если взять M=2

Ну вот например, при i,j = 2 самое верхнее уравнение будет только одно. Переменная y_0^0 – первая по счету в уравнении, потом y_0^1 – вторая и т.д.

Спустя некоторое время потихоньку матрица начала заполняться необходимыми коэффициентами, правда пока что немного не во всех местах совпадает как при написании системы вручную…
В конце взял и выкинул из программы, где-то вышел за границы массива.

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

В самых нижних уравнениях оставил циклы по строкам (просто иду по каждому уравнению), но все таки сомневаюсь в правильности вычислений позиций переменных в матрице.

Во-первых у вас неправильное обозначение столбца и строки. Правильно y[j][i] - это обще признанное обозначение.
Если кратко, то не надо. Есть такая вещь как граничные условия. Примером такого граничного условия может служить начальное распределение y[0][i] и вы из него находите все остальные используя циклы.

К примеру как тут

А зачем вам их вычислять?
Есть два основных подхода это МКР и МКЭ. Так вот вы сейчас их пытаетесь скрестить разумеется у вас получается ерунда.

Просто в этом то и проблема , что у меня гран.условия не самые простые… Поэтому я и пытаюсь составить систему в общем виде…
Как начать по таким гран.условиям я долго не мог сообразить… Поэтому решил попробовать сделать это в общем виде… Для волнового уравнения там гораздо проще – там немного по другому выглядят гран.условия…

Я уже отчаялся полностью, не знаю как считать по такой схеме. Преподаватели говорят составляйте систему и решайте… Так вот я и пытаюсь составить… И не выходит… Другого пути я не вижу пока что

Такое уравнение возникло в результате выполнения курсовой… Я сижу с ним уже с весны… И сколько бы я не пытался, найти путной информации о том, как решить такое уравнение — мне не удалось… Поэтому пишу на разные форумы

Видимо Вы через чур обобщили. Как звучит исходная задача?

Вам именно исходную механическую формулировку? Или математическую?
Если механическую, то у нас есть плоская задача теории упругости в полупространстве. На плоскость z = 0 действует нагрузка p(x) (она может быть либо постоянной, либо как-то зависящей от координаты x) – нагрузка действует на конечном отрезке. Нужно определить напряжения и деформации при действии такой нагрузки на полуплоскость. Ну поскольку значение напряжений на бесконечности не задашь при численной реализации, то мы с преподавателем решили рассмотреть прямоугольную область [-a,a]х[0,b], b достаточно большое, больше, чем ширина участка действия нагрузки. Постановка этой задачи есть в книге Джонсона по механике контактного взаимодействия (стр 21). Цель работы – исследовать напряженно деформированное состояние в случае различного вида нагрузок и ширины участка нагружения…

Дальше была получена математическая постановка задачи – краевая задача для бигармонического уравнения в прямоугольнике с условиями такого вида. Уравнение было получено из основных соотношений теории упругости (условия совместности деформаций, закона Гука и т.п.) Для удобства была введена функция Эри, через которую и выражается напряжение и деформации…

Дальше нужно решить это уравнение численно. Я выбрал с помощью метода конечных разностей – получить систему уравнений на значение сеточной функции y_i^j в каждой точке сетки. Разностная схема представлена выше. У меня возникла идея, что можно по системе уравнений получить матрицу коэффициентов перед неизвестными и дальше любыми методами получить решение системы. И возникла проблема (даже не проблема скорее всего, а затруднение) – как автоматически заполнить эту матрицу. Т.е., если я знаю двойную индексацию переменной, получить ее позицию в столбце матрицы A[k][m] – k-номер уравнения, m – номер столбца.

Я полагаю, что если я знаю , какая у меня переменная (точнее ее индексацию) y[i[j], то в каждой k строке ее позиция есть Mj + i + 1…

Поэтому я и спросил, как все таки явно получить эту зависимость, да даже проще говорять, по системе уравнений получить матрицу системы относительно y[0][0], y[1][0] , … y[M][N]. Если эту систему записать на бумажке в случае сетки из 25 точек…

Если у вас будут идеи по поводу того, как заполнить эту матрицу – буду очень рад! Просто уже долго сижу с этой проблемой и пока оптимального решения не нашел…

Ниже представил код своих попыток – в некоторых местах позиция коэффициентов совпадает… Но не везде…

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][i + 1 + m * j+1] = 20 / pow(h_x, 4);

				A[k][i + 2 + m * j + 1] = -8 / pow(h_x, 4);
				A[k][i + 1 + m * (j + 1)] = -8 / pow(h_x, 4);
				A[k][i + 1 + m * (j - 1)] = -8 / pow(h_x, 4);

				A[k][i + 2 + m * (j+1)] = 2 / pow(h_x, 4);
				A[k][i  + m * (j+1)] = 2 / pow(h_x, 4);
				A[k][i  + m * (j-1)] = 2 / pow(h_x, 4);
				A[k][i + 2 + m * (j-1)] = 2 / pow(h_x, 4);

				A[k][i + 3 + m * j] = 1 / pow(h_x, 4);
				A[k][i + 1 + m * (j + 2)] = 1 / pow(h_x, 4);
				A[k][i - 1 + m * j] = 1 / pow(h_x, 4);
				A[k][i + 1 + m * (j - 2)] = 1 / pow(h_x, 4);

			}
		}

		
	}


	out_of_system(A);



	

	cout << "Формирование матрицы системы - 2 этап" << endl;

	int q0 = 1;
	int r0 = 2;

	for (int k = (m - 3) * (n - 3); k < (m - 3) * (n - 3) + m - 2; k++)
	{
		B[k] = -p(X[k], a, selector);		

				
			A[k][r0 + 1 +  q0] = 1 / pow(h_x, 2);
			A[k][r0 + 1 - 1 +  q0] = -2 / pow(h_x, 2);
			A[k][r0 + 1 - 2 + q0 ] = 1 / pow(h_x, 2);
			
		
		r0++;
		
		
	}


	out_of_system(A);
	
	int q1 = 1;
	int r1 = 1;

	cout << "Формирование матрицы системы - 3 этап" << endl;

	for (int k = (m - 3) * (n - 3) + m - 2; k < (m - 3) * (n - 3) + m - 2 + m - 1; k++)

	{
			
			A[k][r1 + 1 + m*0 + r1 * q1] = 1 / pow(h_x, 2);
			A[k][r1+1 + 1 + m*0 + r1 * q1] = -1 / pow(h_x, 2);
			A[k][r1 + 1 + m*1 + r1 * q1] = -1 / pow(h_x, 2);
			A[k][r1 + 1 + 1 + m*1 + r1 * q1] = 1 / pow(h_x, 2);

		
		r1++;

	}

	out_of_system(A);

	int q2 = 1;
	int r2 = 2;




	cout << "Формирование матрицы системы - 4 этап" << endl;

	for (int k = (m - 3) * (n - 3) + 2*m - 3 ; k < (m - 3) * (n - 3) + 3*m-5; k++)
	{
		
		
			A[k][r2-2 + 1 + n*n + r2*q2] = 1 / pow(h_x, 2);
			A[k][r2-1 + 1 + n*n + r2*q2] = -2 / pow(h_x, 2);
			A[k][r2 + 1 + n*n + r2*q2] = 1 / pow(h_x, 2);
	

		r2++;

	}

	out_of_system(A);

	int q3 = 1;
	int r3 = 0;
	

	cout << "Формирование матрицы системы - 5 этап" << endl;

	for (int k = (m - 3) * (n - 3) + 3 * m - 5; k < (m - 3) * (n - 3) + 4 * m - 5; k++)
	{
		A[k][r3 + 1 + 1 + n*n + r3*q3] = 1 / pow(h_x, 2);
		A[k][r3 + 1 + n*n + r3*q3] = -1 / pow(h_x, 2);
		A[k][r3 + 1+ 1+ (n-1)*(n-1) + r3*q3] = -1 / pow(h_x, 2);
		A[k][r3 +1 + (n-1)*(n-1) + r3*q3] = 1 / pow(h_x, 2);

	}

	out_of_system(A);

	/*cout << "Формирование матрицы системы - 6 этап" << endl;

	int q4 = 1;
	int r4 = 0;

	for (int k = (m - 3) * (n - 3) + 4 * m - 5; k < (m - 3) * (n - 3) + 4 * m - 5 + n-1; k++)
	{

		A[k][0 + (r4+1)*n + r4*q4] = 1 / pow(h_x, 2);
		A[k][0 + r4 * n + r4 * q4] = -2 / pow(h_x, 2);
		A[k][0 + (r4 - 1) * n + r4 * q4] = 1 / pow(h_x, 2);
	}

	int q5 = 1;
	int r5 = 0;

	for (int k = (m - 3) * (n - 3) + 4 * m - 5 + n - 1; k < (m - 3) * (n - 3) + 4 * m - 5 + 2*(n - 1); k++)
	{
		A[k][1 + (r5+1)*n + r5*q5] = 1 / pow(h_x, 2);
		A[k][0 + (r5 + 1) * n + r5 * q5] = -1 / pow(h_x, 2);
		A[k][1 + (r5 + 0) * n + r5 * q5] = -1 / pow(h_x, 2);
		A[k][0 + (r5 + 0) * n + r5 * q5] = 1 / pow(h_x, 2);
	}

	int q6 = 1;
	int r6 = 0;


	for (int k = (m - 3) * (n - 3) + 4 * m - 5 + 2 * (n - 1); k < (m - 3) * (n - 3) + 4 * m - 5 + 3 * (n - 1); k++)
	{
		A[k][m + (r6 + 1) * n + r6 * q6] = 1 / pow(h_x, 2);
		A[k][m + r6 * n + r6 * q6] = -2 / pow(h_x, 2);
		A[k][m + (r6 - 1) * n + r6 * q6] = 1 / pow(h_x, 2);
	}

	int q7 = 1;
	int r7 = 0;

	for (int k = (m - 3) * (n - 3) + 4 * m - 5 + 3 * (n - 1); k < (m - 3) * (n - 3) + 4 * m - 5 + 4 * (n - 1); k++)
	{
		A[k][m + (r7 + 1) * n + r7 * q7] = 1 / pow(h_x, 2);
		A[k][m-1 + (r7 + 1) * n + r7 * q7] = -1 / pow(h_x, 2);
		A[k][m + (r7 + 0) * n + r7 * q7] = -1 / pow(h_x, 2);
		A[k][m-1 + (r7 + 0) * n + r7 * q7] = 1 / pow(h_x, 2);
	}

Скорее всего я наверное перемудрил со способом решения… Но пока это единственное, что я придумал. Преподаватели, с кем я консультировался – говорят делать то же самое – тупо записать систему уравнений…

Почему именно Си, не MATLAB?

Два дня писал Вам ответ. Всё зачеркнул.

Постановка этой задачи есть в книге Джонсона по механике контактного взаимодействия (стр 21).

Там неправильная постановка от слов совсем. И книгу следует выкинуть в мусорное ведро.

Короче берём книгу там
Зенкевич О. Методы конечных элементов в технике 1975
Там дана верная постановка.

Поэтому я и спросил, как все таки явно получить эту зависимость, да даже проще говорять, по системе уравнений получить матрицу системы относительно

Что касается вычисления коэффициентов. То лучше всего это описано в книге
Matthew N.O. Sadiku - Numerical Techniques in Electromagnetics-CRC Press (2000)

1 Симпатия

О как, это уже становится более понятным. Так понимаю, на современный лад прозвучит как: Конечные автоматы.
Или ошибаюсь?

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

1 Симпатия

Ошибаетесь. Метод конечных элементов(МКЭ) если перевести на английский Finite Element Method (FEM).
А конечные автоматы это FSM.
Метод МКЭ предназначен для решения дифференциальных уравнений возникающий в технике и физике.

К примеру нужно решить задачу о разрушении цельно сварного стула или напечатанного на 3D-принтере. Условие разрушение предельное растяжение элементов стула. Тут сила приложена к стулу показана стелкой.

Прежде чем решать такие сложные задачи педологию рассмотреть более простую …

Вот тут метод не плохо описан.
http://bibl.nngasu.ru/electronicresources/uch-metod/physics/847226.pdf

1 Симпатия

Смог решить проблему через Wolfram Mathematica и сформировать матрицу коэффициентов