Исправленные методы А.Ю. Виноградова: решения краевых задач, в том числе жестких краевых задач

Размещено на http://www.allbest.ru/

Исправленные методы А.Ю. Виноградова: решения краевых задач, в том числе жестких краевых задач

к.ф.-м.н. Виноградов А.Ю.

Введение

На примере системы дифференциальных уравнений цилиндрической оболочки ракеты - системы обыкновенных дифференциальных уравнений 8-го порядка (после разделения частных производных).

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

Y(x) = A(x) • Y(x) + F(x),

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

Здесь и далее вектора обозначаем жирным шрифтом вместо черточек над буквами

Краевые условия имеют вид:

U•Y(0) = u,

V•Y(1) = v,

Где Y(0) - значение искомой вектор-функции на левом крае х=0 размерности 8х1, U - прямоугольная горизонтальная матрица коэффициентов краевых условий левого края размерности 4х8, u - вектор внешних воздействий на левый край размерности 4х1,

Y(1) - значение искомой вектор-функции на правом крае х=1 размерности 8х1, V - прямоугольная горизонтальная матрица коэффициентов краевых условий правого края размерности 4х8, v - вектор внешних воздействий на правый край размерности 4х1.

В случае, когда система дифференциальных уравнений имеет матрицу с постоянными коэффициентами A=const, решение задачи Коши имеет вид [Гантмахер]:

Y(x) = e• Y(x) + e• e• F(t) dt,

где

e= E + A(x-x) + A (x-x)/2! + A (x-x)/3! + …,

где E это единичная матрица.

Матричная экспонента ещё может называться матрицей Коши или матрициантом и может обозначаться в виде:

K(x<x) = K(x - x) = e.

Тогда решение задачи Коши может быть записано в виде:

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

где

Y*(x<x) = e• e• F(t) dt

это вектор частного решения неоднородной системы дифференциальных уравнений

1. Случай переменных коэффициентов

Этот вариант рассмотрения переменных коэффициентов проверялся в кандидатской диссертации.

Из теории матриц [Гантмахер] известно свойство перемножаемости матричных экспонент (матриц Коши):

e= e• e • … • e • e,

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

В случае, когда система дифференциальных уравнений имеет матрицу с переменными коэффициентами A=A(x), решение задачи Коши предлагается искать при помощи свойства перемножаемости матриц Коши. То есть интервал интегрирования разбивается на малые участки и на малых участках матрицы Коши приближенно вычисляются по формуле для постоянной матрицы в экспоненте. А затем матрицы Коши, вычисленные на малых участках, перемножаются:

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

где матрицы Коши приближенно вычисляются по формуле:

K(x<x) = e, где ?x= x- x.

2. Формула для вычисления вектора частного решения неоднородной системы дифференциальных уравнений

Эта очень простая формула еще не обсчитана на компьютерах. Вместо неё обсчитывалась (а может быть и не обсчитывалась, не знаю) значительно ранее выведенная (выведенная моим отцом) и гораздо более сложная формула, приведенная в:

Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf

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

Y*(x<x) = e• e• F(t) dt

предлагается использовать следующую формулу для каждого отдельного участка интервала интегрирования:

Y*(x<x) = Y*(x- x) = K(x- x) •K(x- t) • F(t) dt .

Правильность приведенной формулы подтверждается следующим:

Y*(x- x) = e•e• F(t) dt ,

Y*(x- x) = e•e• F(t) dt ,

Y*(x- x) = e• F(t) dt ,

Y*(x- x) = e• F(t) dt ,

Y*(x- x) = e• e• F(t) dt ,

Y*(x<x) = e• e• F(t) dt,

что и требовалось подтвердить.

Вычисление вектора частного решения системы дифференциальных уравнений производиться при помощи представления матрицы Коши под знаком интеграла в виде ряда и интегрирования этого ряда поэлементно:

Y*(x<x) = Y*(x- x) = K(x- x) •K(x- t) • F(t) dt =

= K(x- x) • (E + A(x- t) + A (x- t)/2! + … ) • F(t) dt =

