Методы решения жестких краевых задач, включая новые методы и программы на С++ для реализации приведенных методов

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;

}

5. Второй вариант метода «переноса краевых условий» в произвольную точку интервала интегрирования.

Этот вариант метода еще не обсчитан на компьютерах.

Предложено выполнять интегрирование по формулам теории матриц [Гантмахер] сразу от некоторой внутренней точки интервала интегрирования к краям:

Y(0) = K(0<x) • Y(x) + Y*(0<x) ,

Y(1) = K(1<x) • Y(x) + Y*(1<x) .

Подставим эти формулы в краевые условия и получим:

U•Y(0) = u,

U•[ K(0<x) • Y(x) + Y*(0<x) ] = u,

[ U• K(0<x) ] • Y(x) = u - U•Y*(0<x) .

и

V•Y(1) = v,

V•[ K(1<x) • Y(x) + Y*(1<x) ] = v,

[ V• K(1<x) ] • Y(x) = v - V•Y*(1<x) .

То есть получаем два матричных уравнения краевых условий, перенесенные в рассматриваемую точку x:

[ U• K(0<x) ] • Y(x) = u - U•Y*(0<x) ,

[ V• K(1<x) ] • Y(x) = v - V•Y*(1<x) .

Эти уравнения аналогично объединяются в одну систему линейных алгебраических уравнений с квадратной матрицей коэффициентов для нахождения решения Y(x) в любой рассматриваемой точке x:

 • Y(x) = .

В случае «жестких» дифференциальных уравнений предлагается следующий алгоритм.

Используем свойство перемножаемости матриц Коши:

K(x<x) = K(x<x) • K(x<x) • … • K(x<x) • K(x<x)

и запишем выражения для матриц Коши, например, в виде:

K(0<x) = K(0<x) • K(x<x) • K(x<x),

K(1<x) = K(1<x) • K(x<x) • K(x<x) • K(x<x),

Тогда перенесенные краевые условия можно записать в виде:

[ U• K(0<x) • K(x<x) • K(x<x) ] • Y(x) = u - U•Y*(0<x) ,

[ V• K(1<x) • K(x<x) • K(x<x) • K(x<x) ] • Y(x) = v - V•Y*(1<x)

или в виде:

[ U• K(0<x) • K(x<x) • K(x<x) ] • Y(x) = u* ,

[ V• K(1<x) • K(x<x) • K(x<x) • K(x<x) ] • Y(x) = v* .

Тогда рассмотрим левое перенесенное краевое условие:

[ U• K(0<x) • K(x<x) • K(x<x) ] • Y(x) = u* ,

[ U• K(0<x) ] • { K(x<x) • K(x<x) • Y(x) } = u* ,

[ матрица ] • { вектор } = вектор .

Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[ U• K(0<x) ] • { K(x<x) • K(x<x) • Y(x) } = u* .

Далее последовательно можно записать:

[[ U• K(0<x) ] • K(x<x) ] • { K(x<x) • Y(x) } = u* ,

[ матрица ] • { вектор } = вектор .

Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[[ U• K(0<x) ] • K(x<x) ] • { K(x<x) • Y(x) } = u* ,

Далее аналогично можно записать:

[[[ U• K(0<x) ] • K(x<x) ] • K(x<x) ] • { Y(x) } = u*, [ матрица ] • { вектор} = вектор .

Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[[[ U• K(0<x) ] • K(x<x) ] • K(x<x) ] • Y(x) = u*

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

Далее проортонормированные уравнения краевых условий:

[ U• K(0<x) ] • Y(x) = u* ,

[ V• K(1<x) ] • Y(x) = v*

как и ранее объединяются в одну обычную систему линейных алгебраических уравнений с квадратной матрицей коэффициентов для нахождения искомого вектора Y(x) :

• Y(x) = .

6. Метод дополнительных краевых условий

Этот метод еще не обсчитан на компьютерах.

Запишем на левом крае ещё одно уравнение краевых условий:

M • Y(0) = m

