Методы решения жестких и нежестких краевых задач

void orto_norm_4x8(double A[4][8], double b[4], double C[4][8], double d[4]){

double NORM;

double mult0,mult1,mult2,mult3,mult4,mult5,mult6,mult7;

//Получаем 1-ю строку уравнения C*x=d:

NORM=norma(A,0);

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

C[0][k]=A[0][k]/NORM;

}

d[0]=b[0]/NORM;

//Получаем 2-ю строку уравнения C*x=d:

mult0=mult(A,1,C,0);

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

C[1][k]=A[1][k]-mult0*C[0][k];

}

NORM=norma(C,1);

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

C[1][k]/=NORM;

}

d[1]=(b[1]-mult0*d[0])/NORM;

//Получаем 3-ю строку уравнения C*x=d:

mult0=mult(A,2,C,0); mult1=mult(A,2,C,1);

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

C[2][k]=A[2][k]-mult0*C[0][k]-mult1*C[1][k];

}

NORM=norma(C,2);

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

C[2][k]/=NORM;

}

d[2=(b[2]-mult0*d[0]-mult1*d[1])/NORM;

//Получаем 4-ю строку уравнения C*x=d:

mult0=mult(A,3,C,0); mult1=mult(A,3,C,1); mult2=mult(A,3,C,2);

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

C[3][k]=A[3][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k];

}

NORM=norma(C,3);

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

C[3][k]/=NORM;

}

d[3]=(b[3]-mult0*d[0]-mult1*d[1]-mult2*d[2])/NORM;

}

//Произведение матрицы A1 размерности 4х8 на матрицу А2 размерности 8х8. Получаем матрицу rezult размерности 4х8:

void mat_4x8_on_mat_8x8(double A1[4][8], double A2[8][8], double rezult[4][8]){

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

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

rezult[i][j]=0.0;

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

rezult[i][j]+=A1[i][k]*A2[k][j];

}

}

}

}

//Умножение матрицы A на вектор b и получаем rezult.

void mat_on_vect(double A[8][8], double b[8], double rezult[8]){

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

rezult[i]=0.0;

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

rezult[i]+=A[i][k]*b[k];

}

}

}

//Умножение матрицы A размерности 4х8 на вектор b размерности 8 и получаем rezult размерности 4.

void mat_4x8_on_vect_8(double A[4][8], double b[8], double rezult[4]){

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

rezult[i]=0.0;

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

rezult[i]+=A[i][k]*b[k];

}

}

}

//Вычитание из вектора u1 вектора u2 и получение вектора rez=u1-u2. Все вектора размерности 4.

void minus(double u1[4], double u2[4], double rez[4]){

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

rez[i]=u1[i]-u2[i];

}

}

//Вычисление матричной экспоненты EXP=exp(A*delta_x)