= K(x- x) • (EF(t) dt + A•(x- t) • F(t) dt + A/2! •(x- t) • F(t) dt + … ) .

Эта формула справедлива для случая системы дифференциальных уравнений с постоянной матрицей коэффициентов A=const.

Вектор F(t) может рассматриваться на участке (x- x) приближенно в виде постоянной величины F(х)=constant, что позволят вынести его из под знака интеграла, что приводит к совсем простому ряду для вычислений на рассматриваемом участке.

Для случая дифференциальных уравнений с переменными коэффициентами в приведенной выше формуле для каждого участка может использоваться осредненная матрица А: A=А(х) коэффициентов системы дифференциальных уравнений.

Приведем (итерационные или рекуррентные) формулы вычисления вектора частного решения, например, Y*(x<x) на рассматриваемом участке (x<x) через вектора частного решения Y*(x<x), Y*(x<x), Y*(x<x) соответствующих подучастков (x<x), (x<x), (x<x).

Имеем:

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

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

Y*(x<x) = Y*(x- x) = K(x- x) •K(x- t) • F(t) dt.

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

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

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

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

Y(x) = K(x<x) • [K(x<x) • Y(x) + Y*(x<x)] + Y*(x<x) =

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

Сравним полученное выражение с формулой:

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

и получим, очевидно, что

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

и, самое главное здесь - для частного вектора получаем формулу:

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

То есть вектора подучастков Y*(x<x) и Y*(x<x) не просто складываются друг с другом, а с участием матриц Коши подучастков.

Аналогично запишем:

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

И подставим сюда формулу для Y(x) и получим:

Y(x) = K(x<x) • [K(x<x) • K(x<x) • Y(x) + K(x<x) • Y*(x<x) + Y*(x<x)] + Y*(x<x) = K(x<x) • K(x<x) • K(x<x) • Y(x) + K(x<x) • K(x<x) • Y*(x<x) + K(x<x) • Y*(x<x) + Y*(x<x).

Сравнив полученное выражение с формулой:

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

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

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

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

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

3. Метод "переноса краевых условий" в произвольную точку интервала интегрирования

Метод обсчитан на компьютерах. По нему уже сделано 3 кандидатских физ-мат диссертации.

Метод подходит для любых краевых задач. А для "жестких" краевых задач показано, что метод считает быстрее, чем метод С.К.Годунова до 2-х порядков (в 100 раз), а для некоторых "жестких" краевых задач не требует ортонормирования вовсе. Численный метод переноса краевых условий для жестких дифференциальных уравнений строительной механики Журнал "ММ", Том: 14 (2002), Номер: 9, 3 стр. 1409-003r.pdf

Полное решение системы дифференциальных уравнений имеет вид

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

Или можно записать:

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

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

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) .

Или получаем краевые условия, перенесенные в точку x:

U• Y(x) = u ,

где

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

Далее запишем аналогично

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

И подставим это выражение для Y(x) в перенесенные краевые условия точки x

U• Y(x) = u,

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

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

Или получаем краевые условия, перенесенные в точку x:

U• Y(x) = u ,

где

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

Покажем перенос краевых условий с правого края.

Можно записать:

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

Подставляем это выражение для Y(1) в краевые условия правого края и получаем:

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:

V• Y(x) = v ,

где

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

Далее запишем аналогично

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

И подставим это выражение для Y(x) в перенесенные краевые условия точки x

V• Y(x) = v ,

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

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

Или получаем краевые условия, перенесенные в точку x:

V• Y(x) = v ,

где

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

И так в точку x переносим матричное краевое условие с левого края и таким же образом переносим матричное краевое условие с правого края и получаем:

U• Y(x) = u ,

V• Y(x) = v .

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

• Y(x) = .

А в случае "жестких" дифференциальных уравнений предлагается применять построчное ортонормирование матричных краевых условий в процессе их переноса в рассматриваемую точку. Для этого формулы ортонормирования систем линейных алгебраических уравнений можно взять в [Березин, Жидков]. То есть, получив

U• Y(x) = u,

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

U• Y(x) = u.

