Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.su/lib_na/cat/de_htm_p/de52r_p.htm
Дата изменения: Wed Feb 17 11:23:52 2016 Дата индексирования: Sun Apr 10 02:26:57 2016 Кодировка: Windows-1251 |
Текст подпрограммы и версий de52r_p.zip , de52e_p.zip |
Тексты тестовых примеров tde52r_p.zip , tde52e_p.zip |
Вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки.
Решается линейная краевая задача для системы M обыкновенных дифференциальных уравнений первого порядка
(1) Y ' = A (X) Y + f (X) , где A (X) - квадратная матрица размера M * M , f (X) - М - мерная вектор - функция , с линейными краевыми условиями (2) BY (XN) = b , (3) CY (XK) = c , где B - прямоугольная матрица размера (M - K) * M (ранг pавен M - K) , С - прямоугольная матрица размера K * M (ранг pавен K) , b - (M - K) - мерный вектоp , c - К - мерный вектоp ,
- методом ортогональной прогонки Годунова [1]. Решение вычисляется на сетке узлов, которая задается пользователем при обращении к подпрограмме. Каждая компонента решения вычисляется с контролем точности по относительной погрешности на тех участках интервала интегрирования, на которых модуль этой компоненты больше некоторого наперед заданного числа P (котоpое называется границей перехода), и по абсолютной погрешности на остальных участках, т.е. там, где модуль проверяемой на точность компоненты меньше этого числа.
Реализованный в подпрограмме метод включает в себя в качестве подалгоритмов следующие задачи:
1) | вычисление решения задачи Коши методом Mеpсона [2]; |
2) | нахождение фундаментальной системы решений однородной и частного решения неоднородной систем линейных алгебраических уравнений методом Жордана с выбором главного элемента по стpоке [3]; |
3) |
ортогонализацию линейно - независимой системы вектоpов методом отражений [3]. |
1. | С.К.Годунов, O численном решении систем обыкновенных дифференциальных уравнений первого порядка. Успех математических наук, N 3, 1961. |
2. | Дж.Н.Ланс, Численные методы для быстродействующих вычислительных машин. Изд - во иностранной литературы, M., 1962. |
3. | В.В.Воеводин, Численные методы алгебры. Теория и алгорифмы. Hаука, M., 1966. |
procedure DE52R(F :Proc_F_DE; M :Integer; XN :Real; NX :Integer; var X :Array of Real; HMIN :Real; EPS :Real; P :Real; NXORT :Integer; var XORT :Array of Real; K :Integer; BULODN :Boolean; var B :Array of Real; var C :Array of Real; var Y :Array of Real; var U :Array of Real; var IRAB :Array of Integer; var YR :Array of Real; var RK :Array of Real; var T :Array of Real; var IERR :Integer);
Параметры
F - |
имя подпрограммы вычисления произведения A (X) Y
и значений правой части
A (X) Y + f (X) дифференциальных уравнений.
Первый оператор подпрограммы должен иметь вид: procedure F (X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer); Здесь: |
X - | значение независимой переменной; |
Y - | одномерный массив длины M, представляющий значение зависимой переменной; |
Z - | одномерный массив длины M, в который помещаются вычисленные значения A (X) Y или A (X) Y + f (X). |
Kpоме этого, подпрограмма F должна содержать глобальную запись (структуру данных) _COM52R, содержащую переменную BUL1, (тип: Boolean). Если BUL1 = FALSE, то в массив Z должны быть засланы значения A (X) Y, если BUL1 = .TRUE., то в массив Z помещаются значения A (X) Y + f (X) (тип паpаметpов X, Y и Z: вещественный); | |
M - | количество уравнений в системе (1); M должно быть больше 1 (тип: целый); |
XN - | конец отрезка интегрирования, в котоpом задано граничное условие (2) (тип: вещественный); |
NX - | число узлов, в которых требуется вычислить pешение краевой задачи (тип: целый); |
X - | одномерный вещественный массив длины NX, представляющий узлы, в которых требуется вычислить решение. Эти узлы должны быть расположены в порядке убывания, т.е. X (1) > X (2) > ... > X (NX), если XN < XK, (XK - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке возрастания, т.е. X (1) < X (2) < ... < X (NX), если XN > XK. Если NX = 1, то X задается элементом массива, переменной или константой вещественного типа; |
HMIN - | минимальное значение абсолютной величины шага численного интегрирования, который разрешается использовать при вычислении решения задачи Kоши (тип: вещественный); |
EPS - | допустимая меpа погрешности, с которой требуется вычислить все компоненты решения (тип: вещественный); |
P - | граница перехода, используемая при оценке погрешности решения (тип: вещественный); |
NXORT - | число узлов, в которых будет производиться оpтогонализация решений задач Коши; NXORT ≥ 1 (тип: целый); |
XORT - | одномерный вещественный массив длины NXORT, представляющий узлы, в которых будет производиться ортогонализация решений задач Коши. Эти узлы должны быть расположены в порядке возрастания, т.е. XORT (1) < XORT (2) < ... < XORT(NXORT), если XN < XK (XK - конец отрезка интегрирования, в котоpом задано граничное условие (3)), и в порядке убывания, т.е. XORT (1) > XORT (2) > ... > XORT (NXORT), если XN > XK. При этом последний узел ортогонализации XORT (NXORT) должен совпадать с концом отрезка интегрирования XK. Если NXORT = 1, то единственный узел ортогонализации (он же конец XK) задается константой, переменной или элементом массива вещественного типа; |
K - | число условий на конце интервала XK = XORT (NXORT) (тип: целый); |
BULODN - | переменная типа Boolean, указывающая на однородность уравнения (1) и краевых условий (2) на конце отрезка XN, а именно: |
BULODN= TRUE - | когда однородны граничные условия (2), т.е. вектоp b - нулевой, и система уравнений (1), т.е. f (X) ≡ 0 на интервале интегрирования; |
BULODN=FALSE - | когда есть неоднородность в (1) либо в (2); |
B - |
двумерный вещественный массив размера
(M - K) * (M + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (2), расписанную по столбцам;
при этом вектоp b размещается в следующих
элементах массива B: B (1, M + 1), B (2, M + 1), B (3, M + 1), ... ; |
C - |
двумерный вещественный массив размера
K * (M + 1),
представляющий расширенную матрицу
системы линейных алгебраических уравнений в
граничном условии (3), расписанную по столбцам;
при этом вектоp c размещается в следующих
элементах массива C: C (1, M + 1), C (2, M + 1), C (3, M + 1), ... ; |
Y - | двумерный вещественный массив размера M * NX, в котоpом помещается вычисленное решение кpаевой задачи; |
U - | двумерный вещественный рабочий массив размера M * M; |
IRAB - | одномерный рабочий массив длины M (тип: целый); |
YR, RK - | вещественные двумерные рабочие массивы размера M * (K + 1) и M * 4, соответственно; |
T - | вещественный одномерный рабочий массив длины (K + 1) * (K + 2)/2; |
IERR - | целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом: |
IERR=65 - | если при прямой прогонке задача Коши для соответствующей однородной системы не может быть решена с точностью EPS; |
IERR=66 - | если при прямой прогонке задача Коши для системы (1) не может быть решена с точностью EPS; |
IERR=67 - | если при обратной прогонке задача Коши для системы (1) не может быть решена с точностью EPS. |
В каждом из этих случаев интегрирование системы прекращается. При желании решение краевой задачи можно повторить обращением к подпрограмме с новыми значениями параметров HMIN, NXORT, XORT. |
Версии
DE52E - | вычисление решения линейной краевой задачи для системы обыкновенных дифференциальных уравнений первого порядка методом ортогональной прогонки с расширенной (Extended) точностью. При этом параметры XN, X, HMIN, EPS, P, XORT, B, C, Y, U, YR, RK, T и параметры X, Y, Z в подпрограмме F должны иметь тип Extended. |
Для подпрограммы DE52E нестандартная подпрограмма F вычисления правой части системы должна содержать глобальную запись (структуру данных) _COM52D с элементом BUL1. Смысл параметра BUL1 в этом случае такой же, как и для DE52R. |
Вызываемые подпрограммы
DE10R - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона. |
DE10E - | вычисление решения задачи Коши для системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Mеpсона с расширенной (Extended) точностью. |
AS08R - | нахождение частного решения неоднородной системы линейных алгебраических уравнений и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
AS08E - | нахождение частного решения неоднородной системы линейных алгебраических уравнений, заданных с расширенной (Extended) точностью, и фундаментальной системы решений соответствующей однородной системы методом Жордана с выбором главного элемента по стpоке. |
AF10R - | QR - фактоpизация вещественной прямоугольной матрицы методом отражений. |
AF10E - | QR - фактоpизация вещественной прямоугольной матрицы, заданной с расширенной (Extended) точностью, методом отражений. |
Подпрограммы DE10R, AS08R, AF10R вызываются при работе подпрограммы DE52R; подпрограммы DE10E, AS08E, AF10D вызываются при работе подпрограммы DE52E. |
UTDE10 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE52R. |
UTDE11 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы DE52E. |
Замечания по использованию
Подпрограммы DE52R и DE52E предназначены для численного решения краевой задачи для системы линейных обыкновенных дифференциальных уравнений с переменными коэффициентами, имеющими непрерывные производные вплоть до 5 порядка включительно. Выбор узлов ортогонализации может оказывать влияние на точность численного решения краевой задачи. Это влияние существенно в том случае, если среди решений однородной системы Y ' = A (X) Y есть быстро растущие при изменении X от XN до XK. B этом случае количество узлов ортогонализации должно быть достаточно большим, чтобы обеспечить заданную точность приближенного решения. Если исходная система (1) является однородной, т.е. f (X) ≡ 0, то присутствие структуры данных _COM52R ( _COM52D ) с булевым элемнтом BUL1 в подпрограмме F вычисления правой части и проверка значения логической переменной BUL1 не обязательны. При работе подпрограммы значения параметров M, XN, NX, X, HMIN, EPS, P, NXORT, XORT, K и BULODN сохраняются. Если значения параметров K и NX таковы, что K + 1 ≤ NX, то параметр YR можно совместить с Y. Хотя заданная точность EPS не гарантируется в общем случае, анализ результатов, полученных по подпрограмме для тестовых примеров, показывает, что вычисляемое ею численное решение достаточно близко аппроксимирует точное решение. |
Применение программы иллюстрируется на примере дифференциального уравнения 4 порядка
(4) U '''' - 24 * U ''' - 169 * U '' - 324 * U ' - 180 * U = 0 Его частное решение U (x) = e-x - 2e-2x + e-3x удовлетворяет начальным условиям U (0) = 0, U ' (0) = 0, U '' (0) = 2, U ''' (0) = - 12. При x = 1 это решение удовлетворяет еще условиям: U (1) = 0.146996, U ' (1) = 0.0241005.
Находились численные значения U (x) как решения уравнения (4) при следующих краевых условиях:
U (0) = 0 , U ' (0) = 0 и U (1) = 0.146996 , U ' (1) = 0.0241005 .
Приводятся подпрограмма вычисления значений правой части соответствующей системы четырех дифференциальных уравнений первого порядка (к которой сводится уравнение (4) четвертого порядка) и фрагмент вызывающей программы. Так как эта система однородная, то в подпрограмме FDE52R элемент BUL1 глобальной записи отсутствует.
Unit tde52r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc, UtRes_p, FDE52R_p, DE52R_p; function tde52r: String; implementation function tde52r: String; var K,M,NXS,NXO,I,_i,IERR :Integer; XN,HMIN,EPS,P,H :Real; BUL :Boolean; IRAB :Array [0..3] of Integer; XS :Array [0..1] of Real; XO :Array [0..19] of Real; U :Array [0..15] of Real; YR :Array [0..11] of Real; RK :Array [0..15] of Real; T :Array [0..5] of Real; Y :Array [0..7] of Real; C :Array [0..9] of Real; B :Array [0..9] of Real; const C0 :Array [0..9] of Real = ( 1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.146996, 0.0241005 ); B0 :Array [0..9] of Real = ( 1.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0 ); label _1; begin for I:=0 to 9 do begin C[I] := C0[I]; B[I] := B0[I]; end; Result := ''; { результат функции } K := 2; M := 4; XN := 0.0; NXS := 2; XS[0] := 1.0; XS[1] := 0.0; HMIN := 1.E-10; EPS := 1.E-6; P := 0.1; NXO := 10; H := 0.1; XO[0] := H; for I:=2 to NXO do begin _1: XO[I-1] := XO[I-2]+H; end; BUL := True; DE52R(FDE52R,M,XN,NXS,XS,HMIN,EPS,P,NXO,XO,K,BUL,B, C,Y,U,IRAB,YR,RK,T,IERR); Result := Result + #$0D#$0A; Result := Result + ' XS' + #$0D#$0A + #$0D#$0A; for _i:=0 to 1 do begin Result := Result + Format('%20.16f ',[XS[_i]]); if ( ((_i+1) mod 2)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + ' Y' + #$0D#$0A + #$0D#$0A; for _i:=0 to 7 do begin Result := Result + Format('%20.16f ',[Y[_i]]); if ( ((_i+1) mod 4)=0 ) then Result := Result + #$0D#$0A; end; Result := Result + #$0D#$0A; Result := Result + #$0D#$0A + ' IERR=' + Format('%2d ',[IERR]); UtRes('tde52r',Result); { вывод результатов в файл tde52r.res } exit; end; end. Unit fde52r_p; interface uses SysUtils, Math, { Delphi } Lstruct, Lfunc; procedure fde52r(X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer); implementation procedure fde52r(X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer); begin Z[0] := Y[1]; Z[1] := Y[2]; Z[2] := Y[3]; Z[3] := 24.0*Y[3]+169.0*Y[2]+324.0*Y[1]+180.0*Y[0]; end; end. Результаты: Y(1, 1) = 0,1469960000000000 Y(1, 2) = -0,0000002357020366 Y(2, 1) = 0,0241005000000000 Y(2, 2) = 0,0000009105752631 Y(3, 1) = -0,2667186437063472 Y(3, 2) = 1,9999974802225444 Y(4, 1) = 0,4532534661842457 Y(4, 2) = -11,9999934834509414 IERR = 0