Основы высшей математики для инженеров 2009

Как видим, все эти три процедуры являются циклическими. При этом цикл вычисления прогоночных коэффициентов называют пря мой прогонкой (Direct race), а цикл вычисления значений искомой функции обратной прогонкой (Return race).

Взаключение следует отметить, что в пределах нескольких стра ниц мы изложили метод решения краевых задач с дифференциальным уравнением второго порядка с переменными коэффициентами. Фор тран программа практической реализации решений займет еще одну страницу. Данный пример есть наглядная иллюстрация того, какие ог ромные возможности дает нам современный компьютер.

Вэтой связи интересно также отметить, что решению одного толь ко дифференциального уравнения Бесселя

в свое время были посвящены отдельные монографии.

7.8.3. Фортран>программаM62решения краевых задач для дифференциальных

уравнений второго порядка методом прогонки

Данная программа написана как Difur.f90, но Фортран 77 обладает всеми возможностями, которые здесь использованы. Если потребуется оформить ее как Difur.for, т. е. с расширением .for, то для этого потре буется только заменить символ комментария, символы продолжения строки и сдвинуть некоторые операторы вправо, которые написаны не

с7 й позиции, а выходят за пределы поля операторов.

Кчислу исходных данных здесь отнесены:

N — число узловых точек;

r —длина отрезка ?0, l@, на котором строится решение; g —коэффициенты граничных условий.

Причем ввод исходных данных осуществляется оператором READ из внешнего файла с файловым указателем (ФУ) 11.

Коэффициенты уравнения p(x), q(x) и правая часть f (x) вычисляют ся как заданные функции. Эти программные единицы находятся в конце программы. В качестве тестового варианта заданы: p 0, q 1, f 1 при h1 1; g1 0; e1 0; h2 0; g2 1; e2 2. Тестовое уравнение решается при всех других вариантах расчета, и результаты его решения записываются во внешний файл с ФУ 12 в качестве второй таблицы результатов, где первая таблица соответствует заданному варианту ос новной программы.

7.8. Линейное уравнение второго порядка с переменными коэффицинтами... 343

Таким образом, чтобы получить решение заданного уравнения, достаточно вписать в программу выражения заданных функций p(x), q(x),f(x) и заполнить две строки файла исходных данных «ZADANO»:

OPEN(11,FILE=’ZADANO.txt’)

OPEN(12,FILE=’RESULT_Y.txt’)

!********************************************************M62

READ(11,*)N,eL,h1,g1,e1,h2,g2,e2 !* WRITE(*,’(2x,A,I5,2x,A,F4.0,/,&

2x,A,F4.0,2x,A,F4.0,2x,A,F4.0,/,&

2x,A,F4.0,2x,A,F4.0,2x,A,F4.0)’),& ‘N=’,N,’eL=’,eL, ‘h1=’,h1,’g1=’,g1,’e1=’,e1,& ‘h2=’,h2,’g2=’,g2,’e2=’,e2 !* WRITE(12,’(2x,A,I5,2x,A,F4.0,/,&

2x,A,F4.0,2x,A,F4.0,2x,A,F4.0,/,&

2x,A,F4.0,2x,A,F4.0,2x,A,F4.0)’),& ‘N=’,N,’eL=’,eL, ‘h1=’,h1,’g1=’,g1,’e1=’,e1,& ‘h2=’,h2,’g2=’,g2,’e2=’,e2

!*********************************************************

ALLOCATE (S(N),r(N),Xk(N),Yk(N))

dx=eL/(N-1);A1=g1-3.*h1*0.5/dx;B1=-2.*h1/dx;C1=-0.5*h1/dx d1=e1;x=dx;Xk(1)=0.

Ak=1.-p(x)*DX/2.; Ck=1+p(x)*DX/2.; Bk=2.+q(x)*DX**2 dk=f(x)*DX**2

S(1)=(B1-C1*Bk/Ck)/(A1-C1*Ak/Ck) r(1)=(d1-C1*dk/Ck)/(A1-C1*Ak/Ck)

DO K=2,N-1 Xk(K)=Xk(K-1)+dx; x=Xk(K)