void exponent(double A[8][8], double delta_x, double EXP[8][8]) {

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

int n=100, m;

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

printf( "The file 'C:/test.txt' was not opened\n" );

else

printf( "The file 'C:/test.txt' was opened\n" );

for(nn=1;nn<=nn_last;nn++){ //ЦИКЛ ПО ГАРМОНИКАМ, НАЧИНАЯ С 1-ОЙ ГАРМОНИКИ (БЕЗ НУЛЕВОЙ ГАРМОНИКИ)

double expo_from_minus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (0-x1)

double expo_from_plus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (x1-0)

double mat_row_for_minus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (0-x1)

double mat_row_for_plus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (x1-0)

double MATRIXS[100+1][8][8]={0};//Массив из матриц размерности 8x8 для решения СЛАУ в каждой точке интервала интегрирования

double VECTORS[100+1][8]={0};//Массив векторов правых частей размерности 8 соответствующих СЛАУ

double U[4][8]={0};//Матрица краевых условий левого края размерности 4х8

double u_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий левого края

double V[4][8]={0};//Матрица краевых условий правого края размерности 4х8

double v_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий правого края

double Y[100+1][8]={0};//Массив векторов-решений соответствующих СЛАУ (в каждой точке интервала между краями): MATRIXS*Y=VECTORS

double A[8][8]={0};//Матрица коэффициентов системы ОДУ

double FF[8]={0};//Вектор частного решения неоднородной ОДУ на участке интервала интегрирования

double Ui[4][8]={0};//Вспомогательная матрица коэффициентов переносимых краевых условий с левого края

double ui_[4]={0};//Правые части переносимых краевых условий с левого края

double ui_2[4]={0};//вспомогательный вектор (промежуточный)

double UiORTO[4][8]={0};//Ортонормированная переносимая матрица с левого края

double ui_ORTO[4]={0};//Соответственно правые части ортонормированного переносимого уравнения с левого края

double Vj[4][8]={0};//Вспомогательная матрица коэффициентов переносимых краевых условий с правого края

double vj_[4]={0};//Правые части переносимых краевых условий с правого края

double vj_2[4]={0};//Вспомогательный вектор (промежуточный)

double VjORTO[4][8]={0};//Ортонормированная переносимая матрица с правого края

double vj_ORTO[4]={0};//Соответственно правые части ортонормированного переносимого уравнения с правого края

double MATRIX_2[8][8]={0};//Вспомогательная матрица

double VECTOR_2[8]={0};//Вспомогательный вектор

double Y_2[8]={0};//Вспомогательный вектор

double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Возведенный в соответствующие степени номер гармоники nn

nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4;

//Здесь надо первоначально заполнить ненулевыми значениями матрицы и вектора краевых условий U*Y[0]=u_ (слева) и V*Y[100]=v_ (справа) :

U[0][0]=nju/tan(start_angle);

U[0][1]=1.0;

U[0][2]=nju*nn/sin(start_angle);

U[0][4]=(1+nju);

u_[0]=0.0;//Сила T1 на левом крае равна нулю

U[1][0]=-(1-nju)/2/sin(start_angle);

U[1][2]=-(1-nju)/2/tan(start_angle);

U[1][3]=(1-nju)/2;

U[1][4]=-c2*nn*(1-nju)/sin(start_angle)/tan(start_angle);

U[1][5]=c2*nn*(1-nju)/sin(start_angle);

u_[1]=0.0;//Сила S* на левом краю равна нулю

U[2][4]=-nju*nn2/sin(start_angle)/sin(start_angle);

U[2][5]=nju/tan(start_angle);

U[2][6]=1.0;

u_[2]=0;//Момент M1 на левом краю равен нулю

U[3][4]=-(3-nju)*nn2/sin(start_angle)/sin(start_angle)/tan(start_angle);

U[3][5]=nju+1.0/tan(start_angle)/tan(start_angle)-(nju-2)*nn2/sin(start_angle)/sin(start_angle);

U[3][6]=-1.0/tan(start_angle);

U[3][7]=-1.0;

u_[3]=-sin(nn*gamma)/(nn*gamma);//Сила Q1* на левом крае распределена на угол -gamma +gamma

V[0][0]=1.0; v_[0]=0.0;//Перемещение u на правом крае равно нулю

V[1][2]=1.0; v_[1]=0.0;//Перемещение v на правом крае равно нулю

V[2][4]=1.0; v_[2]=0.0;//Перемещение w на правом крае равно нулю

V[3][5]=1.0; v_[3]=0.0;//Угол поворота на правом крае равен нулю

//Здесь заканчивается первоначальное заполнение U*Y[0]=u_ и V*Y[100]=v_

orto_norm_4x8(U, u_, UiORTO, ui_ORTO);//Первоначальное ортонормирование краевых условий

orto_norm_4x8(V, v_, VjORTO, vj_ORTO);

//Первоначальное заполнение MATRIXS и VECTORS матричными уравнениями краевых условий соответственно

//UiORTO*Y[0]=ui_ORTO и VjORTO*Y[100]=vj_ORTO:

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

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

MATRIXS[0][i][j]=UiORTO[i][j];//Левый край; верхнее матричное уравнение

MATRIXS[100][i+4][j]=VjORTO[i][j];//Правый край (точка номер 101 с индексом 100 - отсчет идет с нуля); нижнее матричное уравнение

}

VECTORS[0][i]=ui_ORTO[i];//Левый край; верхнее матричное уравнение

VECTORS[100][i+4]=vj_ORTO[i];//Правый край (точка номер 101 с индексом 100 - отсчет идет с нуля); нижнее матричное уравнение

}

//Цикл по точкам ii интервала интегрирования заполнения ВЕРХНИХ частей матричных уравнений MATRIXS[ii]*Y[ii]=VECTORS[ii],

//начиная со второй точки - точки с индексом ii=1

angle=start_angle;//начальное значение угловой координаты

for(int ii=1;ii<=100;ii++){

angle+=step;//Угловая координата

A_perem_coef(nju, c2, nn, angle, A);//Вычисление матрицы А коэффициентов системы ОДУ при данной угловой координате angle

exponent(A,(-step),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты)

mat_row_for_partial_vector(A, step, mat_row_for_minus_expo);

mat_4x8_on_mat_8x8(UiORTO,expo_from_minus_step,Ui);//Вычисление матрицы Ui=UiORTO*expo_from_minus_step

//partial_vector(FF);//Вычисление НУЛЕВОГО вектора частного решения системы ОДУ на шаге

partial_vector_real(expo_from_minus_step, mat_row_for_minus_expo, angle, (-step),FF);// - для движения слева на право

mat_4x8_on_vect_8(UiORTO,FF,ui_2);//Вычисление вектора ui_2=UiORTO*FF

minus(ui_ORTO, ui_2, ui_);//Вычисление вектора ui_=ui_ORTO-ui_2

orto_norm_4x8(Ui, ui_, UiORTO, ui_ORTO);//Ортонормирование для текущего шага по ii

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

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

MATRIXS[ii][i][j]=UiORTO[i][j];

}

