Программа (код) на С++ решения жесткой краевой задачи методом А.Ю.Виноградова
| Примечание | от автора: Программа окончательно написана и проверена счетом 30 июня 2011; от редактора: приведен только текст, формулы и графики статьи находятся в архивном файле |
| Загрузить архив: | |
| Файл: ref-31702.zip (271kb [zip], Скачиваний: 1) скачать |
06 сентября 2011 (расчет ЦИЛИНДРИЧЕСКОЙ ОБОЛОЧКИ)
Алексей Юрьевич Виноградов
Кандидат физико-математических наук (1996 года защиты)
Дата рождения: 12 апреля 1970 (а то в интернете много моих полных тезок)
Мои сайты по методам решения краевых задач в интернете:
СОДЕРЖАНИЕ:
//from_A_Yu_Vinogradov.cpp: главный файл проекта.
//Решение краевой задачи - цилиндрической оболочки.
//Интервал интегрирования разбит на 100 участков: левый край - точка 0 и правый край - точка 100
#include "stdafx.h"
#include
#include
using namespace std;
//Скалярное произведение векторов - 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 (m 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 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; } } //Решение СЛАУ размерности 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 "< } 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 step=0.02; //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 openedn" ); else printf( "The file 'C:/test.txt' was openedn" ); 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 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; //Заполнение ненулевых элементов матрицы А коэффициентов системы ОДУ 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_ 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 exponent(A,(-step),expo_from_minus_step);//Шаг отрицательный (значение шага меньше нуля из-за направления вычисления матричной экспоненты) x=0.0;//начальное значение координаты - для расчета частного вектора mat_row_for_partial_vector(A, step, mat_row_for_minus_expo); for(int ii=1;ii<=100;ii++){ x+=step;//Координата для расчета частного вектора на шаге 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, x, (-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-- (уменьшение индекса точки) exponent(A,step,expo_from_plus_step);//Шаг положительный (значение шага больше нуля из-за направления вычисления матричной экспоненты) x=step*100;//Координатаправогокрая mat_row_for_partial_vector(A, (-step), mat_row_for_plus_expo); for(int ii=(100-1);ii>=0;ii--){ x-=step;//Движение справа на лево - для расчета частного вектора 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, x,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]; } } //Вычисление момента во всех точках между краями 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,"%fn",Moment[ii]); } fclose(fp); printf( "PRESS any key to continue...n" ); _getch(); return 0; } Метод Алексея Юрьевича Виноградова «переноса краевых условий» для решения краевых задач, включая «жесткие» краевые задачи. 1. Введение. Метод проверен и его эффективность выше, чем эффективность (скорость счета) метода ортогональной прогонки С.К.Годунова - до 2 порядков (в 100 раз). Для разных типов задач преимущество по скорости разное. А для некоторых типов задач метод А.Ю.Виноградова не требует применять ортонормирование вовсе. Рассмотрим метод на примере системы дифференциальных уравнений цилиндрической оболочки ракеты – системы обыкновенных дифференциальных уравнений 8-го порядка (после разделения частных производных). Система линейных обыкновенных дифференциальных уравнений имеет вид: где Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами Краевые условия имеют вид: где В случае, когда система дифференциальных уравнений имеет матрицу с постоянными коэффициентами где Матричная экспонента ещё может называться матрицей Коши или матрициантом и может обозначаться в виде: Тогда решение задачи Коши может быть записано в виде: где 2. Случай переменных коэффициентов. Из теории матриц [1] известно свойство перемножаемости матричных экспонент (матриц Коши):
В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами где матрицы Коши приближенно вычисляются по формуле: 3. Вычисление вектора частного решения неоднородной системы дифференциальных уравнений. Вместо формулы для вычисления вектора частного решения неоднородной системы дифференциальных уравнений в виде [1]: предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования: Правильность приведенной формулы подтверждается следующим: что и требовалось подтвердить. Вычисление вектора частного решения неоднородной системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно: Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов Вектор Для случая дифференциальных уравнений с переменными коэффициентами в приведенной выше формуле для каждого участка может использоваться осредненная матрица Рассмотрим вариант, когда шаги интервала интегрирования выбираются достаточно малыми, что позволяет рассматривать вектор Известно, что при T=(at+b) имеем В нашем случае имеем Тогда получаем Тогда получаем ряд для вычисления вектора частного решения неоднородной системы дифференциальных уравнений на малом участке Приведем формулы вычисления вектора частного решения, например, Имеем Также имеем формулу для отдельного подучастка: Можем записать: Подставим Сравним полученное выражение с формулой: и получим, очевидно, что: и для частного вектора получаем формулу: То есть вектора подучастков Аналогично запишем Сравнив полученное выражение с формулой: очевидно, получаем, что: и вместе с этим получаем формулу для частного вектора: То есть именно так и вычисляется частный вектор – вектор частного решения неоднородной системы дифференциальных уравнений, то есть так вычисляется, например, частный вектор 4. Метод «переноса краевых условий» в произвольную точку интервала интегрирования. Полное решение системы дифференциальных уравнений имеет вид Или можно записать: Подставляем это выражение для Или получаем краевые условия, перенесенные в точку где Далее запишем аналогично И подставим это выражение для Или получаем краевые условия, перенесенные в точку где И так в точку Покажем шаги переноса краевых условий с правого края. Можем записать: Подставляем это выражение для Или получаем краевые условия правого края, перенесенные в точку где Далее запишем аналогично И подставим это выражение для Или получаем краевые условия, перенесенные в точку где И так во внутреннюю точку Из этих двух матричных уравнений с прямоугольными горизонтальными матрицами коэффициентов очевидно получаем одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов: 5. Случай «жестких» дифференциальных уравнений. В случае «жестких» дифференциальных уравнений предлагается применять построчное ортонормирование матричных краевых условий в процессе их переноса в рассматриваемую точку. Для этого формулы ортонормирования систем линейных алгебраических уравнений можно взять в [2]. То есть, получив применяем к этой группе линейных алгебраических уравнений построчное ортонормирование и получаем эквивалентное матричное краевое условие: И теперь уже в это проортонормированное построчно уравнение подставляем И получаем Или получаем краевые условия, перенесенные в точку где Теперь уже к этой группе линейных алгебраических уравнений применяем построчное ортонормирование и получаем эквивалентное матричное краевое условие: И так далее. И аналогично поступаем с промежуточными матричными краевыми условиями, переносимыми с правого края в рассматриваемую точку. В итоге получаем систему линейных алгебраических уравнений с квадратной матрицей коэффициентов, состоящую из двух независимо друг от друга поэтапно проортонормированных матричных краевых условий. Эта система решается методом Гаусса с выделением главного элемента для получения решения 6. Контроль точности вычислений. Для однородной системы дифференциальных уравнений имеем: Можем записать: Подставляем одну формулу в другую и получаем: то есть получаем: но последнее возможно только когда то есть матрицы То есть доказано, что то есть Таким образом, получаем, что точность вычислений можно контролировать при помощи определения точности вычисления матричных экспонент (или иначе говоря - матриц Коши или матрициантов), что можно проверять, удостоверяясь, что на участках интегрирования выполняется условие взаимной обратности соответствующих матричных экспонент: 7. Применяемые формулы ортонормирования. Взято из [2]. Пусть дана система линейных алгебраических уравнений порядка n: Здесь над векторами поставим черточки вместо их обозначения жирным шрифтом. Будем рассматривать строки матрицы Ортонормируем эту систему векторов. Первое уравнение системы При этом получим: где Второе уравнение системы заменяется на: где Аналогично поступаем дальше. Уравнение с номером i примет вид: где Здесь ( Процесс будет осуществим, если система линейных алгебраических уравнений линейно независима. В результате мы придем к новой системе 8. Вычислительные эксперименты. В качестве проверочных задач использовалась схема консольно закрепленных цилиндрической и сферической оболочек с параметрами R/h=50, 100, 200. Длина цилиндрической оболочки рассматривалась L/R=2, а угловые координаты сферической оболочки рассматривались от p/4 до 3p/4. На свободном крае рассматривалось нормальное к поверхности оболочек погонное усилие, равномерно распределенное в интервале [-p/4, p/4]. В качестве среды программирования использовалась система MicrosoftVisualStudio 2010 (VisualC++). Первоначально метод был предложен и обсчитывался в кандидатской диссертации А.Ю.Виноградова в 1993-1995 годах. Тогда оказалось, что без использования ортонормирования в рамках метода успешно решаются задачи осесимметрично нагруженных оболочек вращения. Расчеты тогда выполнялись на компьютере поколения 286. Задачи же неосесимметричного нагружения оболочек вращения можно было решать на компьютерах поколения 286 только с применением процедур построчного ортонормирования - как это и предлагалось в рамках метода. Без процедур ортонормирования в неосесимметричных случаях выдавались только ошибочные графики, представлявшие собой хаотично скачущие большие отрицательные и большие положительные значения, например, изгибающего обезразмеренного момента М1. Современные компьютеры имеют значительно более совершенное внутреннее устройство и более точные внутренние операции с числами, чем это было в 1993-1995 годах. Поэтому было интересно рассмотреть возможность расчета неосесимметрично нагруженных оболочек, например, цилиндров, на современном аппаратном и программном обеспечении в рамках предложенного метода «переноса краевых условий» совсем без использования процедур построчного ортонормирования. Оказалось, что неосесимметрично нагруженные цилиндры при некоторых параметрах на современных компьютерах уже можно решать в рамках предложенного метода «переноса краевых условий» совсем без применения операций построчного ортонормирования. Это, например, при параметрах цилиндра L/R=2 и R/h=100. При параметрах цилиндра L/R=2 и R/h=200 все же оказываются необходимыми процедуры ортонормирования. Но на современных персональных компьютерах уже не наблюдаются сплошные скачки значений от больших отрицательных до больших положительных по всему интервалу между краями цилиндра - как это было на компьютерах поколения 286. В частном случае L/R=2 и R/h=200 наблюдаются лишь незначительные скачки в районе максимума изгибающего обезразмеренного момента М1 на левом крае и небольшой скачек обезразмеренного момента М1 на правом крае. Приводятся графики изгибающего обезразмеренного момента М1: - слева приводятся графики, полученные при использовании операций построчного ортонормирования на каждом из 100 шагов, на которые разделялся участок между краями цилиндра, - справа приводятся графики, полученные совсем без применения операций построчного ортонормирования. Следует сказать, что в качестве расчетной среды использовалась 32-х битная операционная система WindowsXP и среда программирования MicrosoftVisualStudio 2010 (VisualC++) использовалась в тех же рамках 32-х битной организации операций с числами. Параметры компьютера такие: ноутбук ASUSM51V (CPUDuoT5800). Компьютеры будут и дальше развиваться такими же темпами как сейчас и это означает, что в самое ближайшее время для подобных расчетов типа расчета неосесимметрично нагруженных оболочек вращения совсем не потребуется применять ортонормирование в рамках предложенного метода «переноса краевых условий», что существенно упрощает программирование метода и увеличивает скорость расчетов не только по сравнению с другими известными методами, но и по сравнению с собственными характеристиками метода «переноса краевых условий» предыдущих лет. СПИСОК ЛИТЕРАТУРЫ P.S. Метод Алексея Юрьевича Виноградова «переноса краевых условий» первоначально был опубликован в Межвузовском сборнике МИРЭА (кажется в 1995 году). МИРЭА это Московский институт радиотехники, электроники и автоматики. Точное название и год выхода стать можно посмотреть в Ленинской библиотеке в списке литературы моей диссертации. Там у меня только одна статья в МИРЭА. К сожалению, на руках у меня нет экземпляра моей кандидатской диссертации, поэтому не могу привести точное название статьи, но называется она, кажется, что-то вроде «Метод приведения краевых задач к задаче Коши».
,
– искомая вектор-функция задачи размерности 8х1, 
– производная искомой вектор-функции размерности 8х1, 
– квадратная матрица коэффициентов дифференциального уравнения размерности 8х8, 
– вектор-функция внешнего воздействия на систему размерности 8х1.
– значение искомой вектор-функции на левом крае х=0 размерности 8х1, 
– прямоугольная горизонтальная матрица коэффициентов краевых условий левого края размерности 4х8, 
– вектор внешних воздействий на левый край размерности 4х1,
– значение искомой вектор-функции на правом крае х=1 размерности 8х1, 
– прямоугольная горизонтальная матрица коэффициентов краевых условий правого края размерности 4х8, 
– вектор внешних воздействий на правый край размерности 4х1.=const, решение задачи Коши имеет вид [1]:
,
, где 
- это единичная матрица.
.
,
это вектор частного решения неоднородной системы дифференциальных уравнений.
, решение задачи Коши предлагается искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются:,
, где 
.
.
,
,
,
,
,
,
=const.
может рассматриваться на участке 
приближенно в виде постоянной величины 
, что позволяет вынести его из под знака интеграла, что приводит к совсем простому ряду для вычислений на рассматриваемом участке.
коэффициентов системы дифференциальных уравнений.на участке приближенно в виде постоянной величины , что позволяет вынести этот вектор из под знаков интегралов:
.:
на рассматриваемом участке 
через вектора частного решения 
, 
, 
соответствующих подучастков 
, 
, 
...
,
.
в 
и получим:
.
.
не просто складываются друг с другом, а с участием матрицы Коши подучастка.
и подставим сюда формулу для и получим:
на рассматриваемом участке через вычисленные частные вектора , , соответствующих подучастков , , ..
.в краевые условия левого края и получаем:
,
,
.
:
,
и
.
в перенесенные краевые условия точки :
,
,
.
:
,
и
.
переносим матричное краевое условие с левого края и таким же образом переносим матричное краевое условие с правого края.
в краевые условия правого края и получаем:
,
,
:
,
и
.
в перенесенные краевые условия точки :,
,
.
:
,
и
.интервала интегрирования переносим матричное краевое условие, как показано, и с левого края и таким же образом переносим матричное краевое условие с правого края и получаем:
,
.
.
..
,
.:,
и
.
в рассматриваемой точке :
.
.
и
,
,
- единичная матрица,
и 
взаимообратны.
,
..
=
.системы как векторы:
=(
,
,…,
).=делим на 
.
+
+…+
=
,
=(,,…,),
=
, =
, 
=1.
+
+…+
=
,
=(,,…,),
=
, =
,
=
-(
,),
=
-(,).
+
+…+
=
,
=(,,…,),
=
, =
,
=
-(,)-(,)-…-(,
)
,
=
-(,)-(,)-…-(,)
.,) – это скалярное произведение вектора на вектор и т.п.
, где матрица 
будет с ортонормированными строками, то есть обладает свойством 
, где - это единичная матрица.