И теперь уже в это проортонормированное построчно уравнение подставляем

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

И получаем

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

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

Или получаем краевые условия, перенесенные в точку x:

U• Y(x) = u ,

где

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

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

U• Y(x) = u.

И так далее.

И аналогично поступаем с промежуточными матричными краевыми условиями, переносимыми с правого края в рассматриваемую точку.

В итоге получаем систему линейных алгебраических уравнений с квадратной матрицей коэффициентов, состоящую из двух независимо друг от друга поэтапно проортонормированных матричных краевых условий. Эта система решается методом Гаусса с выделением главного элемента для получения решения Y(x) в рассматриваемой точке x:

• Y(x) = .

4. Программа на С++ расчета цилиндрической оболочки

В качестве проверочных задач использовалась схема консольно закрепленных цилиндрической и сферической оболочек с параметрами R/h=50, 100, 200. Длина цилиндрической оболочки рассматривалась L/R=2, а угловые координаты сферической оболочки рассматривались от /4 до 3/4. На свободном крае рассматривалось нормальное к поверхности оболочек погонное усилие, равномерно распределенное в интервале [-/4, /4]. В качестве среды программирования использовалась система Microsoft Visual Studio 2010 (Visual C++).

Первоначально метод был предложен и обсчитывался в кандидатской диссертации А.Ю.Виноградова в 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-х битная операционная система Windows XP и среда программирования Microsoft Visual Studio 2010 (Visual C++) использовалась в тех же рамках 32-х битной организации операций с числами. Параметры компьютера такие: ноутбук ASUS M51V (CPU Duo T5800).

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

4.1 Программа на С++ (расчет цилиндра)

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

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

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

#include "stdafx.h"

#include <iostream>

#include <conio.h>

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

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

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

double A[8][8];

double b[8];

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

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

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

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

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

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

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

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