VECTORS[ii][i]=ui_ORTO[i];

}

}//Цикл по шагам ii (ВЕРХНЕЕ заполнение)

//Цикл по точкам ii интервала интегрирования заполнения НИЖНИХ частей матричных уравнений MATRIXS[ii]*Y[ii]=VECTORS[ii],

//начиная с предпоследней точки - точки с индексом ii=(100-1) используем ii-- (уменьшение индекса точки)

angle=finish_angle;//Угловая координата правого края

for(int ii=(100-1);ii>=0;ii--){

angle-=step;//Движение справа на лево

A_perem_coef(nju, c2, nn, angle, A);//Вычисление матрицы А коэффициентов системы ОДУ при данной угловой координате angle

exponent(A,step,expo_from_plus_step);//Шаг положительный (значение шага больше нуля из-за направления вычисления матричной экспоненты)

mat_row_for_partial_vector(A, (-step), mat_row_for_plus_expo);

mat_4x8_on_mat_8x8(VjORTO,expo_from_plus_step,Vj);//Вычисление матрицы Vj=VjORTO*expo_from_plus_step

//partial_vector(FF);//Вычисление НУЛЕВОГО вектора частного решения системы ОДУ на шаге

partial_vector_real(expo_from_plus_step, mat_row_for_plus_expo, angle, step,FF);// - для движения справа на лево

mat_4x8_on_vect_8(VjORTO,FF,vj_2);//Вычисление вектора vj_2=VjORTO*FF

minus(vj_ORTO, vj_2, vj_);//Вычисление вектора vj_=vj_ORTO-vj_2

orto_norm_4x8(Vj, vj_, VjORTO, vj_ORTO);//Ортонормирование для текущего шага по ii

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

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

MATRIXS[ii][i+4][j]=VjORTO[i][j];

}

VECTORS[ii][i+4]=vj_ORTO[i];

}

}//Цикл по шагам ii (НИЖНЕЕ заполнение)

//Решение систем линейных алгебраических уравнений

for(int ii=0;ii<=100;ii++){

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

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

MATRIX_2[i][j]=MATRIXS[ii][i][j];//Вспомогательное присвоение для соответствия типов в вызывающей функции GAUSS

}

VECTOR_2[i]=VECTORS[ii][i];//Вспомогательное присвоение для соответствия типов в вызывающей функции GAUSS

}

GAUSS(MATRIX_2,VECTOR_2,Y_2);

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

Y[ii][i]=Y_2[i];

}

}

//Вычисление момента во всех точках между краями

angle=start_angle;//начальное значение угловой координаты

for(int ii=0;ii<=100;ii++){

Moment[ii]+=Y[ii][4]*(-nju*nn2/sin(angle)/sin(angle))+Y[ii][5]*(nju/tan(angle))+Y[ii][6]*(1.0);//Момент M1 в точке [ii]

angle+=step;

//U[2][4]=-nju*nn2/sin(start_angle)/sin(start_angle);

//U[2][5]=nju/tan(start_angle);

//U[2][6]=1.0; Момент

}

}//ЦИКЛ ПО ГАРМОНИКАМ ЗДЕСЬ ЗАКАНЧИВАЕТСЯ

for(int ii=0;ii<=100;ii++){

fprintf(fp,"%f\n",Moment[ii]);

}

fclose(fp);

printf( "PRESS any key to continue...\n" );

_getch();

return 0;

}



Приложение 3

Программа на С++ расчета цилиндра - для метода из главы 7.

Вычислительные эксперименты проводились в сравнении с методом переноса краевых условий Алексея Юрьевича Виноградова. В этом методе используется построчное ортонормирование.

Без ортонормирования в методе переноса краевых условий А.Ю.Виноградова успешно решается задача, например, нагружения цилиндрической оболочки, которая консольно заделана по правому краю и нагружена по левому краю силой, равномерно распределенной по дуге окружности, с отношением длинны к радиусу L/R=2 и с отношением радиуса к толщине R/h=100. Для отношения R/h=200 задача без ортонормирования в методе переноса краевых условий уже не решается, так как выдаются ошибки из-за неустойчивости счета. С применением же ортонормирования в методе переноса краевых условий решаются успешно задачи и для параметров, например, R/h=300, R/h=500, R/h=1000.

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

Для тестовых расчетов задач с вышеуказанными параметрами новым предлагаемым методом интервал интегрирования разделялся на 10 участков, а между узлами, как и сказано выше, решение находилось как решение задачи Коши. Для решения задач удерживалось 50 гармоник рядов Фурье, так как результат при 50 гармониках уже не отличался от случая удержания 100 гармоник.