В качестве строк матрицы M можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры краевых условий левого края L или линейно независимы с ними. Это вполне возможно, так как у краевых задач столько независимых физических параметров какова размерность задачи, а в параметры краевых условий входит только половина физических параметров задачи. То есть, например, если рассматривается задача об оболочке ракеты, то на левом крае могут быть заданы 4 перемещения. Тогда для матрицы М можно взять параметры сил и моментов, которых тоже 4, так как полная размерность такой задачи - 8. Вектор m правой части неизвестен и его надо найти и тогда можно считать, что краевая задача решена, то есть сведена к задаче Коши, то есть найден вектор Y(0) из выражения:

• Y(0) = ,

то есть вектор Y(0) находится из решения системы линейных алгебраических уравнений с квадратной невырожденной матрицей коэффициентов, состоящей из блоков U и M.

Аналогично запишем на правом крае ещё одно уравнение краевых условий:

N • Y(0) = n ,

где матрица N записывается из тех же соображений дополнительных линейно независимых параметров на правом крае, а вектор n неизвестен.

Для правого края тоже справедлива соответствующая система уравнений:

• Y(1) = .

Запишем Y(1) = K(1<0) •Y(0) + Y*(1<0) и подставим в последнюю систему линейных алгебраических уравнений:

• [ K(1<0) •Y(0) + Y*(1<0) ] = ,

• K(1<0) •Y(0) = - • Y*(1<0),

• K(1<0) •Y(0) = ,

• K(1<0) •Y(0) = .

Запишем вектор Y(0) через обратную матрицу:

Y(0) = •

и подставим в предыдущую формулу:

• K(1<0) • • = .

Таким образом, мы получили систему уравнений вида:

В • = ,

где матрица В известна, векторы u и s известны, а векторы m и t неизвестны.

Разобьем матрицу В на естественные для нашего случая 4 блока и получим:

• = ,

откуда можем записать, что

В11 • u + B12 • m = s,

B21 • u + B22 • m = t.

Следовательно, искомый вектор m вычисляется по формуле:

m = B12 • (s - B11• u).

А искомый вектор n вычисляется через вектор t:

t = B21 • u + B22 • m,

n = t + N • Y*(1<0).

В случае «жестких» дифференциальных уравнений предлагается выполнять поочередное построчное ортонормирование.

Запишем приведенную выше формулу

 • K(1<0) • • =

в виде:

• K(1<x2) • K(x2<x1) • K(x1<0) • • = .

Эту формулу можно записать в виде разделения левой части на произведение матрицы на вектор:

[ • K(1<x2) ] • { K(x2<x1) • K(x1<0) • • } =

[ матрица ] • { вектор } = вектор

Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[ • K(1<x2) ] • { K(x2<x1) • K(x1<0) • • } =

Здесь следует сказать, что подвектор t подвергать преобразованию не нужно, так как невозможно, так как его первоначальное значение не известно. Но подвектор t нам оказывается и не нужен для решения задачи.

Далее запишем:

[[ • K(1<x2) ] • K(x2<x1)] • { K(x1<0) • • } =

[ матрица ] • { вектор } = вектор

Аналогично и эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[[ • K(1<x2) ] • K(x2<x1)] • { K(x1<0) • • } = .

И так далее.

В результате поочередного ортонормирования получим:

В • = ,

• = .

Следовательно, искомый вектор m вычисляется по формуле:

m = B12 • (s - B11• u).

7. Формула для начала счета методом прогонки С.К. Годунова

Эта формула обсчитана на компьютерах в кандидатской диссертации.

Рассмотрим проблему метода прогонки С.К.Годунова.

Предположим, что рассматривается оболочка ракеты. Это тонкостенная труба. Тогда система линейных обыкновенных дифференциальных уравнений будет 8-го порядка, матрица A(x) коэффициентов будет иметь размерность 8х8, искомая вектор-функция Y(x) будет иметь размерность 8х1, а матрицы краевых условий будут прямоугольными горизонтальными размерности 4х8.

Тогда в методе прогонки С.К.Годунова для такой задачи решение ищется в следующем виде:

Y(x) = Y(x) c + Y(x) c + Y(x) c + Y(x) c + Y*(x),

или можно записать в матричном виде:

