не могу найти ошибку в коде С++ - Помощь студентам

kirillkucelap · 06.04.2014, 18:56

я писал программу по методу конечных элементов, в этом методе по исходным формируется матрица gg и столбец ff некоторой размерности m*m и m соответственно, далее к ним применяется метод гаусса (решение системы линейных уравнений), сам гаусс начинается на 167 строчке кода. я вводил в качестве исходных данных ввожу последовательно числа: 3-1-1-1-1-1-1-1-2-1-1-1-1-2-2-1-1-1-1-1-1-2-3. В результате правильно формируется система уравнений, но она неправильно решается. В комментарии на 171 строчке приведены те же правильные значения элементов массивов gg и ff, и если убрать //- символ комментариев, то ту же систему уравнений программа решает правильно!!!! объясните мне в чем дело. Код вынужден разбить на несколько частей (в одном сообщении не более 5000 символов)

Код:

#include <iostream>
#include<math.h>
#include<clocale>
using namespace std;
const int sizeuzel = 51;
const int sizeelement = 51;
struct uzels{
long double x;// координата Х для данного узла
long double y;// координата Y
long double teploprovod;//величина теплопроводности лямда
long double teplootdacha;//величина теплоотдачи альфа
long double qv;
long double qs;
} uzel[sizeuzel];
struct elements{
	int first; 
	int second;
	int third;
	long double S;// площадь треугольника
long double teploprovod;//величина теплопроводности лямда
long double teplootdacha;//величина теплоотдачи альфа
long double qv;
long double qs;
long double a1, a2, a3, b1, b2, b3, c1, c2, c3;
} element[sizeelement];
int M; //число узлов
int N; //число элементов
int q, w, e, r,u;//для циклов
long double L[sizeuzel][sizeuzel];//расстояние между узлами
long double g [sizeelement][4][4];
long double f [sizeelement][4];
long double gg[sizeuzel][sizeuzel];//глобальные матрицы
long double ff[sizeuzel];
long double pp[sizeuzel];
int p,kk;
int T[sizeuzel];// для метода гаусса, отслеживает перестановки в матрице
long double z,k;
int main ()
{std:setlocale(LC_ALL,"Russian_Russia.1251");
cout<<"Максимальное число узлов "<<sizeuzel-1<<"\n";
cout<<"Максимальное число элементов "<<sizeelement-1<<"\n";
mmm:
cout<<"Введите число узлов";//ввод числа узлов
cin>>M;
if (M>sizeuzel-1){ 
	cout<<"ошибка";
	goto mmm;
}//предотвращение ощибки большого числа узлов(доделать)
cout<<"Введите координаты и параметры узлов"<<"\n";// введение координат узлов
for (q=1;q<M+1;q++){
cout<<"Введите координату Х узла"<<q<<"\n";
cin>>uzel[q].x;
cout<<"\n"<<"Ведите координату Y узла"<<q<<"\n";
cin>>uzel[q].y;
cout<<"\n"<<"Введите теплопроводность в узле"<<q<<"\n";// введение физических параметров в узлах
cin>>uzel[q].teploprovod;
cout<<"\n"<<"Введите теплоотдачу в узле"<<q<<"\n";
cin>>uzel[q].teplootdacha;
cout<<"\n"<<"Введите qv в узле "<<q<<"\n";
cin>>uzel[q].qv;
cout<<"\n"<<"Введите qs в узле "<<q<<"\n";
cin>>uzel[q].qs;
}
q=0;
mm:
cout<<"Ведите число элементов";//ввод числа элементов
cin>>N;
if (N>sizeelement-1) {
	cout<<"ошибка";
	goto mm;
}//предотвращение ощибки большого числа элементов(доделать)
for (q=1;q<N+1;q++){// ввод номеров узлов в элементы
	cout<<"Введите номер первого узла элемента "<<q<<"\n";
	cin>>element[q].first;
	cout<<"Введите номер второго узла элемента "<<q<<"\n";
	cin>>element[q].second;
	cout<<"Введите номер третьего узла элемента "<<q<<"\n";
	cin>>element[q].third;

kirillkucelap · 06.04.2014, 18:58

Код:

element[q].teploprovod=(uzel[element[q].first].teploprovod+uzel[element[q].second].teploprovod+uzel[element[q].third].teploprovod)/3;
    element[q].teplootdacha=(uzel[element[q].first].teplootdacha+uzel[element[q].second].teplootdacha+uzel[element[q].third].teplootdacha)/3;
	element[q].qv=(uzel[element[q].first].qv+uzel[element[q].second].qv+uzel[element[q].third].qv)/3;
	element[q].qs=(uzel[element[q].first].qs+uzel[element[q].second].qs+uzel[element[q].third].qs)/3;//далее высчитывается площадь теугольников
	element[q].S=fabs((uzel[element[q].second].x*uzel[element[q].third].y-uzel[element[q].third].x*uzel[element[q].second].y+uzel[element[q].first].x*uzel[element[q].second].y-uzel[element[q].first].x*uzel[element[q].third].y+uzel[element[q].third].x*uzel[element[q].first].y-uzel[element[q].second].x*uzel[element[q].first].y))/2;
    element[q].a1=(uzel[element[q].second].x*uzel[element[q].third].y-uzel[element[q].third].x*uzel[element[q].second].y)/(2*element[q].S);
	element[q].a2=(uzel[element[q].third].x*uzel[element[q].first].y-uzel[element[q].first].x*uzel[element[q].third].y)/(2*element[q].S);
	element[q].a3=(uzel[element[q].first].x*uzel[element[q].second].y-uzel[element[q].second].x*uzel[element[q].first].y)/(2*element[q].S);
	element[q].b1=(uzel[element[q].second].y-uzel[element[q].third].y)/(2*element[q].S);
	element[q].b2=(uzel[element[q].third].y-uzel[element[q].first].y)/(2*element[q].S);
	element[q].b3=(uzel[element[q].first].y-uzel[element[q].second].y)/(2*element[q].S);
	element[q].c1=(uzel[element[q].third].x-uzel[element[q].second].x)/(2*element[q].S);
	element[q].c2=(uzel[element[q].first].x-uzel[element[q].third].x)/(2*element[q].S);
	element[q].c3=(uzel[element[q].second].x-uzel[element[q].first].x)/(2*element[q].S);
}
q=0;
for(e=1;e<M+1;e++){//определение внешних ребер и расчет их длин
	for(w=1;w<M+1;w++){
		L[e][w]=0;
	}
}
cout<<"Проверьте список внешних ребер:";
for(e=1;e<M+1;e++){
	for(w=e+1;w<M+1;w++){
		r=0;
		for(q=1;q<N+1;q++){
			if(((e==element[q].first)||(e==element[q].second)||(e==element[q].third))&&((w==element[q].first)||(w==element[q].second)||(w==element[q].third))) r=r+1;
		}
		if(r>2) cout<<"ошибка - 1 ребро принадлежит 3 элементам";
		if(r==1) {
		L[e][w]=sqrt((uzel[e].x-uzel[w].x)*(uzel[e].x-uzel[w].x)+(uzel[e].y-uzel[w].y)*(uzel[e].y-uzel[w].y));
		L[w][e]=sqrt((uzel[e].x-uzel[w].x)*(uzel[e].x-uzel[w].x)+(uzel[e].y-uzel[w].y)*(uzel[e].y-uzel[w].y));
		cout<<e<<" "<<w<<" его длина ="<<L[e][w]<<"\n";
		}
	}
}
q=0;e=0;w=0;// обнуление локальных массивов g и f
for (q=0;q<sizeelement;q++){
	for(w=0;w<4;w++){
		f[q][w]=0;
    	for(e=0;e<4;e++) g[q][w][e]=0;
	}
}
q=0;w=0;e=0;

kirillkucelap · 06.04.2014, 19:00

Код:

for (q=1;q<N+1;q++){// составление локальных матриц f и g
  g[q][1][1]=(element[q].b1*element[q].b1+element[q].c1*element[q].c1)*element[q].teploprovod*element[q].S+(((uzel[element[q].first].teplootdacha+uzel[element[q].second].teplootdacha)/2)*L[element[q].first][element[q].second]+((uzel[element[q].first].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].first][element[q].third])/3;
  g[q][1][2]=(element[q].b1*element[q].b2+element[q].c1*element[q].c2)*element[q].teploprovod*element[q].S+((uzel[element[q].first].teplootdacha+uzel[element[q].second].teplootdacha)/2)*L[element[q].first][element[q].second]/6;
  g[q][1][3]=(element[q].b1*element[q].b3+element[q].c1*element[q].c3)*element[q].teploprovod*element[q].S+((uzel[element[q].first].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].first][element[q].third]/6;
  g[q][2][1]=(element[q].b1*element[q].b2+element[q].c1*element[q].c2)*element[q].teploprovod*element[q].S+((uzel[element[q].first].teplootdacha+uzel[element[q].second].teplootdacha)/2)*L[element[q].first][element[q].second]/6;
  g[q][2][2]=(element[q].b2*element[q].b2+element[q].c2*element[q].c2)*element[q].teploprovod*element[q].S+(((uzel[element[q].first].teplootdacha+uzel[element[q].second].teplootdacha)/2)*L[element[q].first][element[q].second]+((uzel[element[q].second].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].second][element[q].third])/3;
  g[q][2][3]=(element[q].b2*element[q].b3+element[q].c2*element[q].c3)*element[q].teploprovod*element[q].S+((uzel[element[q].second].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].second][element[q].third]/6;
  g[q][3][1]=(element[q].b1*element[q].b3+element[q].c1*element[q].c3)*element[q].teploprovod*element[q].S+((uzel[element[q].first].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].first][element[q].third]/6;
  g[q][3][2]=(element[q].b2*element[q].b3+element[q].c2*element[q].c3)*element[q].teploprovod*element[q].S+((uzel[element[q].second].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].second][element[q].third]/6;
  g[q][3][3]=(element[q].b3*element[q].b3+element[q].c3*element[q].c3)*element[q].teploprovod*element[q].S+(((uzel[element[q].first].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].first][element[q].third]+((uzel[element[q].second].teplootdacha+uzel[element[q].third].teplootdacha)/2)*L[element[q].second][element[q].third])/3;
  f[q][1]=element[q].qv*element[q].S/3+(((uzel[element[q].first].qs+uzel[element[q].second].qs)/2)*L[element[q].first][element[q].second]+((uzel[element[q].first].qs+uzel[element[q].third].qs)/2)*L[element[q].first][element[q].third])/2;
  f[q][2]=element[q].qv*element[q].S/3+(((uzel[element[q].first].qs+uzel[element[q].second].qs)/2)*L[element[q].first][element[q].second]+((uzel[element[q].second].qs+uzel[element[q].third].qs)/2)*L[element[q].second][element[q].third])/2;
  f[q][3]=element[q].qv*element[q].S/3+(((uzel[element[q].first].qs+uzel[element[q].third].qs)/2)*L[element[q].first][element[q].third]+((uzel[element[q].second].qs+uzel[element[q].third].qs)/2)*L[element[q].second][element[q].third])/2;
}
q=0;w=0;e=0;//обнуление глобальных матриц 
for(q=0;q<sizeuzel;q++){
	ff[q]=0;
	for(w=0;w<sizeuzel;w++){
		gg[q][w]=0;
	}
}
q=0;w=0;e=0;
for(q=1;q<N+1;q++){
   gg[element[q].first][element[q].first]=gg[element[q].first][element[q].first]+g[q][1][1];
   gg[element[q].first][element[q].second]=gg[element[q].first][element[q].second]+g[q][1][2];
   gg[element[q].first][element[q].third]=gg[element[q].first][element[q].third]+g[q][1][3];
   gg[element[q].second][element[q].first]=gg[element[q].second][element[q].first]+g[q][2][1];
   gg[element[q].second][element[q].second]=gg[element[q].second][element[q].second]+g[q][2][2];
   gg[element[q].second][element[q].third]=gg[element[q].second][element[q].third]+g[q][2][3];
   gg[element[q].third][element[q].first]=gg[element[q].third][element[q].first]+g[q][3][1];
   gg[element[q].third][element[q].second]=gg[element[q].third][element[q].second]+g[q][3][2];
   gg[element[q].third][element[q].third]=gg[element[q].third][element[q].third]+g[q][3][3];
   ff[element[q].first]=ff[element[q].first]+f[q][element[q].first];
   ff[element[q].second]=ff[element[q].second]+f[q][element[q].second];
   ff[element[q].third]=ff[element[q].third]+f[q][element[q].third];
}
q=0;w=0;e=0;
cout<<"система уравнений:"<<"\n";
for(q=1;q<M+1;q++){
	for(w=1;w<M+1;w++){
		cout<<gg[q][w]<<"*"<<"u"<<w<<"+";
	}
	cout<<"="<<ff[q]<<"\n"; 
}

kirillkucelap · 06.04.2014, 19:01

Код:

q=0;w=0;e=0;//метод ГАУССА
cout<<gg[1][1]<<" "<<gg[1][2]<<" "<<gg[1][3]<<" "<<ff[1]<<"\n";
cout<<gg[2][1]<<" "<<gg[2][2]<<" "<<gg[2][3]<<" "<<ff[2]<<"\n";
cout<<gg[3][1]<<" "<<gg[3][2]<<" "<<gg[3][3]<<" "<<ff[3]<<"\n"<<"\n";
//gg[1][1]=1.30474;gg[1][2]=-0.333333;gg[1][3]=0.235702;gg[2][1]=-0.333333;gg[2][2]=1.66667;gg[2][3]=-0.333333;gg[3][1]=0.235702;gg[3][2]=-0.333333;gg[3][3]=1.30474;ff[1]=1.37377;ff[2]=1.66667;ff[3]=1.37377;
cout<<gg[1][1]<<" "<<gg[1][2]<<" "<<gg[1][3]<<" "<<ff[1]<<"\n";
cout<<gg[2][1]<<" "<<gg[2][2]<<" "<<gg[2][3]<<" "<<ff[2]<<"\n";
cout<<gg[3][1]<<" "<<gg[3][2]<<" "<<gg[3][3]<<" "<<ff[3]<<"\n"<<"\n";
for(q=0;q<M+1;q++){//иницианализация вектора R
	T[q]=q;
}
q=0;w=0;e=0;u=0;
for(q=1;q<M+1;q++) {// цикл по уравнениям в методе гаусса прямой ход
	e=0;
	z=0;
	for(w=q;w<M+1;w++){//нахождение максимального элемента в строке
    if(fabs(gg[q][w])>fabs(z)) {
		z=gg[q][w];
		e=w;
	}
	}
	kk=T[q];
	T[q]=T[e];
	T[e]=kk;
	
	for(u=1;u<M+1;u++){//перестановка столбцов
		pp[u]=gg[u][q];//сохраняем первый (а потом второй и т д ) столбец
		gg[u][q]=gg[u][e];
		gg[u][e]=pp[u];
        }
	cout<<gg[1][1]<<" "<<gg[1][2]<<" "<<gg[1][3]<<" "<<ff[1]<<"\n";
cout<<gg[2][1]<<" "<<gg[2][2]<<" "<<gg[2][3]<<" "<<ff[2]<<"\n";
cout<<gg[3][1]<<" "<<gg[3][2]<<" "<<gg[3][3]<<" "<<ff[3]<<"\n"<<"\n";
	w=0;
	k=gg[q][q];
	ff[q]=ff[q]/k;
	for(w=q;w<M+1;w++){
		gg[q][w]=gg[q][w]/k;//переписываем кутую строчку
	}
	cout<<gg[1][1]<<" "<<gg[1][2]<<" "<<gg[1][3]<<" "<<ff[1]<<"\n";
cout<<gg[2][1]<<" "<<gg[2][2]<<" "<<gg[2][3]<<" "<<ff[2]<<"\n";
cout<<gg[3][1]<<" "<<gg[3][2]<<" "<<gg[3][3]<<" "<<ff[3]<<"\n"<<"\n";
	
	w=0;e=0;
	for(e=q+1;e<M+1;e++){//переписываем другие строчки
		k=gg[e][q];
		ff[e]=ff[e]-k*ff[q];
		for(w=q;w<M+1;w++){
			
		gg[e][w]=gg[e][w]-k*gg[q][w];
		}
	}
	cout<<gg[1][1]<<" "<<gg[1][2]<<" "<<gg[1][3]<<" "<<ff[1]<<"\n";
cout<<gg[2][1]<<" "<<gg[2][2]<<" "<<gg[2][3]<<" "<<ff[2]<<"\n";
cout<<gg[3][1]<<" "<<gg[3][2]<<" "<<gg[3][3]<<" "<<ff[3]<<"\n"<<"\n";
}
q=0;w=0;e=0;u=0;
for(q=M;q>0;q--){    //обратный ход
	for(w=1;w<q;w++){
		ff[w]=ff[w]-ff[q]*gg[w][q];
	}
    }
cout<<"\n"<<"температуры в узлах равны: ";
q=0;w=0;
for(q=1;q<M+1;q++){
	cout<<"u"<<q<<"=";
	for(w=1;w<M+1;w++){
		if(T[w]==q) {
			cout<<ff[w]<<"\n";
		}
	}
}



cout<<"Для выхода введите любое число";
cin>>p;
return 0;  
}