Ak=1.-p(x)*DX/2.; Ck=1+p(x)*DX/2.; Bk=2.+q(x)*DX**2 dk=f(x)*DX**2

344 Глава 7. Обыкновенные дифференциальные уравнения

! Here Ab = AN+1; Bb = BN+1; Cb = CN+1 Ab=h2*0.5/dx; Bb=2.*h2/dx; Cb=g2+1.5*h2/dx; db=e2 Chisl=db+Bb*r(N-1)-Ab*r(N-2)-Ab*S(N-2)*r(N-1) Yk(N)=Chisl/(Ab*S(N-2)*S(N-1)-Bb*S(N-1)+Cb); Xk(N)=eL

DO K=1,N-1 Yk(N-K)=S(N-K)*Yk(N-K+1)+r(N-K)

Enddo! Reverse RACE m=(N-1)/10

WRITE(12,*)’It is Yk(x)’ write(12,’(2x,A,F16.4,3x,A,F16.4)’)&

((‘x=’,Xk(1+m*(K-1)),’y=’,Yk(1+m*(K-1))),k=1,11)

7.9. Системы линейных уравнений второго порядка с переменными коэффициентами. Метод матричной прогонки

Постановка задачи. Уравнения и граничные условия. Конечно разност ная аппроксимация уравнений и граничных условий. Решение конечно разностных уравнений методом матричной прогонки. Фортран подпро грамма вычисления обратной матрицы.

7.9.1. Постановка задачи. Уравнения и граничные условия

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

...........................................

где pij, qij , fi — заданные функции на отрезке ?0, l@.

Требуется найти решение данной системы уравнений на отрезке ?0, l@ при заданных условиях на границах.

Для более компактной записи данных уравнений и последующих преобразований введем вектор столбец Y с компонентами искомых функций, матрицы P, Q коэффициентов, вектор правых частей F и за

346 Глава 7. Обыкновенные дифференциальные уравнения

где

где Hi, Gi (i 1, 2) — квадратные матрицы п го порядка, составленные из коэффициентов заданных граничных условий; а e1, e2 — векторы правых частей.

7.9.2.Конечно>разностная аппроксимация уравнений и граничных условий

так же, как и одно дифференциальноеM62с его граничными условиями. Причем здесь мы имеем не только внешнее сходство этих двух задач, но и совершенно одинаковые алгоритмы их решения.

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

представленные в векторно матричной форме (1), (2), выглядят точно

В данном случае под словом алгоритм нужно иметь в виду после довательность операций, которые нужно выполнить для получения ре шения. Последовательность операций одинаковая. Разными являются лишь элементы, которыми здесь являются векторы и матрицы, и с ко торыми мы оперируем.

Для матриц и векторов тоже определены операции сложения и вы читания, операции умножения матрицы на матрицу и матрицы на вектор, а вычислению обратной величины 1b соответствует вычисле ние обратной матрицы B 1.

Поэтому дифференциальное уравнение (1) можно также записать в конечных разностях для всех внутренних узловых точек. Тогда при ис пользовании стандартной трехточечной схемы центральных разностей для первых и вторых производных получим

Совершенно аналогичные выражения можно получить и для пер вой производной на границах:

После некоторых преобразований граничных условий и внутрен них конечно разностных уравнений получаем систему (N 1) алгебраи ческих уравнений, которые совершенно аналогичны предыдущим (см. п. 7.8.2), но с матричными коэффициентами, для определения (N+1) неизвестных векторов Y1, Y2 , Y3, …, YN , YN 1:

AN 1YN 1 BN M621YN CN 1YN 1 dN 1,

где первая и последняя строка — граничные условия, а вторая строка это внутренние уравнения (3). Для коэффициентов этих уравнений получены следующие выражения:

k 2, 3, 4, …, N.

Здесь следует обратить внимание, что в коэффициенты Ak , Ck раз ностных уравнений входит матрица P дифференциального уравнения, она без индекса, но вычисляется при x xk . Аналогично в матрицу Bk разностных уравнений входит матрица Q(xk) дифференциального урав нения. Кроме того, здесь I — единичная матрица, dx — шаг между уз ловыми точками.

7.9.3. Решение конечно>разностных уравнений методом матричной прогонки