Y(x) = Y(x) • c + Y*(x),

где векторы Y(x), Y(x), Y(x), Y(x) - это линейно независимые вектора-решения однородной системы дифференциальных уравнений, а вектор Y*(x) - это вектор частного решения неоднородной системы дифференциальных уравнений.

Здесь Y(x)=|| Y(x), Y(x), Y(x), Y(x) || это матрица размерности 8х4, а c это соответствующий вектор размерности 4х1из искомых констант c,c,c,c.

Но вообще то решение для такой краевой задачи с размерностью 8 (вне рамок метода прогонки С.К.Годунова) может состоять не из 4 линейно независимых векторов Y(x), а полностью из всех 8 линейно независимых векторов-решений однородной системы дифференциальных уравнений:

Y(x)=Y(x)c+Y(x)c+Y(x)c+Y(x)c+

+Y(x)c+Y(x)c+Y(x)c+Y(x)c+Y*(x),

И как раз трудность и проблема метода прогонки С.К.Годунова и состоит в том, что решение ищется только с половиной возможных векторов и констант и проблема в том, что такое решение с половиной констант должно удовлетворять условиям на левом крае (стартовом для прогонки) при всех возможных значениях констант, чтобы потом найти эти константы из условий на правом крае.

То есть в методе прогонки С.К.Годунова есть проблема нахождения таких начальных значений Y(0), Y(0), Y(0), Y(0), Y*(0) векторов Y(x), Y(x), Y(x), Y(x), Y*(x), чтобы можно было начать прогонку с левого края x=0, то есть чтобы удовлетворялись условия U•Y(0) = u на левом крае при любых значениях констант c,c,c,c.

Обычно эта трудность «преодолевается» тем, что дифференциальные уравнения записываются не через функционалы, а через физические параметры и рассматриваются самые простейшие условия на простейшие физические параметры, чтобы начальные значения Y(0), Y(0), Y(0), Y(0), Y*(0) можно было угадать. То есть задачи со сложными краевыми условиями так решать нельзя: например, задачи с упругими условиями на краях.

Ниже предлагается формула для начала вычислений методом прогонки С.К.Годунова.

Выполним построчное ортонормирование матричного уравнения краевых условий на левом крае:

U•Y(0) = u,

где матрица U прямоугольная и горизонтальная размерности 4х8.

В результате получим эквивалентное уравнение краевых условий на левом крае, но уже с прямоугольной горизонтальной матрицей U размерности 4х8, у которой будут 4 ортонормированные строки:

U•Y(0) = u,

где в результате ортонормирования вектор u преобразован в вектор u.

Как выполнять построчное ортонормирование систем линейных алгебраических уравнений можно посмотреть в [Березин, Жидков].

Дополним прямоугольную горизонтальную матрицу U до квадратной невырожденной матрицы W:

W = ,

где матрица М размерности 4х8 должна достраивать матрицу U до невырожденной квадратной матрицы W размерности 8х8.

В качестве строк матрицы М можно взять те краевые условия, то есть выражения тех физических параметров, которые не входят в параметры левого края или линейно независимы с ними. Это вполне возможно, так как у краевых задач столько независимых физических параметров какова размерность задачи, то есть в данном случае их 8 штук и если 4 заданы на левом крае, то ещё 4 можно взять с правого края.

Завершим ортонормирование построенной матрицы W, то есть выполним построчное ортонормирование и получим матрицу W размерности 8х8 с ортонормированными строками:

W = .

Можем записать, что

Y(0) = (М)транспонированная = М.

Тогда, подставив в формулу метода прогонки С.К.Годунова, получим:

Y(0) = Y(0) •с + Y*(0)

или

Y(0) = М•с + Y*(0).

Подставим эту последнюю формулу в краевые условия U•Y(0) = u и получим:

U• [ М•с + Y*(0) ]= u.

Отсюда получаем, что на левом крае константы c уже не на что не влияют, так как

U• М = 0 и остается только найти Y*(0) из выражения:

U• Y*(0) = u.

Но матрица U имеет размерность 4х8 и её надо дополнить до квадратной невырожденной, чтобы найти вектор Y*(0) из решения соответствующей системы линейных алгебраических уравнений:

