Вычислительная математика и программирование - Боглаев Ю.П.
ISBN 5-06-00623-9
Скачать (прямая ссылка):
с -----------------------------------------------------
DATA Т/0.6180339/
XI = В—(В —А)*Т Х2 = А+(В—А)*Т Y1 = F(X1)
Y2 = F(X2)
1 IF (Y2.LT.Y1) GO TO 2
Y = F(C)
IF (Y.EQ.O) GO TO 2
B=X2
108
V
' V
Х2=Х1
Y2=Y1
XI =В—(В—А)*Т Y1=F(X1)
GO ТО 3 А=Х1 XI =Х2
ввод: max число шагов спуска К
ШАГ СПУСКА h, НАЧАЛЬНАЯ
ТОЧКА X, точность ?
СЧЕТЧИК ЧИСЛА ШАГОВ О
ОвО
ВЫЧИСЛЕНИЕ
У-Ф(х)
""" ЗЕТ"
ВЫЧИСЛЕНИЕ grad Ф(Х)
О
ШАГ СПУСКА Х,=Х- h grad Ф(х) *
*
ВЫЧИСЛЕНИЕ У,-Ф(Х,)
|
<^(х1)<Ф(х^ ^
HETV
СТОП ПО ЧИСЛУ ИТЕРАЦИЙ
X<-Xi
У*~Уі
Рис. 3.1
109
Y1=Y2
X2=A+(B—A)*T Y2=F(X2)
3 IF ((B-A).GT.E) GO TO 1 C=(A+B)*0.5 RETURN END
Алгоритм метода градиентного спуска поиска минимума функции Ф (л'!, нескольких переменных задается формулой
(9.6.2). Дополним его процедурой дробления шага (h/2), если Ф(х(*+1))>Ф(х*). Схема алгоритма изображена на рис. 3.1.
Программа, соответствующая схеме, приведенной на рис. 3.1, имеет следующий вид:
SUBROUTINE B5S1(N,X,Y,G,H,F,F1 ,X1 ,J,K)
real; x(n),g(n),y,h,e,xi(n),yi,s
INTEGER N,K,J С ВХОДНЫЕ ПАРАМЕТРЫ
С N—РАЗМЕРНОСТЬ МАССИВОВ X,G
С X—МАССИВ, СОДЕРЖИТ ЭЛЕМЕНТЫ НАЧАЛЬНОГО
С ВЕКТОРА
С Н—НАЧАЛЬНЫЙ ШАГ СПУСКА
С Е—ТОЧНОСТЬ, ОПРЕДЕЛЯЮЩАЯ КОНЕЦ ПРОЦЕДУ-
С РЫ СПУСКА
С F—ИМЯ ВНЕШНЕЙ ПОДПРОГРАММЫ, ВЫЧИСЛЯ-
С ЮЩЕЙ ЗНАЧЕНИЕ ЦЕЛЕВОЙ ФУНКЦИИ
С SUBROUTINE F(X,Y)
С REAL X(N),Y
С F1—ИМЯ • ВНЕШНЕЙ ПОДПРОГРАММЫ, ВЫЧИС-
С ЛЯЮЩЕЙ ЗНАЧЕНИЕ ВЕКТОРА ГРАДИЕНТА ЦЕ-
С ЛЕВОЙ ФУНКЦИИ
С SUBROUTINE FI (X,G)
С REAL X(N),G(N)
С XI—РАБОЧИЙ МАССИВ
С К—МАКСИМАЛЬНОЕ ЧИСЛО ШАГОВ СПУСКА
С ВЫХОДНЫЕ ПАРАМЕТРЫ
С Y—ЗНАЧЕНИЕ ЦЕЛЕВОЙ ФУНКЦИИ
С G—ЗНАЧЕНИЕ ГРАДИЕНТА ЦЕЛЕВОЙ ФУНКЦИИ
С X—МАССИВ, СОДЕРЖИТ ЭЛЕМЕНТЫ ВЕКТОРА,
С МИНИМИЗИРУЮЩЕГО ЦЕЛЕВУЮ ФУНКЦИЮ
С J—ЧИСЛО ШАГОВ СПУСКА
С ..................................................
J=0
CALL F(X,Y)
CALL FI (X,G)
С ВЫЧИСЛЕНИЕ КВАДРАТА ЕВКЛИДОВОЙ НОРМЫ G
s=o
DO 2 1 = 1,N
110
\
ч.
2 S = S + G(I)*G(I)
С -------------------------
IF (S.LT.E) GO TO 7
J=J + 1 С ШАГ СПУСКА
3 DO 4 1=1,N
4 XI (l)=X(I)—H *G(I)
С -------------------------
CALL F(X1,Y1)
IF (Yl.LE.Y) GO TO 5 H = H*0.5 GO TO 3
5 IF (J.GT.K) GO TO 7
DO 6 1 = 1,N
6 X(I)=X1(I)
Y=Y1
GO TO 1
7 RETURN END
В качестве иллюстрации применения программы B5S1 рассмотрим минимизацию целевой функции
Ф(*1> х2) = (х1-2)4 + Х2-
Положим максимальное число шагов спуска к =100, точность е = 10-5, начальный вектор дс=(1, 2), начальный шаг А=0,5. Программа этой задачи может быть представлена в форме
REAL X(2),Y,G(2),H,E,X1 (2)
INTEGER N,K,J FXTFRNAT F FI
DATA X/1 .,2./,Е,Н/1 .E—5,0.5/,K,N/100,2/
CALL B5S1 (N,X,Y,G,H,E,F,F1,X1,J,K)
WRITE (5,1) X,Y,J 1 FORMAT (2X,2E13.6,2X,E13.6,2X,I3)
END
SUBROUTINE F(X,Y)
REAL X(2),Y
Y=(X(l)—2.)* *4+X(2)*X(2)
RETURN
END
SUBROUTINE F1(X,G)
REAL X(2),G(2)
G(l)=4.*(X(l)—2.)**3 G(2)=2. *X(2)
RETURN
END
111
3.3.13. Решение обыкновенных дифференциальных уравнений; задача Коши. Рассмотрим два алгоритма решения задачи Коши: метод Эйлера и метод Рунге—Кутта 4-го порядка. Один шаг метода Эйлера задается формулой (10.2.5). Пусть у(х)—вектор с п компонентами, заданный в точке л:; тогда один шаг интегрирования методом Эйлера—это вычисление вектора y(x+h) в точке х + h по формуле
y(x+h)=y(x)+hf{x, у(х)),
где f(x, у) — вектор-функция правых частей дифференциального уравнения (в векторной форме)
• =-*• >>•
Программа, одного шага метода Эйлера:
SUBROUTINE B6S0 (N,Y,X,H,R,D)
REAL Y (N),X,H,D (N)
INTEGER N С ВХОДНЫЕ ПАРАМЕТРЫ
С N—РАЗМЕРНОСТЬ ВЕКТОРА Y,D
С Y—МАССИВ, СОДЕРЖАЩИЙ ВЕКТОР Y(X)
С X—НАЧАЛЬНОЕ ЗНАЧЕНИЕ АРГУМЕНТА
С Н—ШАГ ИНТЕГРИРОВАНИЯ
С R—ИМЯ ВНЕШНЕЙ ПОДПРОГРАММЫ, ВЫЧИСЛЯ-
С ЮЩЕЙ ПРАВЫЕ ЧАСТИ УРАВНЕНИЙ
С SUBROUTINE R(X,Y,D)
С REAL Х,1й (N),D (N)
С D—РАБОЧИЙ МАССИВ, СОДЕРЖИТ ПРАВЫЕ ЧАСТИ
С УРАВНЕНИЙ
С ВЫХОДНЫЕ ПАРАМЕТРЫ
С X—ЗНАЧЕНИЕ АРГУМЕНТА Х+Н
С Y—МАССИВ, СОДЕРЖАЩИЙ ВЕКТОР Y(X+H)
С ...................................................
CALL R(X,Y,D)
DO 1 1 = 1,N 1 Y (I) = Y(I)+H *D(I)
X=X+H
RETURN
END
Применим программу B6S0 для интегрирования методом Эйлера с постоянным Шагом А=0,1, на интервале 0<х<5 системы дифференциальных уравнений
^±=sm(yl+y2), У1 (0)= 1,
112
= COS (у 1 Т Ху2 ), у2 (0) = 2. \
’ V
Чтобы решить эту задачу, следует выполнить 50 шагов (50 обращений к B6S0) для получения векторов
ЛМ0» !)\ /ЛЛ°>2)\ /ЛЛ5>°)\.
О/ Um/ U(5>о)л
«стартуя» с начального вектора, на каждом 5-м шаге будем выводить на терминал значения (х, >',(х)). Соответствующая программа:
REAL Y (2),X,H,D (2)
INTEGER K,N,M EXTERNAL R
DATA X/0./,Y/l.,2./,H/0.1/,K,N/50,2/
DO 1 1 = 1,К
CALL B6S0 (N,Y,X,H,R,D)
M = MOD (1,5)
IF (M.NE.O) GO TO 1 WRITE (5,2) X,Y
1 CONTINUE
2 FORMAT (2X/X=',E9.2,2X,2E13.6)
END
SUBROUTINE R(X,Y,D)
REAL X,Y (2),D (2)
D (l)=SIN (Y (1)+Y (2))
D(2)=COS (Y (l)—X* Y(2))
RETURN
END
По аналогии с программой B6S0 составим программу одного шага метода Рунге—Кутта 4-го порядка, который задается для скалярного случая формулой (10.2.12). Описание входных и выходных параметров здесь почти такое же, как в B6S0, поэтому текст программы приводится без дополнительных комментариев:
