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

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++){

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]

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

Подставим эту последнюю формулу в краевые условия 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• [ 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)

Литература

1. Гантмахер Ф.Р. Теория матриц. - М.: Наука, 1988. - 548 с.

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

3. www.vinogradov-design.narod.ru/math.html

4. www.vinogradov-best.narod.ru

5. www.AlexeiVinogradov.narod.ru

6. www.VinogradovAlexei.narod.ru

7. www.Vinogradov-Alexei.narod.ru

8. www.Vinogradov-math.narod.ru

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

...