Вычислительная математика и программирование - Боглаев Ю.П.
ISBN 5-06-00623-9
Скачать (прямая ссылка):
2 A(I) = Q*H/(3.*PI)
С ВЫЧИСЛЕНИЕ МАССИВА В
DO 4 1 = 1, М Q=0.
X=-PI DO 3 J=2,N,2 R1=F(X>*SIN(I*X)
R2 = F(X+H)*SIN(I*(X+H))
R3 = F(X+2.*H)*SIN(I*(X+2.*H))
3 Q=Q + R1+4.*R2-|-R3
4 B(I) = Q*H/(3.*PI)
RETURN END
100
\
’ V
Так как в программе A6S0 отсутствует контроль точности вычисления интегралов, а следовательно, и значений ат, Ьт, то необходимо будет дважды обратиться к A6S0 с различными значениями к и 2к и воспользоваться для оценки погрешности правилом Рунге, например
CALL A6S0 (F, 10,11,10, А, В CALL A6S0 (F, 10,11,20, А, В
Погрешность (массивы ЕА(1), ЕВ(1)) затем определяется в циклах типа
DO 1 1=1,Ml 1 ЕА (I) = ABS (Al (I) - A (I))/15
Алгоритмы Фурье-анализа при больших значениях т (и соответственно п), основанные на общих квадратурных формулах, например формуле Симпсона, оказываются неэффективными.
В этом случае применяются разработанные специальные быстрые алгоритмы, учитывающие специфику тригонометрических функций. Эта группа алгоритмов обычно называется быстрым преобразованием Фурье. Если обычные алгоритмы имеют порядок роста числа арифметических операций 0(т2), оо, то быстрые
алгоритмы имеют порядок роста 0(mlog2m), m->oо, что указы-
вает на их более высокую эффективность. Текст фортран-про-граммы быстрого Фурье-преобразования приводится, например, в [30, с. 166].
3.3.7. Численное дифференцирование. Пусть функция j (л;) задана таблично с постоянным шагом h
y(xi), 1 Xi+1-Xi = h.
Тогда для приближенного вычисления первой и второй производной можно использовать формулы
dx[Xi}' 2 h ’ 2^n h
d2y\_.y(xi+i)-2y(xi)+y(Xj-i) dx
Программа вычисления производных y{x) по приведенным выше формулам является простой иллюстрацией использования операторов цикла:
SUBROUTINE A7S0(N,N2,Y,H,D1,D2)
REAL Y(N),D1(N2),D2(N2),H,H1,H2 INTEGER N, N2
С -------------------------------------------------
с ВХОДНЫЕ ПАРАМЕТРЫ
101
N —РАЗМЕРНОСТЬ МАССИВА У N2—РАЗМЕРНОСТЬ МАССИВА 01,02,ЬГ1=М-1 У —МАССИВ, СОДЕРЖАЩИЙ ЗНАЧЕНИЯ ФУНКЦИИ ВЫХОДНЫЕ ПАРАМЕТРЫ
01—МАССИВ, СОДЕРЖАЩИЙ ЗНАЧЕНИЯ ПЕРВОЙ ПРОИЗВОДНОЙ
02—МАССИВ, СОДЕРЖАЩИЙ ЗНАЧЕНИЯ ВТОРОЙ ПРОИЗВОДНОЙ
Н1 =2.*Н Н2=Н*Н DO 1 1=2,N1
D1(I)=(Y(I+1)-Y(I-1))/H1 D2(l)=(Y(I+l)-2.*Y(I)+Y(I-l))/H2 RETURN END ;
3.3.8. Операции с матрицами и векторами. В этом пункте приводятся некоторые программы операций с матрицами и векторами.
1) Сложение матриц: citj=aitj+bitj, 1<г</и, 1 <п.
SUBROUTINE A8S0(A,B,M,N,C)
REAL А (М, N), В (М, N), С (М, N)
INTEGER M,N DO 1 1 = 1,М DO 1 J=1,N 1 С (I, J)=A (I, J)+В (I, J)
RETURN .
END
m
2) Умножение матриц: cik= Y aijbJik, 1 </’</, l^k^n.
k= 1 ’
SUBROUTINE A8S1 (A,B,L,M,N,C)
REAL A (L, M), В (M, N), С (L,N), S INTEGER M,N,L DO 2 1=1, L DO 2 K=1,N S=0.
DO 1 J=1,M
1 S=S+A(I,J)*B(J,K)
2 C(I,K)=S RETURN END
3) Транспонирование квадратной матрицы А с размещением транспонированной А' на месте A: a'ij=aj>i, 1<г,
SUBROUTINE A8S2(A,N)
REAL A(N,N),S
102 \
V..
с
с
с
с
с
с
с
с
с
1
INTEGER N,N1,11
N1=N —1
DO 1 1 = 1,N1
11=1 + 1
DO 1 J = I1,N
S=A(I,J)
A(I,J)=A(J,I)
1 A(J,I)=S RETURN END
4) Транспонирование прямоугольной матрицы А переходом к одномерным массивам.
SUBROUTINE A8S3 (А, М, N, MN, А1)
REAL A(MN),A1 (MN)
INTEGER M,N,MN,L1,L2 С MN—РАЗМЕРНОСТЬ ОДНОМЕРНЫХ МАССИВОВ,
С СОДЕРЖАЩИХ МАТРИЦУ А
С И ТРАНСПОНИРОВАННУЮ А1
DO 1 1=1,М DO 1 J=1,N LI =(I — l)*N+J L2=(J— 1)*M+I 1 A1(L2)=A(L1)
RETURN END
5) Скалярное произведение вектора x на вектор у
s= Z w
i= 1
REAL FUNCTION A8S4(N,X,Y)
REAL X(N),Y(N),S INTEGER N S = 0.
DO 1 1 = 1,N 1 S=S+X(I)*Y (I)
A8S4=S RETURN END
6) Вычисление нормы матрицы A
M ||= max ? \ai,j\-
1j= 1
REAL FUNCTION A8S5(N,A,S)
REAL A(N,N),S(N),R c s—РАБОЧИЙ массив
INTEGER N
103
)
DO 1 1=1, N R=0.
DO 1 J= 1, N R=R+ABS(A(I, J))
1 S(I)=R
С ВЫБОР МАКСИМАЛЬНОГО ЭЛЕМЕНТА ИЗ S(I)
R = S(1)
DO 2 1 = 2, N IF (R.LT.S(I)) R=S(I)
2 CONTINUE A8S5=,R RETURN END .1
3.3.9. Решение систем линейных алгебраических уравнений. Одним из методов рршений системы уравнений вида
x=Bx+d,
где норма матрицы В меньше 1, является метод простой итерации (см. гл. 8)
xik+1) = Bx{k)+d, к = О, 1, 2, ..., *(О) = 0.
Итерации прекращаются, когда
?l=l|x<*+i>_x<*>|KllzHM?,
где в—заданная точность приближения к решению, норма вектора
||*||= max |*{|. г!рограмма может иметь следующий вид (для 1^1
сокращения программы условие прекращения итераций принято В виде ?!<?):
SUBROUTINE A9S0 (N, В, D, Е, X, Y)
REAL B(N,N),D(N),X(N),Y(N),E,S INTEGER N С ВХОДНЫЕ ПАРАМЕТРЫ
С N—РАЗМЕРНОСТЬ МАССИВОВ B,D,X,Y
С В —МАССИВ, СОДЕРЖАЩИЙ ЭЛЕМЕНТЫ МАТРИЦЫ
С В
С D—МАССИВ, СОДЕРЖАЩИЙ ЭЛЕМЕНТЫ СВОБОД-
С НОГО ВЕКТОРА
С Е —ТОЧНОСТЬ
С Y—РАБОЧИЙ МАССИВ
С ВЫХОДНЫЕ ПАРАМЕТРЫ
С X—ПРИБЛИЖЕННОЕ РЕШЕНИЕ
с ------------------------------------------------------
С ПОЛУЧЕНИЕ ПЕРВОГО ПРИБЛИЖЕНИЯ X
DO 1 1 = 1, N 1 X(I)=D(I)
104 ^ _
V.
С УМНОЖЕНИЕ МАТРИЦЫ В НА ВЕКТОР X
100 DO 3 1=1,N
S=0.
DO 2 J= 1, N
2 S=S+B(I,J)*X(J)
3 Y (I) = S
С СЛОЖЕНИЕ BX+D—ПОЛУЧЕНИЕ НОВОГО ПРИ-
С БЛИЖЕНИЯ
DO 4 1 = 1, N
4 Y(I)=Y(I) + D(I)
С ОПРЕДЕЛЕНИЕ НОРМЫ РАЗНОСТИ Y—X
