Исследование RC-генератора синусоидальных колебаний
шена вдвое. После
вывода двух шагов признак KLP сбрасывается в ноль (24). Выполненный шаг
может быть последним на интервале (КР=1), тогда осуществляется выход из
подпрограммы(26, 27). В блоке 28 выполняется проверка, может ли быть
выполнен удвоенный шаг без выхода за пределы интервала? Если нет, то в
(29) «взводится» признак конца интервала и устанавливается величина
удвоенного шага равной оставшейся части интервала. В блоке 33 и цикле 34-35
последняя вычисленная точка делается начальной для выполнения двух малых
шагов Н и контрольного удвоенного НВ. Соответственно, в 36 устанавливается
признак двух шагов (KLP=1) и осуществляется возврат на блок 6 . Если
«дошагивание» не нужно, то в 30 проверяется, является ли точность расчетов
завышенной и в 31 можно ли удвоить малый шаг? При завышенной точности шаг
можно удвоить, если он не превзойдет максимального НМ и удвоенный шаг не
выведет за пределы интервала интегрирования. Если увеличение шага
допустимо, то блок 32 это выполняет и далее все производится как при
дошагивании, но без взвода признака конца. Если увеличение шага
недопустимо, то в цикле 37, 38 выполняется подготовка к продолжению
расчетов с прежним шагом. Из трех последних точек средняя делается
начальной для выполнения контрольного шага удвоенной величины НВ, а
последняя , - начальной для очередного малого шага Н.
6 Блок - схема алгоритма метода Рунге - Кутта с автоматическим выбором шага
[pic] Рисунок 3
7 Подпрограмма метода Рунге - Кутта с автоматическим выбором шага
Подпрограмма ARK позволяет решать произвольную систему N-го порядка
с автоматическим выбором шага интегрирования. Эта подпрограмма обращается:
к подпрограмме одного шага- SH,
к подпрограмме вычисления правых частей системы,
к подпрограмме вывода.
Подпрограмма SH записана в универсальном виде и приведена выше.
Головной вызывающий модуль, а также подпрограммы правых частей и вывода
Пользователь должен составить самостоятельно.
В главном модуле оператором EXTERNAL должны быть объявлены имена
подпрограмм правых частей и вывода. Оператором DIMENSION должны быть
объявлены массивы - фактические параметры подпрограммы ARK. Эти массивы, по
желанию, могут объявляться как одномерные или как двухмерные. Размеры
массивов (N,4),(N,3),(N,4), где N-порядок системы. Формальные имена этих
массивов в подпрограмме ARK, соответственно, X,R,F. В главном модуле первые
N элементов массива, соответствующего X, заполняются начальными условиями,
а следующие N элементов заполняются весовыми коэффициентами погрешности. В
подпрограммах правых частей и вывода в первых N элементах массива,
соответствующего X, при входе содержатся текущие значения всех N переменных
системы и не должны там переопределяться. Первые N элементов массива,
соответствующего F, должны заполняться в подпрограмме правых частей
вычисляемыми там значениями правых частей системы.
Формальными параметрами в подпрограмме правых частей должны быть
параметры (T,X,F,N), где T-независимая переменная системы.
Формальными параметрами подпрограммы вывода должны быть параметры
(T,X,F,N,IER), где IER- код ошибки, определяемый в подпрограмме ARK:
IER=0,- ошибки нет;
IER=1,- знак заданного начального шага не соответствует движению от начала
интервала интегрирования к его концу;
IER=2,- начальный шаг или/и длина интервала интегрирования ошибочно заданы
равными нулю;
IER=3,- шаг в процессе счёта стал более чем в 1000 раз меньше начального.
Массивы X и F в подпрограммах правых частей и вывода можно объявлять
как одномерные, с регулируемым размером X(N),F(N).
В главном модуле для подпрограммы ARK должны задаваться максимальный
(он же и начальный) шаг интегрирования HM, начало TN и конец TK интервала
интегрирования, а также значение требуемой абсолютной погрешности решения
E.
Подпрограмма ARK вычисляет решение системы и в каждой точке,
удовлетворяющей условиям точности, обращается к подпрограмме вывода,
передавая ей значения параметров T,X,F,IER. Пользователь может
запрограммировать здесь печать необходимых переменных или накопление их в
дополнительных массивах для последующей обработки. (В последнем случае
дополнительные массивы следует передавать в главный модуль через общую
область памяти с помощью оператора COMMON). После возврата из подпрограммы
вывода, ARK продолжает вычисление следующей точки решения.
SUBROUTINE ARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)
C Подпрограмма автоматического выбора шага.
C HM -Задаваемый максимальный шаг.
C TN,TK -Начало и конец отрезка интегрирования.
C N -Порядок системы.
C E -Задаваемое значение абсолютной погрешности.
EXTERNAL PRAV,OUT
C PRAV и OUT имена составляемых Пользователем подпрограмм правых
частей и вывода.
C IER -Выходной код ошибки.
DIMENSION X(N,4),R(N,3),F(N,4)
C Первый столбец массива X должен при входе содержать начальные
условия,
С на выходе в нем содержится решение.
C Второй столбец массива X должен при входе содержать весовые
коэффициенты погрешности.
C Первый столбец массива F должен заполняться вычисляемыми
C в подпрограмме PRAV значениями правых частей системы уравнений.
C Остальные элементы массивов X,R,F -рабочие.
DO 3 K=1,N
3 F(K,4)=X(K,2)
T=TN
HB=2*HM
IER=0
KP=0
KLP=1
CALL PRAV(T,X,F,N)
C Вызов составленной Пользователем подпрограммы правых частей
системы уравнений.
C T -Независимая переменная системы.
IF((TK-TN)*HM)4,5,60
4 IER=1
GO TO 60
5 IER=2
60 CALL OUT(T,X,F,N,IER)
C Вызов составленной Пользователем подпрограммы вывода результатов
шага
IF(IER.NE.0)RETURN
6 H=HB/2
CALL SH(T,H,X,X(1,2),F(1,2),PRAV,N,R)
8 T1=T+H
CALL SH(T1,H,X(1,2),X(1,3),F(1,3),PRAV,N,R)
T2=T+HB
CALL SH(T,HB,X,X(1,4),F,PRAV,N,R)
EH=0
DO 14 K=1,N
14 EH=EH+ABS((X(K,3)-X(K,4))/15*F(K,4))
IF(EH-E)21,21,16
16 HB=HB/2
IF(HB.LT.HM/512)IER=3
IF(IER.EQ.3)RETURN
KP=0
GO TO 6
21 T=T1
IF(KLP.EQ.1)CALL OUT(T1,X(1,2),F(1,2),N,IER)
KLP=0
CALL OUT(T2,X(1,3),F(1,3),N,IER)
IF(KP.EQ.1)RETURN
IF(ABS(TK-T2)-ABS(HB))29,29,30
29 KP=1
HB=TK-T2
GO TO 33
30 IF(EH-E/50)31,37,37
31 IF(2*H.LE.-HM.AND.ABS(TK-T2).GE.ABS(2*HB))GO TO 32
37 DO 38 K=1,N
X(K,1)=X(K,2)
38 X(K,2)=X(K,3)
GO TO 8
32 HB=2*HB
33 T=T2
DO 35 K=1,N
35 X(K,1)=X(K,3)
KLP=1
GO TO 6
END
8 Тестовая задача
Решим дифференциальное уравнение
[pic]
с начальными условиями [pic] . Легко видеть, что решением этой задачи
является функция [pic] . Вычислим решение на интервале [pic] , что
составит почти [pic] периода этой функции. Если при таком длительном
интегрировании амплитуда косинусоиды существенно не изменится, то алгоритм
численно устойчив. Можно также сравнить решение в конечной точке [pic]
Подпрограмма правых частей для этого уравнения будет такой.
SUBROUTINE PRAV(T,X,F,N)
DIMENSION X(N),F(N)
F(1)=X(2)
F(2)= -X(1)
RETURN
END
В подпрограмме вывода предусмотрим заполнение результатами массива D
для построения графиков на интервале в пять периодов, а также заполнение
массива С положительными максимумами вычисляемой функции на всем интервале
интегрирования. Эти массивы передадим в главный модуль через общую область.
SUBROUTINE OUT(T,X,F,N,IER)
DIMENSION X(N),F(N),D(3,1000),C(300)
COMMON K,L,KP,D,C
IF(T.LT.31.4)THEN
K=K+1
D(1,K)=T
D(2,K)=X(1)
D(3,K)=X(2)
ENDIF
IF(X(1).LT.0.AND.KP.EQ.1)THEN
L=L+1
C(1)=X(1)
KP=0
ENDIF
IF(X(1).GT.C(L).AND.X(1).GT.0)THEN
C(L)=X(1)
KP=1
ENDIF
IF(T.EQ.270)PRINT*,’T=270’,’ X(270)=’,X(1)
RETURN
END
В главном модуле введем исходные данные, обратимся к подпрограмме
метода, отпечатаем полученные через общую область максимумы функции и
обратимся к подпрограмме построения графика.
EXTERNAL PRAV,OUT
DIMENSION X(2,4),F(8),R(2,3),D(3,1000),C(300)
COMMON K,L,KP,D,C
READ *,N,TN,TK,HM,((X(K,J),K=1,N),J=1,2),E
K=0
L=1
C(1)=1
KP=1
CALL ARK(HM,TN,TK,X,R,F,N,E,PRAV,OUT,IER)
PRINT 1, (C(J),J=1,L)
1 FORMAT(I4/(5E15.7))
CALL KRIS(D,3,K,2,0,0.,0.)
END
9 Результаты тестирования
Графики вычисленных путем решения дифференциального уравнения функций
приведены на рисунке 4. Видно, что они близки к функциям [pic] и [pic].
[pic]
Рисунок 4
Амплитуды колебаний равны единице, период [pic] .
Выходной файл решения приведен ниже.
T=270 X(270)= 9.810482E-01
0
.1000000E+01 .9994009E+00 .9976879E+00 .9948635E+00 .9930583E+00
.9963406E+00 .9985125E+00 .9995713E+00 .9995162E+00 .9983473E+00
.9960660E+00 .9926749E+00 .9945613E+00 .9972748E+00 .9988768E+00
.9993657E+00 .9987408E+00 .9970031E+00 .9941545E+00 .9925186E+00
.9957730E+00 .9979174E+00 .9989495E+00 .9988685E+00 .9976745E+00
.9953687E+00 .9919540E+00 .9940073E+00 .9966935E+00 .9982686E+00
.9987311E+00 .9980807E+00 .9963180E+00 .9934454E+00 .9919787E+00
.9952052E+00 .9973223E+00 .9983279E+00 .
| | скачать работу |
Исследование RC-генератора синусоидальных колебаний |