Скорость же расчета тестовых задач новым предлагаемым методом не меньше, чем методом переноса краевых условий, так как оба метода в тестовых задачах при удержании 50 гармоник рядов Фурье выдавали готовое решение мгновенно после запуска программы на выполнение (на ноутбуке ASUS M51V CPU Duo T5800). В тоже время программирование нового предложенного здесь метода существенно проще, так как нет необходимости программировать процедуры ортонормирования.

ПРОГРАММА НА С++ (ЦИЛИНДР):

// sopryazhenie.cpp: главный файл проекта.

//Решение краевой задачи - цилиндрической оболочки.

//Интервал интегрирования разбит на 10 сопрягаемых участков: левый край - точка 0 и правый край - точка 10

//БЕЗ ОРТОНОРМИРОВАНИЯ

#include "stdafx.h"

#include <iostream>

#include <conio.h>

using namespace std;

//Умножение матрицы A на вектор b и получаем rezult.

void mat_on_vect(double A[8][8], double b[8], double rezult[8]){

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

rezult[i]=0.0;

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

rezult[i]+=A[i][k]*b[k];

}

}

}

//Вычисление матричной экспоненты EXP=exp(A*delta_x)

void exponent(double A[8][8], double delta_x, double EXP[8][8]) {

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

int n=100, m;

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

}

}

}

}

}

//Задание вектора внешних воздействий в системе ОДУ - вектора POWER: Y'(x)=A*Y(x)+POWER(x):

void power_vector_for_partial_vector(double x, double POWER[8]){

POWER[0]=0.0;

POWER[1]=0.0;

POWER[2]=0.0;

POWER[3]=0.0;

POWER[4]=0.0;

POWER[5]=0.0;

POWER[6]=0.0;

POWER[7]=0.0;

}

//Вычисление 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;

}

}

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

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

double A[8*11][8*11];

double b[8*11];

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

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

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

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

}

}

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

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

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

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

for(int i=jj;i<(8*11);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*11);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*11);i++){//Приведение к верхнетреугольной матрице

for(int j=(jj+1);j<(8*11);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*11)-1]=b[(8*11)-1]/A[(8*11)-1][(8*11)-1];//Первоначальное определение последнего элемента искомого решения х (87-го)

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

t=0;