Для решения полученной системы уравнений (4) мы представим метод матричной прогонки, который аналогичен изложенному выше методу скалярной прогонки. Последовательность выполнения всех операций здесь повторяется полностью. И мы последовательно их пе реносим в этот более сложный алгоритм с учетом их преобразования применительно к векторам и матрицам.

Преобразуем первое граничное условие, исключив из него вектор y3 с помощью первого внутреннего уравнения, из которого можно по лучить

y3 C2 1(B2 y2 A2 y1 d2).

Подставим это выражение в первое граничное условие и преобра зуем его к виду

Теперь можно подставить y1 в первое внутреннее уравнение и по лучить аналогичное соотношение между двумя следующими перемен

ными:

y2 S2 y3 r2 ,

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

(9)

k 2, 3, …,N.

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

Осталось не использованным пока второе граничное условие, из которого с помощью двух последних уравнений (8) можно исключить

yN 1 SN 1(SN yN 1 rN ) rN 1

ивычислить

!(dN 1 BN 1rN AN 1rN 1 AN 1SN 1rN ).

После этого из уравнений (8) в обратной последовательности при k N, N 1, N 1, …, 2 вычисляются значения искомой функции во всех точках. Таким образом, алгоритм решения задачи в целом состо ит из следующих трех блоков:

вычисление коэффициентовM62разностных уравнений (5),(6);

вычисление первых прогоночных коэффициентов (7) и всех по следующих Sk , rk по формулам (9);

вычисление значения искомой функции (11) на границе и всех

значений yk с помощью уравнений (8).

Как видим, все эти три процедуры являются циклическими. При этом цикл вычисления прогоночных коэффициентов называют пря мой прогонкой (Direct race), а цикл вычисления значений искомой функции обратной прогонкой (Return race).

Все операции здесь нужно выполнять с векторами и матрицами. Сложение, вычитание и умножение матриц можно запрограммировать при использовании современного фортрана достаточно просто, даже не прибегая к использованию стандартных подпрограмм. Что же каса ется операции вычисления обратной матрицы, то здесь необходимо использовать подпрограмму. В связи с этим ниже приведена такая подпрограмма.

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

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

Возможность обобщения метода прогонки на метод матричной прогонки по представленной здесь схеме применительно к решению краевых задач теории оболочек и исследованию устойчивости была впервые опубликована в работах [21, 22]. При этом следует отметить, что в отдельных случаях для обеспечения возможности применения метода матричной прогонки, необходимо преодолеть некоторые труд ности. Они возникают, например, при конечно разностной аппрокси мации граничных условий и в связи с тем, что на одной границе все условия заданы для одной части вектора искомых функций, а на дру гой границе для другой части вектора, как, например, при решении краевых задач переноса тепла излучением [26].

7.9.4. Фортран>подпрограммаM62вычисления обратной матрицы

Вызов данной подпрограммы осуществляется оператором

CALL invertM(A,N,Det,V1,V2),

где А — данная матрица, обратная матрица возвращается в этот же массив; N — порядок матрицы А; Det — определитель матрицы А, ко торый вычисляется при выполнении операции обращения матрицы А; V1, V2 — рабочие одномерные массивы размерности N.

До вызова данной подпрограммы размерности массивов A, V1, V2 должны быть уже объявлены.

Ниже представлен текст данной подпрограммы вместе с примером ее использования для обращения матрицы 4 го порядка, элементы ко торой вводятся операторами присваивания.

Program Example_1

DIMENSION A(4,4),V1(4),V2(4)

OPEN (11,FILE=’A.txt’)

A(1,1)=1.; A(1,2)=2.; A(1,3)=3.; A(1,4)=4.

A(2,1)=0.; A(2,2)=2.; A(2,3)=3.; A(2,4)=4.

A(3,1)=0.; A(3,2)=0.; A(3,3)=3.; A(3,4)=4.

A(4,1)=0.; A(4,2)=0.; A(4,3)=0.; A(4,4)=4.

WRITE(11,’(2x,4F10.4)’)((A(i,j),j=1,4),I=1,4)

WRITE(*,*)’ It is given Matrix A’

WRITE(*,’(2x,4F10.4)’)((A(i,j),j=1,4),I=1,4)