• Y*(0) = ,

где 0 - любой вектор, в том числе вектор из нулей.

Отсюда получаем при помощи обратной матрицы:

Y*(0) = • ,

Тогда итоговая формула для начала вычислений методом прогонки С.К.Годунова имеет вид:

Y(0) = М•с + • .

8. Второй алгоритм для начала счета методом прогонки С.К. Годунова

Этот алгоритм обсчитан на компьютерах в кандидатской диссертации.

Этот алгоритм требует дополнения матрицы краевых условий U до квадратной невырожденной:

Начальные значения Y(0), Y(0), Y(0), Y(0), Y*(0) находятся из решения следующих систем линейных алгебраических уравнений:

• Y*(0) = ,

• Y(0) = , где i = , , , ,

где 0 - вектор из нулей размерности 4х1.

9. Замена метода численного интегрирования Рунге-Кутта в методе прогонки С.К. Годунова

Эта замена формул Рунге-Кутта на формулу теории матриц обсчитана на компьютерах в кандидатской диссертации.

В методе С.К.Годунова как показано выше решение ищется в виде:

Y(x) = Y(x) • c + Y*(x)

На каждом конкретном участке метода прогонки С.К.Годунова между точками ортогонализации можно вместо метода Рунге-Кутта пользоваться теорией матриц и выполнять расчет через матрицу Коши:

Y(x) = K(x- x) •Y(x)

Так выполнять вычисления быстрее, особенно для дифференциальных уравнений с постоянными коэффициентами.

И аналогично через теорию матриц можно вычислять и вектор Y*(x) частного решения неоднородной системы дифференциальных уравнений. Или для этого вектора отдельно можно использовать метод Рунге-Кутта, то есть можно комбинировать теорию матриц и метод Рунге-Кутта.

10. Метод половины констант

Этот метод пока не обсчитан на компьютерах.

Выше было показано, что решение системы линейных обыкновенных дифференциальных уравнений можно искать в виде только с половиной возможных векторов и констант. Была приведена формула для начала вычислений:

Y(0) = М•с + • .

Из теории матриц известно, что если матрица ортонормирована, то её обратная матрица есть её транспонированная матрица. Тогда последняя формула приобретает вид:

Y(0) = М•с + U•u

или

Y(0) = U•u + М•с

или

Y(0) = • ,

Таким образом записана в матричном виде формула для начала счета с левого края, когда на левом крае удовлетворены краевые условия.

Далее запишем V•Y(1) = v и Y(1) = K(1<0) •Y(0) + Y*(1<0) совместно:

V• [ K(1<0) •Y(0) + Y*(1<0) ] = v

V• K(1<0) •Y(0) = v - V•Y*(1<0)

и подставим в эту формулу выражение для Y(0):

V• K(1<0) • • = v - V•Y*(1<0).

V• K(1<0) • • = p.

Таким образом мы получили выражение вида:

D • = p,

где матрица D имеет размерность 4х8 и может быть естественно представлена в виде двух квадратных блоков размерности 4х4:

• = p.

Тогда можем записать:

D1• u + D2 • c = p.

Отсюда получаем, что:

c = D2 • ( p - D1• u )

Таким образом, искомые константы найдены.

Далее показано как применять этот метод для решения «жестких» краевых задач.

Запишем

V• K(1<0) • • = p.

совместно с K(1<0) = K(1<x2) • K(x2<x1) • K(x1<0) и получим:

V• K(1<x2) • K(x2<x1) • K(x1<0) • • = p.

Эту систему линейных алгебраических уравнений можно представить в виде:

[ V• K(1<x2) ] • { K(x2<x1) • K(x1<0) • • } = p.

[ матрица ] • { вектор } = вектор

Эту группу линейных алгебраических уравнений можно подвергнуть построчному ортонормированию, которое сделает строчки [матрицы] ортонормированными, {вектор} затронут не будет, а вектор получит преобразование. То есть получим:

[V• K(1<x2) ] • { K(x2<x1) • K(x1<0) • • } = p.

И так далее.