kirillkucelap · 06.04.2014, 19:03

собственно метод гаусса приведен в 4 сообщении, по логике если в 5 строчке 4 сообщения убрать //, то ответ измениться не должен, но он меняется на правильный. В чем же дело ?

Опции темы	Поиск в этой теме
Версия для печати Отправить по электронной почте	Поиск в этой теме: Расширенный поиск

Похожие темы
Тема	Автор	Раздел	Ответов	Последнее сообщение
Не могу найти ошибку в коде	youngster	Паскаль, Turbo Pascal, PascalABC.NET	3	17.11.2013 22:26
Не могу найти ошибку в коде	afirat	Общие вопросы Delphi	7	28.12.2012 21:13
Не могу найти ошибку в коде	sashmedv	Паскаль, Turbo Pascal, PascalABC.NET	2	25.03.2012 08:31
не могу найти ошибку в коде	pavelstraut	Общие вопросы C/C++	5	24.07.2009 23:20
Не могу найти ошибку в коде!	Natasha666	Помощь студентам	1	20.05.2009 09:27

06.04.2014, 19:03	#5
kirillkucelap Регистрация: 06.04.2014 Сообщений: 8	собственно метод гаусса приведен в 4 сообщении, по логике если в 5 строчке 4 сообщения убрать //, то ответ измениться не должен, но он меняется на правильный. В чем же дело ?