for(int j=1;j<((8*11)-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 step=0.05; //step=(L/R)/100 - величина шага расчета оболочки - шага интервала интегрирования (должна быть больше нуля, т.е. положительная)

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

h_div_R=1.0/100;

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

printf( "The file 'C:/test.txt' was not opened\n" );

else

printf( "The file 'C:/test.txt' was opened\n" );

for(nn=1;nn<=nn_last;nn++){ //ЦИКЛ ПО ГАРМОНИКАМ, НАЧИНАЯ С 1-ОЙ ГАРМОНИКИ (БЕЗ НУЛЕВОЙ ГАРМОНИКИ)

double x=0.0;//Координата от левого края - нужна для случая неоднородной системы ОДУ для вычисления частного вектора FF

double expo_from_minus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (0-x1)

double expo_from_plus_step[8][8]={0};//Матрица для расположения в ней экспоненты на шаге типа (x1-0)

double mat_row_for_minus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (0-x1)

double mat_row_for_plus_expo[8][8]={0};//вспомогательная матрица для расчета частного вектора при движении на шаге типа (x1-0)

double U[4][8]={0};//Матрица краевых условий левого края размерности 4х8

double u_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий левого края

double V[4][8]={0};//Матрица краевых условий правого края размерности 4х8

double v_[4]={0};//Вектор размерности 4 внешнего воздействия для краевых условий правого края

double Y[100+1][8]={0};//Массив векторов-решений соответствующих СЛАУ (в каждой точке интервала между краями): MATRIXS*Y=VECTORS

double A[8][8]={0};//Матрица коэффициентов системы ОДУ

double FF[8]={0};//Вектор частного решения неоднородной ОДУ на участке интервала интегрирования

double Y_many[8*11]={0};// составной вектор из векторов Y(xi) в 11-ти точках с точки 0 (левый край Y(0)) до точки 10 (правый край Y(x10))

double MATRIX_many[8*11][8*11]={0};//матрица СЛАУ

double B_many[8*11]={0};// вектор правых частей СЛАУ: MATRIX_many*Y_many=B_many

double Y_vspom[8]={0};//вспомогательный вектор

double Y_rezult[8]={0};//вспомогательный вектор

double nn2,nn3,nn4,nn5,nn6,nn7,nn8;//Возведенный в соответствующие степени номер гармоники nn

nn2=nn*nn; nn3=nn2*nn; nn4=nn2*nn2; nn5=nn4*nn; nn6=nn4*nn2; nn7=nn6*nn; nn8=nn4*nn4;

//Заполнение ненулевых элементов матрицы А коэффициентов системы ОДУ

A[0][1]=1.0;

A[1][0]=(1-nju)/2*nn2; A[1][3]=-(1+nju)/2*nn; A[1][5]=-nju;

A[2][3]=1.0;

A[3][1]=(1+nju)/(1-nju)*nn; A[3][2]=2*nn2/(1-nju); A[3][4]=2*nn/(1-nju);

A[4][5]=1.0;

A[5][6]=1.0;

A[6][7]=1.0;

A[7][1]=-nju/c2; A[7][2]=-nn/c2; A[7][4]=-(nn4+1/c2); A[7][6]=2*nn2;

//Здесь надо первоначально заполнить ненулевыми значениями матрицы и вектора краевых условий U*Y[0]=u_ (слева) и V*Y[100]=v_ (справа) :

U[0][1]=1.0; U[0][2]=nn*nju; U[0][4]=nju; u_[0]=0.0;//Сила T1 на левом крае равна нулю

U[1][0]=-(1-nju)/2*nn; U[1][3]=(1-nju)/2; U[1][5]=(1-nju)*nn*c2; u_[1]=0.0;//Сила S* на левом краю равна нулю

U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Момент M1 на левом краю равен нулю

U[3][5]=(2-nju)*nn2; U[3][7]=-1.0;

u_[3]=-sin(nn*gamma)/(nn*gamma);//Сила Q1* на левом крае распределена на угол -gamma +gamma

V[0][0]=1.0; v_[0]=0.0;//Перемещение u на правом крае равно нулю

V[1][2]=1.0; v_[1]=0.0;//Перемещение v на правом крае равно нулю

V[2][4]=1.0; v_[2]=0.0;//Перемещение w на правом крае равно нулю

V[3][5]=1.0; v_[3]=0.0;//Угол поворота на правом крае равен нулю

//Здесь заканчивается первоначальное заполнение U*Y[0]=u_ и V*Y[100]=v_

exponent(A,(-step*10),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты)

//x=0.0;//начальное значение координаты - для расчета частного вектора

//mat_row_for_partial_vector(A, step, mat_row_for_minus_expo);

//Заполнение матрицы коэффициентов СЛАУ MATRIX_many

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

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

MATRIX_many[i][j]=U[i][j];

MATRIX_many[8*11-4+i][8*11-8+j]=V[i][j];

}

B_many[i]=u_[i];

B_many[8*11-4+i]=v_[i];

}

for(int kk=0;kk<(11-1);kk++){//(11-1) единичных матриц и матриц EXPO надо записать в MATRIX_many

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

MATRIX_many[i+4+kk*8][i+kk*8]=1.0;//заполнение единичными матрицами

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

MATRIX_many[i+4+kk*8][j+8+kk*8]=-expo_from_minus_step[i][j];//заполнение матричными экспонентами

}

}

}

//Решение систем линейных алгебраических уравнений

GAUSS(MATRIX_many,B_many,Y_many);

//Вычисление векторов состояния в 101 точке - левая точка 0 и правая точка 100

exponent(A,step,expo_from_plus_step);

for(int i=0;i<11;i++){//заполнение промежуточных точек во всех 10-ти интервалах (всего получим точки от 0 до 100) между 11 узлами

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

Y[0+i*10][j]=Y_many[j+i*8];//в 11-ти узлах векторы беруться из решения СЛАУ - из Y_many

}

}

for(int i=0;i<10;i++){//заполнение промежуточных точек в 10-ти интервалах

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

Y_vspom[j]=Y[0+i*10][j];//начальный вектор для i-го участка, нулевая точка, точка старта i-го участка

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+1][j]=Y_rezult[j];//заполнение 1-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+2][j]=Y_rezult[j];//заполнение 2-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+3][j]=Y_rezult[j];//заполнение 3-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+4][j]=Y_rezult[j];//заполнение 4-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+5][j]=Y_rezult[j];//заполнение 5-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+6][j]=Y_rezult[j];//заполнение 6-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+7][j]=Y_rezult[j];//заполнение 7-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+8][j]=Y_rezult[j];//заполнение 8-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

mat_on_vect(expo_from_plus_step, Y_vspom, Y_rezult);

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

Y[0+i*10+9][j]=Y_rezult[j];//заполнение 9-ой точки интервала

Y_vspom[j]=Y_rezult[j];//для следующего шага

}

}

//Вычисление момента во всех точках между краями

for(int ii=0;ii<=100;ii++){

Moment[ii]+=Y[ii][4]*(-nju*nn2)+Y[ii][6]*1.0;//Момент M1 в точке [ii]

//U[2][4]=-nju*nn2; U[2][6]=1.0; u_[2]=0;//Момент M1

}

}//ЦИКЛ ПО ГАРМОНИКАМ ЗДЕСЬ ЗАКАНЧИВАЕТСЯ

for(int ii=0;ii<=100;ii++){

fprintf(fp,"%f\n",Moment[ii]);

}