В итоге поочередного вычленений матриц слева из вектора и ортонормирования получим систему:

D • = p,

Отсюда получаем, что:

c = D2 • (p - D1• u)

Таким образом, искомые константы найдены.

11. Применяемые формулы ортонормирования

Эти формулы обсчитаны в кандидатской диссертации.

Взято из: Березин И.С., Жидков Н.П. Методы вычислений, том II, Государственное издательство физико-математической литературы, Москва, 1962 г. 635 стр.

Пусть дана система линейных алгебраических уравнений порядка n:

А=.

Здесь над векторами поставим черточки вместо их обозначения жирным шрифтом.

Будем рассматривать строки матрицы А системы как векторы:

=(,,…,).

Ортонормируем эту систему векторов.

Первое уравнение системы А= делим на .

При этом получим:

++…+=, =(,,…,),

где =, =, =1.

Второе уравнение системы заменяется на:

++…+=, =(,,…,),

где =, =,

=-(,), =-(,).

Аналогично поступаем дальше. Уравнение с номером i примет вид:

++…+=, =(,,…,),

где =, =,

=-(,)-(,)-…-(,),

=-(,)-(,)-…-(,).

Процесс будет осуществим, если система линейных алгебраических уравнений линейно независима.

В результате мы придем к новой системе С=, где матрица С будет с ортонормированными строками, то есть обладает свойством С*С= E, где Е - это единичная матрица.

(Таким образом, решение системы можно записать в виде = С.)

12. Вывод формул, позаимствованный из «Теории матриц» Гантмахера

Система линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами имеет вид:

Y(x) = A Y(x) + F(x). (1)

Разложим Y(x) в ряд Маклорена по степеням x:

Y(x)=Y + Yx + Yx/2! + …, где Y=Y(0), Y= Y(0), … (2).

Из (1) почленным дифференцированием при А=const и F(x)=0 получим:

Y= AY= AY, Y= A Y = AY, (3)

Положив в (3) x=0 и подставив в (2) получим:

Y(x) = Y + Ax Y + A x/2! Y + … = e Y, (4)

где e = E + Ax + A x/2! + …, где Е - единичная матрица. (5)

Если принять x=x, то (4) заменится на

Y(x) = e Y(x), (6)

Рассмотрим случай A=const и F?0.

Введем в рассмотрение вектор-функцию Ya(x) в виде:

Y(x)= eYa(x). (7)

Продиффренцируем (7) и подставим в (1). Получим:

eYa(x) = F(x). (8)

При получении (8) учитывалось, что:

= = A + A x + A x/2! + … = A e.

Из (8) следует, что:

Ya(x) = c + . (9)

Подставим в (7) и получаем:

Y(x) = ec + e. (10)

Положив x=x в (10) получим: c = e Y(x). (11)

Окончательно получаем:

Y(x) = e Y(x) + e. (12)

13. P.P.S. Метод для численного интегрирования дифференциальных уравнений

Далее идёт картинка с текстом и с рисунками численного интегрирования:

Читали нам как-то в бауманке численные методы решения дифференциальных уравнений. И, кажется, приводили аналитический вывод формул одного из авторов. Или это просто мелькнуло в учебнике (я имею в виду вывод формул). Уже не очень помню. Запомнилась только собственная мысль, что людям вообще-то проще всего даются геометрические аналогии и выводы, сделанные на основе понятных геометрических картинок. Ну, вот тогда я и нарисовал один из вариантов численного решения дифференциальных уравнений и помню даже перевёл геометрические картинки в буквенные формулы приближённых вычислений. Сейчас повторно выводить буквенные формулы для численного интегрирования дифференциальных уравнений мне не кажется интересным. А вот привести картинки тех студенческих мыслей вполне можно для обсуждения.

14. Программа на С++ расчета цилиндра

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

Без ортонормирования в методе переноса краевых условий А.Ю.Виноградова успешно решается задача, например, нагружения цилиндрической оболочки, которая консольно заделана по правому краю и нагружена по левому краю силой, равномерно распределенной по дуге окружности, с отношением длинны к радиусу 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;

}

Размещено на Allbest.ru