for(int i=jj;i<8;i++){//Находится номер е строки, где лежит главный (максимальный) элемент в столбце jj и делается взаимозамена строк

if ((A[i][jj]*A[i][jj])>s) {//Вместо перемножения (удаляется возможный знак минуса) можно было бы использовать функцию по модулю abs()

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

if (e<0) {

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

if (e>jj) {//Если главный элемент не в строке с номером jj. а в строке с номером е

main=A[e][jj];

for(int j=0;j<8;j++){//Взаимная замена двух строк - с номерами e и jj

t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t;

t=b[jj]; b[jj]=b[e]; b[e]=t;

for(int i=(jj+1);i<8;i++){//Приведение к верхнетреугольной матрице

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

A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//Перерасчет коэффициентов строки i>(jj+1)

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

A[i][jj]=0.0;//Обнуляемые элементы столбца под диагональным элементом матрицы А

}//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

x[8-1]=b[8-1]/A[8-1][8-1];//Первоначальное определение последнего элемента искомого решения х (7-го)

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

t=0;

for(int j=1;j<(8-i);j++){

t=t+A[i][i+j]*x[i+j];

x[i]=(1/A[i][i])*(b[i]-t);

return 0;

int main()

int nn;//Номер гармоники, начиная с 1-й (без нулевой)

int nn_last=50;//Номер последней гармоники

double Moment[100+1]={0};//Массив физического параметра (момента), что рассчитывается в каждой точке между краями

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

fclose(fp);

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

_getch();

return 0;

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

Программа на С++ (расчет сферы):

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

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 (m<n) {

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

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

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

//Вычисление матрицы MAT_ROW в виде матричного ряда для последующего использования

//при вычислении вектора partial_vector - вектора частного решения неоднородной системы ОДУ на шаге delta_x

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

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

int n=100, m;

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

int i,j,k;

//E - единичная матрица - первый член ряда MAT_ROW

E[0][0]=1.0; E[1][1]=1.0; E[2][2]=1.0; E[3][3]=1.0;

E[4][4]=1.0; E[5][5]=1.0; E[6][6]=1.0; E[7][7]=1.0;

//первоначальное заполнение вспомогательного массива TMP1 - предыдущего члена ряда для следующего перемножения

//и первоначальное заполнение MAT_ROW первым членом ряда

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

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

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

MAT_ROW[i][j]=E[i][j];

//ряд вычисления MAT_ROW, начиная со 2-го члена ряда (m=2;m<=n)

for(m=2;m<=n;m++) {

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

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

TMP2[i][j]=0;

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

TMP2[i][j]+=TMP1[i][k]*A[k][j];

TMP2[i][j]*=delta_x;

TMP2[i][j]/=m;

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

//заполнение вспомогательного массива TMP1 для вычисления следующего члена ряда - TMP2 в следующем шаге цикла по m

if (m<n) {

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

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

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

//Вычисление vector - НУЛЕВОГО (частный случай) вектора частного решения

//неоднородной системы дифференциальных уравнений на рассматриваемом участке:

void partial_vector(double vector[8]){

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

vector[i]=0.0;

//Вычисление vector - вектора частного решения неоднородной системы дифференциальных уравнений на рассматриваемом участке delta_x:

void partial_vector_real(double expo_[8][8], double mat_row[8][8], double x_, double delta_x, double vector[8]){

double POWER_[8]={0};//Вектор внешней нагрузки на оболочку

double REZ[8]={0};

double REZ_2[8]={0};

power_vector_for_partial_vector(x_, POWER_);//Расчитываем POWER_ при координате x_

mat_on_vect(mat_row, POWER_, REZ);//Умножение матрицы mat_row на вектор POWER_ и получаем вектор REZ

mat_on_vect(expo_, REZ, REZ_2);//Умножение матрицы expo_ на вектор REZ и получаем вектор REZ_2

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

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

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

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

double A[8][8];

double b[8];

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

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

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

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

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

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

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

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

for(int i=jj;i<8;i++){//Находится номер е строки, где лежит главный (максимальный) элемент в столбце jj и делается взаимозамена строк

if ((A[i][jj]*A[i][jj])>s) {//Вместо перемножения (удаляется возможный знак минуса) можно было бы использовать функцию по модулю abs()

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

if (e<0) {

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

if (e>jj) {//Если главный элемент не в строке с номером jj. а в строке с номером е

main=A[e][jj];

for(int j=0;j<8;j++){//Взаимная замена двух строк - с номерами e и jj

t=A[jj][j]; A[jj][j]=A[e][j]; A[e][j]=t;

t=b[jj]; b[jj]=b[e]; b[e]=t;

for(int i=(jj+1);i<8;i++){//Приведение к верхнетреугольной матрице

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

A[i][j]=A[i][j]-(1/main)*A[jj][j]*A[i][jj];//Перерасчет коэффициентов строки i>(jj+1)

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

A[i][jj]=0.0;//Обнуляемые элементы столбца под диагональным элементом матрицы А

}//Цикл по столбцам jj преобразования матрицы А в верхнетреугольную

x[8-1]=b[8-1]/A[8-1][8-1];//Первоначальное определение последнего элемента искомого решения х (7-го)

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

t=0;

for(int j=1;j<(8-i);j++){

t=t+A[i][i+j]*x[i+j];

x[i]=(1/A[i][i])*(b[i]-t);

return 0;

int main()

int nn;//Номер гармоники, начиная с 1-й (без нулевой)

int nn_last=50;//Номер последней гармоники

double Moment[100+1]={0};//Массив физического параметра (момента), что рассчитывается в каждой точке между краями

double angle;

double start_angle, finish_angle;

start_angle=3.14159265359/4;

finish_angle=start_angle+(3.14159265359/2);

double step=(3.14159265359/2)/100; //step=(3.14159265359/2)/100 - величина шага расчета оболочки - шага интервала интегрирования (должна быть больше нуля, т.е. положительная)

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

h_div_R=1.0/200;

double c2;

c2=h_div_R*h_div_R/12;//Величина h*h/R/R/12

double nju;

nju=0.3;

double gamma;

gamma=3.14159265359/4;//Угол распределения силы по левому краю

//распечатка в файлы:

FILE *fp;

// Open for write

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

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

else

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

...