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,"%f\n",Moment[ii]);

}

fclose(fp);

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

_getch();

return 0;

}

4.3. Программа на С++ расчета сферической оболочки (переменные коэффициенты).

ПРОГРАММА НА С++ (РАСЧЕТ СФЕРЫ):

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

//Решение краевой задачи с переменными коэффициентами - сфера.

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

#include "stdafx.h"

#include <iostream>

#include <conio.h>

#include <math.h> //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;