fclose(fp);

printf( "PRESS any key to continue...\n" );

_getch();

return 0;

}

Приложение 4

Метод главы 7 и программа для этого метода на С++ на английском языке.

The simplest method of solution of the stiff boundary value problems -

the method of solution of the stiff boundary value problems without orthonormalization

offered by Y. I. Vinogradov and A. Y. Vinogradov.

Alexei Yurievich Vinogradov (Candidate of science /Physics and Mathematics/)

(14 November 2011, Moscow, Russia)

The authors: Yuri Ivanovich Vinogradov (1938 year of birth), doctor of science (Physics and Mathematics) and Alexei Yurievich Vinogradov (1970 year of birth), candidate of science (Physics and Mathematics).

The idea of overcoming the problems of computing by dividing an integration interval into matching sectors belongs to Y. I. Vinogradov (doctor of science /Physics and Mathematics/). His doctorate thesis was based on that material (including that idea).

Proposed realization of this idea (division and matching) through matrix theory formulas i.e. by means of the matrix exponents (or by other words by means of a matrizant - the Cauchy matrix, Cauchy-Krylov functions) belongs to A. Y. Vinogradov (candidate of science /Physics and Mathematics/).

AT PRESENT MOMENT THIS METHOD IS THE BEST ONE. The method was finalized in the mid October 2011 and a C++ software program (that is given in this paper as well) was finalized and checked by comparative computations on 14 November 2011.

1. The simplest method of solution of the stiff boundary value problems.

Let's divide a boundary value problem integration interval e.g. for three segments. We will have points (nodes) including the boundaries:

.

We have the boundary conditions in the following form:

We can write sector matching matrix equations as follows:

It is possible to rewrite these formulas in the more convenient form:



where E - unit matrix.

Then in the combined matrix form we obtain a system of the linear algebraic equations in the following form:

This system is solved by the Gauss method by discrimination of the basic element.

Solution in the points located between the nodes is obtained by solving the Cauchy problem with initial conditions in the ith node:

It appears to be that there is no need to apply orthonormalization for the boundary value problems for stiff ordinary differential equations.

2. Detailed description of the designations used above.

Let's use as an example a system of differential equations of the cylindrical shell of the rocket - a system of the ordinary differential equations of 8th order (after partial derivatives partitioning by the Fourier method).

A system of the linear ordinary differential equations has the following form:

where - an unknown vector-function of the 8x1 dimensionality problem, - a derivative of the unknown vector-function of 8x1 dimensionality, - a square matrix of the coefficients of the differential equation of 8x8 dimensionality, - a vector-function of the external influence on the system of 8x1 dimensionality.

The boundary conditions have the following form:

where - is the value of the unknown vector-function at the left boundary x=0 of the 8x1 dimensionality, - a square horizontal matrix of the coefficients of the boundary conditions of the left boundary of the 4x8 dimensionality, and - an external influence on the left boundary of the 4x1 dimensionality.

is the value of the unknown vector-function at the right boundary x=1 of the 8x1 dimensionality, - a square horizontal matrix of the coefficients of the boundary conditions of the right boundary of the 4x8 dimensionality, and - an external influence on the left boundary of the 4x1 dimensionality.

In case when a system of differential equations has a matrix with constant coefficients = const the solution of the Cauchy problem has the following form:

where ,

where is a unit matrix.

Matrix exponent can be named as the Cauchy matrix or matrizant and be denoted in the following form:

Then Cauchy problem solution can be written in the following form:



where is a vector of the partial solution of the heterogeneous system of the differential equations.

A matrix exponent (the Cauchy matrix) multiplication feature is known from the matrix theory [1]:

In case when a differential equation system has a matrix with variable coefficients it is offered to find a solution of the Cauchy problem by means of the Cauchy matrix multiplication feature. That is the integration integral is split into small segments and the Cauchy matrixes are computed approximately in these small segments by means of the formula for a constant matrix in the exponent. Then the Cauchy matrixes computed in these small sections are multiplied:

where the Cauchy matrixes are computed approximately by the following formula:

where

Instead of the formula for computation of the vector of the partial solution of the heterogeneous system of the differential equations in the form [1]:



it is proposed to use the following formula for each individual segment of the integration interval:

Correctness of the above mentioned formula is justified by the following:

as was to be justified.

Computation of the vector of the partial solution of the heterogeneous system of differential equations is carried out by representing the Cauchy matrix under the integral sign in the form of the series and integrating this series element by element:

This formula is correct for a system of differential equations with the constant coefficient matrix =const.

vector can be considered in the segment approximately in the form of the constant value that allows its taking off from the integral sign that results in very simple series for computation in the considered segment.

In case of the differential equations with variable coefficients the averaged matrix of the coefficients of the differential equation system can be used in the above mentioned formula for each segment.

