Учебные материалы


06 сентября 2011 (расчет сферической оболочки) Алексей Юрьевич Виноградов Кандидат физико-математических наук (1996 года защиты) Дата рождения: 12 апреля 1970 (а то в интернете много моих полных тезок)



Карта сайта diamantesdetequila.com 1 2 3 06 сентября 2011 (расчет СФЕРИЧЕСКОЙ ОБОЛОЧКИ) Алексей Юрьевич Виноградов Кандидат физико-математических наук (1996 года защиты) Дата рождения: 12 апреля 1970 (а то в интернете много моих полных тезок) Мои сайты по методам решения краевых задач в интернете: www.AlexeiVinogradov.narod.ru www.VinogradovAlexei.narod.ru www.Vinogradov-Alexei.narod.ru www.Vinogradov-Math.narod.ru СОДЕРЖАНИЕ:

  • Программа (код) на С++, написанная в среде MS Visual Studio 2010 (Visual C++), – программа решения «жесткой» краевой задачи для системы обыкновенных дифференциальных уравнений с ПЕРЕМЕННЫМИ коэффициентами – СФЕРА. Страницы 2 – 13.
  • Теория метода Алексея Юрьевича Виноградова «переноса краевых условий» для краевых задач, включая «жесткие» краевые задачи, для реализации которого написана приводимая программа. Страницы 13 - 28. //sfera_from_A_Yu_Vinogradov.cpp: главный файл проекта. //Решение краевой задачи с переменными коэффициентами - сфера. //Интервал интегрирования разбит на 100 участков: левый край - точка 0 и правый край - точка 100 #include "stdafx.h" #include #include #include //for tan() using namespace std; //Вычисление для гармоники с номером nn для значения переменной (угла) angle_fi - матрицы A_perem 8x8 коэффициентов системы ОДУ void A_perem_coef(double nju, double c2, int nn, double angle_fi, double A_perem[8][8]){ 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; for(int i=0;i<8;i++){ for(int j=0;j<8;j++){ A_perem[i][j]=0.0;//Первоначальное обнуление матрицы } } //Заполнение ненулевых элементов матрицы А коэффициентов системы ОДУ A_perem[0][1]=1.0; A_perem[1][0]=(1-nju)*nn2/2/sin(angle_fi)/sin(angle_fi)+nju+1.0/tan(angle_fi)/tan(angle_fi); A_perem[1][1]=-1.0/tan(angle_fi); A_perem[1][2]=(3-nju)/2/sin(angle_fi)/tan(angle_fi); A_perem[1][3]=-(1+nju)*nn/2/sin(angle_fi); A_perem[1][5]=-(1+nju); A_perem[2][3]=1.0; A_perem[3][0]=(3-nju)*nn/(1-nju)/sin(angle_fi)/tan(angle_fi); A_perem[3][1]=(1+nju)*nn/(1-nju)/sin(angle_fi); A_perem[3][2]=2*nn2/(1-nju)/sin(angle_fi)/sin(angle_fi)-1.0+1.0/tan(angle_fi)/tan(angle_fi); A_perem[3][3]=-1.0/tan(angle_fi); A_perem[3][4]=(1+nju)*2*nn/(1-nju)/sin(angle_fi); A_perem[4][5]=1.0; A_perem[5][6]=1.0; A_perem[6][7]=1.0; A_perem[7][0]=-(1+nju)/tan(angle_fi)/c2; A_perem[7][1]=-(1+nju)/c2; A_perem[7][2]=-(1+nju)*nn/c2/sin(angle_fi); A_perem[7][4]=nn2/sin(angle_fi)/sin(angle_fi)*(2+(2-nn2)/sin(angle_fi)/sin(angle_fi)+2.0/tan(angle_fi)/tan(angle_fi))-2*(1+nju)/c2; A_perem[7][5]=(-2.0-(2*nn2+1)/sin(angle_fi)/sin(angle_fi))/tan(angle_fi); A_perem[7][6]=-1.0+(2*nn2+1)/sin(angle_fi)/sin(angle_fi); A_perem[7][7]=-2.0/tan(angle_fi); } //Задание вектора внешних воздействий в системе ОДУ - вектора 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; } //Скалярное произведение векторов - i-й строки матрицы А и j-й строки матрицы С. double mult(double A[8][8], int i, double C[8][8], int j){ double result=0.0; for(int k=0;k<8;k++){ result+=A[i][k]*C[j][k]; } return result; } //Вычисление нормы вектора, где вектор это i-я строка матрицы А. double norma(double A[8][8], int i){ double norma_=0.0; for(int k=0;k<8;k++){ norma_+=A[i][k]*A[i][k]; } norma_=sqrt(norma_); return norma_; } //Выполнение ортонормирования. Исходная система A*x=b размерности 8х8 приводиться к системе C*x=d, где строки матрицы С ортонормированы. void orto_norm_8x8(double A[8][8], double b[8], double C[8][8], double d[8]){ 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; //Получаем 5-ю строку уравнения C*x=d: mult0=mult(A,4,C,0); mult1=mult(A,4,C,1); mult2=mult(A,4,C,2); mult3=mult(A,4,C,3); for(int k=0;k<8;k++){ C[4][k]=A[4][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]- mult3*C[3][k]; } NORM=norma(C,4); for(int k=0;k<8;k++){ C[4][k]/=NORM; } d[4]=(b[4]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3])/NORM; //Получаем 6-ю строку уравнения C*x=d: mult0=mult(A,5,C,0); mult1=mult(A,5,C,1); mult2=mult(A,5,C,2); mult3=mult(A,5,C,3); mult4=mult(A,5,C,4); for(int k=0;k<8;k++){ C[5][k]=A[5][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]- mult3*C[3][k]-mult4*C[4][k]; } NORM=norma(C,5); for(int k=0;k<8;k++){ C[5][k]/=NORM; } d[5]=(b[5]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4])/NORM; //Получаем 7-ю строку уравнения C*x=d: mult0=mult(A,6,C,0); mult1=mult(A,6,C,1); mult2=mult(A,6,C,2); mult3=mult(A,6,C,3); mult4=mult(A,6,C,4); mult5=mult(A,6,C,5); for(int k=0;k<8;k++){ C[6][k]=A[6][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]- mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k]; } NORM=norma(C,6); for(int k=0;k<8;k++){ C[6][k]/=NORM; } d[6]=(b[6]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]- mult5*d[5])/NORM; //Получаем 8-ю строку уравнения C*x=d: mult0=mult(A,7,C,0); mult1=mult(A,7,C,1); mult2=mult(A,7,C,2); mult3=mult(A,7,C,3); mult4=mult(A,7,C,4); mult5=mult(A,7,C,5); mult6=mult(A,7,C,6); for(int k=0;k<8;k++){ C[7][k]=A[7][k]-mult0*C[0][k]-mult1*C[1][k]-mult2*C[2][k]- mult3*C[3][k]-mult4*C[4][k]-mult5*C[5][k]-mult6*C[6][k]; } NORM=norma(C,7); for(int k=0;k<8;k++){ C[7][k]/=NORM; } d[7]=(b[7]-mult0*d[0]-mult1*d[1]-mult2*d[2]-mult3*d[3]-mult4*d[4]- mult5*d[5]-mult6*d[6])/NORM; } //Выполнение ортонормирования системы A*x=b с прямоугольной матрицей A коэффициентов размерности 4х8. 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 (ms) {//Вместо перемножения (удаляется возможный знак минуса) можно было бы использовать функцию по модулю abs() e=i; s=A[i][jj]*A[i][jj]; } } if (e<0) { cout<<"Mistake "<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; } 1 2 3


  • edu 2018 год. Все права принадлежат их авторам! Главная