double E[8][8]={0}, TMP1[8][8], TMP2[8][8];

int i,j,k;

//E - единичная матрица - первый член ряда экспоненты

E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;

E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;

//первоначальное заполнение вспомогательного массива TMP1 - предыдущего члена ряда для следующего перемножения

//и первоначальное заполнение экспоненты первым членом ряда

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=E[i][j];

EXP[i][j]=E[i][j];

}

}

//ряд вычисления экспоненты EXP, начиная со 2-го члена ряда (m=2;m<=n)

for(m=2;m<=n;m++) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP2[i][j]=0;

for(k=0;k<8;k++) {

//TMP2[i][j]+=TMP1[i][k]*A[k][j]*delta_x/(m-1);

TMP2[i][j]+=TMP1[i][k]*A[k][j];

}

TMP2[i][j]*=delta_x;//вынесено за цикл произведения строки на столбец

TMP2[i][j]/=(m-1);//вынесено за цикл произведения строки на столбец

EXP[i][j]+=TMP2[i][j];

}

}

//заполнение вспомогательного массива TMP1 для вычисления следующего члена ряда - TMP2 в следующем шаге цикла по m

if (m<n) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=TMP2[i][j];

}

}

}

}

}

//Вычисление матрицы MAT_ROW в виде матричного ряда для последующего использования

//при вычислении вектора partial_vector - вектора частного решения неоднородной системы ОДУ на шаге delta_x

void mat_row_for_partial_vector(double A[8][8], double delta_x, double MAT_ROW[8][8]) {

//n - количество членов ряда в MAT_ROW, m - счетчик членов ряда (m<=n)

int n=100, m;

double E[8][8]={0}, TMP1[8][8], TMP2[8][8];

int i,j,k;

//E - единичная матрица - первый член ряда MAT_ROW

E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;

E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;

//первоначальное заполнение вспомогательного массива TMP1 - предыдущего члена ряда для следующего перемножения

//и первоначальное заполнение MAT_ROW первым членом ряда

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=E[i][j];

MAT_ROW[i][j]=E[i][j];

}

}

//ряд вычисления MAT_ROW, начиная со 2-го члена ряда (m=2;m<=n)

for(m=2;m<=n;m++) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP2[i][j]=0;

for(k=0;k<8;k++) {

TMP2[i][j]+=TMP1[i][k]*A[k][j];

}

TMP2[i][j]*=delta_x;

TMP2[i][j]/=m;

MAT_ROW[i][j]+=TMP2[i][j];

}

}

//заполнение вспомогательного массива TMP1 для вычисления следующего члена ряда - TMP2 в следующем шаге цикла по m

if (m<n) {

for(i=0;i<8;i++) {

for(j=0;j<8;j++) {

TMP1[i][j]=TMP2[i][j];

}

}

}

}

}

//Вычисление vector - НУЛЕВОГО (частный случай) вектора частного решения

//неоднородной системы дифференциальных уравнений на рассматриваемом участке:

void partial_vector(double vector[8]){

for(int i=0;i<8;i++){

vector[i]=0.0;

}

}

//Вычисление vector - вектора частного решения неоднородной системы дифференциальных уравнений на рассматриваемом участке delta_x:

void partial_vector_real(double expo_[8][8], double mat_row[8][8], double x_, double delta_x, double vector[8]){

double POWER_[8]={0};//Вектор внешней нагрузки на оболочку

double REZ[8]={0};

double REZ_2[8]={0};

power_vector_for_partial_vector(x_, POWER_);//Расчитываем POWER_ при координате x_

mat_on_vect(mat_row, POWER_, REZ);//Умножение матрицы mat_row на вектор POWER_ и получаем вектор REZ

mat_on_vect(expo_, REZ, REZ_2);//Умножение матрицы expo_ на вектор REZ и получаем вектор REZ_2

for(int i=0;i<8;i++){

vector[i]=REZ_2[i]*delta_x;

}

}

//Решение СЛАУ размерности 8 методом Гаусса с выделением главного элемента

int GAUSS(double AA[8][8], double bb[8], double x[8]){

double A[8][8];

double b[8];

for(int i=0;i<8;i++){

b[i]=bb[i];//Работать будем с вектором правых частей b, чтобы исходный вектор bb не изменялся при выходе из подпрограммы

for(int j=0;j<8;j++){

A[i][j]=AA[i][j];//Работать будем с матрицей А, чтобы исходная матрица АА не менялась при выходе из подпрограммы

}

}

int e;//номер строки, где обнаруживается главный (максимальный) коэффициент в столбце jj

double s, t, main;//Вспомогательная величина

for(int jj=0;jj<(8-1);jj++){//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

e=-1; s=0.0; main=A[jj][jj];

for(int i=jj;i<8;i++){//Находится номер е строки, где лежит главный (максимальный) элемент в столбце jj и делается взаимозамена строк

if ((A[i][jj]*A[i][jj])>s) {//Вместо перемножения (удаляется возможный знак минуса) можно было бы использовать функцию по модулю abs()

e=i; s=A[i][jj]*A[i][jj];

}

}

if (e<0) {

cout<<"Mistake "<<jj<<"\n"; return 0;

}

if (e>jj) {//Если главный элемент не в строке с номером jj. а в строке с номером е

main=A[e][jj];

for(int j=0;j<8;j++){//Взаимная замена двух строк - с номерами e и jj

t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t;

}

t=b[jj]; b[jj]=b[e]; b[e]=t;

}

for(int i=(jj+1);i<8;i++){//Приведение к верхнетреугольной матрице

for(int j=(jj+1);j<8;j++){

A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//Перерасчет коэффициентов строки i>(jj+1)

}

b[i]=b[i]-(1/main)*b[jj]*A[i][jj];

A[i][jj]=0.0;//Обнуляемые элементы столбца под диагональным элементом матрицы А

}

}//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

x[8-1]=b[8-1]/A[8-1][8-1];//Первоначальное определение последнего элемента искомого решения х (7-го)

for(int i=(8-2);i>=0;i--){//Вычисление елементов решения x[i] от 6-го до 0-го

t=0;

for(int j=1;j<(8-i);j++){

t=t+A[i][i+j]*x[i+j];

}

x[i]=(1/A[i][i])*(b[i]-t);

}

return 0;

}

int main()

{

int nn;//Номер гармоники, начиная с 1-й (без нулевой)

int nn_last=50;//Номер последней гармоники

double Moment[100+1]={0};//Массив физического параметра (момента), что рассчитывается в каждой точке между краями

double angle;

double start_angle, finish_angle;

start_angle=3.14159265359/4;

finish_angle=start_angle+(3.14159265359/2);

double step=(3.14159265359/2)/100; //step=(3.14159265359/2)/100 - величина шага расчета оболочки - шага интервала интегрирования (должна быть больше нуля, т.е. положительная)

double h_div_R;//Величина h/R

h_div_R=1.0/200;

double c2;

c2=h_div_R*h_div_R/12;//Величина h*h/R/R/12

double nju;

nju=0.3;

double gamma;

gamma=3.14159265359/4;//Угол распределения силы по левому краю

//распечатка в файлы:

FILE *fp;

// Open for write

if( (fp = fopen( "C:/test.txt", "w" )) == NULL ) // C4996