Let's consider an option when the integration interval steps are selected as sufficiently small that results in considering vector in the segment approximately in the form of the constant value resulting in taking off this vector from the integral signs:

It is known that at T=(at+b) we have:

In our case we have:

Then we obtain

.

Then we obtain the series for computation of the vector of the partial solution of the heterogeneous system of differential equations in the small segment :

If the segment is not a small one it is possible to divide it by sub-segments and then it would be possible to propose the following recurrent (iteration) formulas for particular vector computing:

We have .

We have the formula for an individual sub-segment as well:

We can write:



,

.

Let's substitute in and obtain the following:

Let's compare the obtained expression with the formula:

and we will obtain evidently that:

so, we obtain the formula for a particular vector:

Hence, the vectors of the sub-segments are not added to each other simply but they are added with involvement of the Cauchy matrix of the sub-segment.

Similarly we can write and having substituted the formula for in this expression we will obtain:



having compared this obtained expression with the formula:

we have obtained evidently that:

and at the same time we have obtained a formula for a particular vector:

Hence, a particular vector - a vector of the partial solution of the heterogeneous system of differential equations is computed by this way exactly, for example, the particular vector is computed similarly in the considered segment through the computed particular vectors , , and of the corresponding subs-segments , and .

Computation accuracy control can be carried out by the following way.

For a homogeneous system of differential equations we have:

.

We may write the following:

and

Having substituted one formula in another one we obtain:

,

i.e. we obtain:

,

however, the latter is possible only in case when

- unit matrix,

i.e. the matrixes and are inverse ones.

In other words it has been proved that

.

Hence, computing accuracy can be controlled by means of determination of the accuracy of matrix exponent computation (in other words the Cauchy matrixes or matrizants), that can be checked by having ascertained that the reciprocal inversion of the corresponding matrix exponents is observed in the integration segments:

.

3. Computational experiments.

Computational experiments have been carried out in comparison with the method of the boundary conditions transfer of Alexey Vinogradov. Line-by-line orthonormalization is used in this method.

Without use of orthonormalization it is possible by means of the boundary conditions transfer method to solve a problem of cylindrical shell loading that is fixed in cantilever fashion at the right boundary and loaded at the left boundary by the force distributed uniformly by the circular arch with the ratio of the length to the radius L/R=2 and with ratio of the radius to the thickness R/h=100. In case of ratio R/h=200 the problem by means of the method of the boundary conditions transfer without orthonormalozation cannot be solved by this time because there are mistakes resulted in due to counting instability. However, in case of use of orthonormalization in the method of the boundary conditions transfer the problems related to the parameters R/h=300, R/h=500, R/h=1000 can be solved.

A new method proposed in this paper allows solving of all above mentioned test problems without use of orthonormalization operation that results in significant simplification of its programming.

In case of the test computations of the problems characterized by the above mentioned parameters by means of this new proposed method the integration interval is divided into 10 segments while between the nodes as aforesaid the solution was found as a solution of the Cauchy problem. 50 harmonics of the Fourier series were used for solving the problem since the result in case of usage of 50 harmonics didn't differ from the case when 100 harmonics were used.

Test problem computing speed by means of the proposed method is not less compared to the boundary conditions transfer method because both methods when used for test problems while using 50 harmonics of the Fourier series produced a final solution instantaneously after launching a program (notebook ASUS M51V CPU Duo T5800). At the same time programming of this newly proposed method is significantly simpler because there is no need in orthonormalization procedure programming.

LIST OF REFERENCES

1. Gantmaher F.R. Matrix theory. - M.: Nauka, 1988. - 548 p.

C++ program

// sopryazhenie.cpp: main file of the project.

//Solution of the boundary value problem - a cylindrical shell problem.

//Integration interval is divided into 10 matching segments: left boundary - point 0 and right boundary - point 10.

//WITHOUT ORTHONORMALIZATION

#include "stdafx.h"

#include <iostream>

#include <conio.h>

using namespace std;

//Multiplication of A matrix by b vector and obtain rezult.

void mat_on_vect(double A[8][8], double b[8], double rezult[8]){

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

rezult[i]=0.0;

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

rezult[i]+=A[i][k]*b[k];

}

}

}

//Computation of the matrix exponent EXP=exp(A*delta_x)

void exponent(double A[8][8], double delta_x, double EXP[8][8]) {

//n - number of the terms of the series in the exponent, m - a counter of the number of the terms of the series (m<=n)

int n=100, m;

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

int i,j,k;

//E - unit matrix - the first term of the series of the exponent

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;

//initial filling-in of the auxiliary array TMP1 - the previous term of the series for follow-up multiplication

//and initial filling-in of the exponent by the first term of the series

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

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

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

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

}

}

//series of EXP exponent computation starting from the 2nd term of the series (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;//taken out beyond the cycle of multiplication of the row by the column

TMP2[i][j]/=(m-1);// taken out beyond the cycle of multiplication of the row by the column

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

}

}

//filling-in of the auxiliary array TMP1 for computing the next term of the series TMP2 in the next step of the cycle by m

if (m<n) {

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

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

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

}

}

}

}

}

//computation of the matrix MAT_ROW in the form of the matrix series for follow-up use

//when computing a vector - partial vector - a vector of the partial solution of the heterogeneous system of the ordinary differential equations at the step delta x

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

//n - number of the terms of the series in MAT_ROW, m - a counter of the number of the terms of the series (m<=n)

int n=100, m;

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

int i,j,k;

//E - unit matrix - the first term of the series 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;

//initial filling-in of the auxiliary array TMP1 - the previous term of the series for following multiplication

//and initial filling-in of MAT_ROW by the first term of the series

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

}

}

//a series of computation of MAT_ROW starting from the second term of the series (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];

}

}

//filling-in of the auxiliary array TMP1 for computing the next term of the series - TMP2 in the next step of the cycle by m

if (m<n) {

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

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

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

}

}

}

}

}

//specifying the external influence vector in the system of ordinary differential equations - POWER vector: Y'(x)=A*Y(x)+POWER(x):

void power_vector_for_partial_vector(double x, double POWER[8]){

POWER[0]=0.0;

POWER[1]=0.0;

POWER[2]=0.0;

POWER[3]=0.0;

POWER[4]=0.0;

POWER[5]=0.0;

POWER[6]=0.0;

POWER[7]=0.0;

}

//computation of the vector - ZERO (particular case) vector of the partial solution

//heterogeneous system of differential equations in the segment of interest:

void partial_vector(double vector[8]){

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

vector[i]=0.0;

}

}

//computation of the vector - the vector of the partial solution of the heterogeneous system of differential equations in the segment of interest 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};//External influence vector on the shell

double REZ[8]={0};

double REZ_2[8]={0};

power_vector_for_partial_vector(x_, POWER_);//Computing POWER_ at coordinate x_

mat_on_vect(mat_row, POWER_, REZ);//Multiplication of the matrix mat_row by POWER vector and obtain REZ vector

mat_on_vect(expo_, REZ, REZ_2);// Multiplication of matrix expo_ by vector REZ and obtain vector REZ_2

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

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

}

}

//Solution of SLAE of 88 dimensionality by the Gauss method with discrimination of the basic element

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

double A[8*11][8*11];

double b[8*11];

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

b[i]=bb[i];//we will operate with the vector of the и right parts to provide that initial vector bb would not change when exiting the subprogram

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

A[i][j]=AA[i][j];//we will operate with A matrix to provide that initial AA matrix would not change when exiting the subprogram

}

}

int e;//number of the row where main (maximal) coefficient in the column jj is found

double s, t, main;//Ancillary quantity

for(int jj=0;jj<((8*11)-1);jj++){//Cycle by columns jj of transformation of A matrix into upper-triangle one

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

for(int i=jj;i<(8*11);i++){//there is a number of у row where main (maximal) element is placed in the column jj and row interchanging is made

if ((A[i][jj]*A[i][jj])>s) {//Instead of multiplication (potential minus sign is deleted) it could be possible to use a function by abs() module

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

}

}

if (e<0) {

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

}

if (e>jj) {//If the main element isn't placed in the row with jj number but is placed in the row with у number

main=A[e][jj];

for(int j=0;j<(8*11);j++){//interchanging of two rows with e and jj numbers

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*11);i++){//reduction to the upper-triangle matrix

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

A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//re-calculation of the coefficients of the row i>(jj+1)

}

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

A[i][jj]=0.0;//nullified elements of the row under diagonal element of A matrix

}

}//Cycle by jj columns of transformation of A matrix into upper-triangle one

x[(8*11)-1]=b[(8*11)-1]/A[(8*11)-1][(8*11)-1];//initial determination of the last element of the desired solution x (87th)

for(int i=((8*11)-2);i>=0;i--){//Computation of the elements of the solution x[i] from 86th to 0th

t=0;

for(int j=1;j<((8*11)-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;//Number of the harmonic starting from the 1st (without zero one)

int nn_last=50;//Number of the last harmonic

double Moment[100+1]={0};//An array of the physical parameter (momentum) that is computed in each point between the boundaries

double step=0.05; //step=(L/R)/100 - step size of shell computation - a step of integration interval (it should be over zero, i.e. be positive)

double h_div_R;//Value of h/R

h_div_R=1.0/100;

double c2;

c2=h_div_R*h_div_R/12;//Value of h*h/R/R/12

double nju;

nju=0.3;

double gamma;

gamma=3.14159265359/4;//The force distribution angle by the left boundary

//printing to files:

FILE *fp;

// Open for write

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

printf( "The file 'C:/test.txt' was not opened\n" );

else

printf( "The file 'C:/test.txt' was opened\n" );

...