Численные методы

5.2. Поиск нулей

Раздел содержит численные методы решения нелинейных алгебраических уравнений

Для половины методов поиск корня представляет собой итерационный процесс движения по оси X от некоторой точки x₀, выбранной в качестве нулевого приближения, к точному значению корня x^*. Новые приближения вычисляются по рекуррентной формуле x_i+1=G(x_i). Для другой половины методов задается начальный интервал поиска [a,b], сужаемый затем процедурой

Методы бывают нулевого, первого, второго порядков в зависимости от того, какой порядок имеет рекуррентная формула. Кроме того, методы различаются условиями сходимости.

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

5.2.1 Прямые методы

Уравнение должно быть выражено относительно x:

Процесс сходится, если в области поиска корня |f'(x)|<1.

Программа 72. Метод простых итераций:

10 INPUT A
20 B=A:A=...f{B}...:IF A¦B THEN 20
30 PRINT A

Размер: 28+, B

Пример: x=sin x+1/4; x₀=1,2.

В программе:

Программа 73. Метод Эйткена - Стеффенсона:

10 INPUT A
20 B=A:GOSUB 50:C=A:GOSUB 50:D=A
30 E=B-C-C+D:IF ABS E<1E^-5;PRINT D:GOTO 10
40 A=B*D/E-C*C/E:GOTO 20
50 A=...f{A}...:RETURN

Размер: 79+, E

Пример: x=sin x+1/4; x₀=1,2.

В программе:

5.2.2 Методы деления отрезка

Уравнение записывается в канонической форме:

Методы деления отрезка накладывают минимум ограничений на вид F(x): если F(x) непрерывна на отрезке поиска [a,b], и Sign F(a) ¦ Sign F(b), то существует корень x внутри [a,b]. Наилучшая сходимость наблюдается при делении отрезка пополам. Требуемое число итераций

Программа 74. Метод дихотомии.

Пока a¦b, повторять:
x=(a+b)/2; если Sign F(x)=Sign F(a), тогда a=x, иначе b=x.

10 INPUT "a",A,"b",B,"e",H:C=A:GOSUB 50:E=D:C=B:GOSUB 50:F=D
20 IF B-A<H;PRINT C:GOTO 10
30 C=(A+B)/2:GOSUB 50:G=1:IF SGN D=SGN E;G=0
40 A(G)=C:E(G)=D:GOTO 20
50 D=...f{C}...:RETURN

Размер: 120+, H

Пример: F(x)=1/(x - 1) - 1; a=-0,5; b=1,5; e=10^-5.

В программе: 50 D=1/(C-1)-2-1:RETURN

5.2.3 Метод Ньютона и его модификации

Согласно методу Ньютона, новое приближение x_i+1 определяется как пересечение оси X с прямой, которая является касательной кривой F(x) в точке x_i.

Программа 75. Метод Ньютона:

Процесс сходится, если в области поиска корня F'(x) ¦ 0 и F(x)F"(x)>0. Метод следует применять тогда, когда выражение F(x)/F'(x) вычисляется быстрее, чем дважды F(x).

10 INPUT A
20 B=A:A=B-...f{B}/f'{B}...:IF A¦B THEN 20
30 PRINT A

Размер: 30+, B

Пример: F(x)=x-sin x-1/4; (F'(x)=1-cos x); x₀=1,2.

В программе:

Программа 76. Модифицированный метод Ньютона:

В программе принято D=x₀Ч10^-6.

Модификация заключается в том, что производная берется не аналитически, а численно.

10 INPUT A:B=A*1E^-6
20 GOSUB 50:D=C:A=A+B:GOSUB 50
30 D=D*B/(C-D):A=A-D-B:IF D+A¦A THEN 20
40 PRINT A:GOTO 10
50 C=...f{A}...:RETURN

Размер: 82+, D

Пример: F(x)=x-sin x-1/4; x₀=1,2.

В программе:

Программа 77. Метод Рыбакова. Алгоритм осуществляет поиск корней на заданном промежутке [a,b]. Чтобы условия сходимости были обеспечены на всем интервале, F'(x_i) заменена некоторой константой M, которая выбирается из условия M ≥ | F'(x)|, a<x<b.

Не следует задавать M с большим запасом, так как в этом случае поиск пойдет слишком медленно.

Чтобы процесс не увязал в точках, где F(x)=0, в алгоритм вводится интервал различимости корней e (в программе задано e =10^-3).

Положить x=a.

Пока x≤b, повторять:

пока F(x)¦0, повторять: x=x+|F(x)|/M;

отпечатать x; положить x=x+e.

10 INPUT "a",A,"b",B,"M",C
20 D=...f{A}...:A=A+ABS D/C:IF A>B THEN 10
30 IF C+D=C;PRINT A:A=A+1E^-3
40 GOTO 20

Размер: 70+, D

Пример: F(x)=sin(x-0,3)+0,2x; a=-1; b=1; M=1,2.

В программе:

5.2.4 Интерполяционные методы

Интерполяционный метод порядка n заменяет функцию F(x) на интервале [a,b] интерполяционным многочленом степени n и вычисляет его корень аналитически, после чего конструирует новый многочлен.

Программа 78. Комбинированный метод секущих - хорд:

Процесс сходится, если в области поиска корня производная не меняет знак.

10 INPUT "a",A,"b",B:C=A:GOSUB 40:E=D:C=B:GOSUB 40
20 C=B-(B-A)/(D-E)*D:IF B=C;PRINT C:GOTO 10
30 A=B:B=C:E=D:GOSUB 40:GOTO 20
40 D=...f{C}...:RETURN

Размер: 97+, E

Пример: F(x)=x-sin x-1/4; x₀=1,2; x₁=1,4.

В программе:

Программа 79. Метод обратной квадратичной интерполяции:

Здесь a и b - две неподвижные опорные точки.

Процесс сходится, если F(x) непрерывна. Сходимость - быстрая.

10 INPUT "a",A,"b",B:C=A:GOSUB 60:E=D:C=B:GOSUB 60:F=D
20 C=(A+B)/2:H=(B-A)/(F-E)
30 G=C:GOSUB 60:I=(C-A)/(D-E)
40 C=A-E*(I+D*(I-H)/(F-D)):IF C¦G THEN 30
50 PRINT C:GOTO 10
60 D=...f{C}...:RETURN

Размер: 140+, I

Пример: F(x)=x-sin x-1/4; x₀=1,2; x₁=3,4.

В программе:

5.2.5 Методы решения систем нелинейных уравнений

Представлены два популярных метода решения систем - метод Зейделя и метод Ньютона - Канторовича. Обе программы снабжены минимумом сервиса. Для контроля сходимости процесса приостанавливайте счет и просматривайте промежуточные результаты в калькуляторном режиме.

Программа 80. Метод Зейделя. Система записана в нормальной форме:

Условие сходимости:

Рекуррентная формула:

Условие окончания счета:

10 A=n:FOR B=1 TO A:PRINT B;:INPUT D(B):NEXT B
20 INPUT "e",C
30 $=CHR 1:FOR B=1 TO A:D=D(B):GOSUB 60+10*B
40 IF ABS (D(B)-D)>C;$=""
50 NEXT B:IF $="" THEN 30
60 FOR B=1 TO A:PRINT B;D(B):NEXT B:GOTO 20
70 E=...f₁{E,F,...,D(A)}...:RETURN
80 F=...f₂{E,F,...,D(A)}...:RETURN
90 ...
... D(A)=...f_n{E,F,...,D(A)}...:RETURN

Размер: 131+, D(n)

Пример:

В программе:

Программа 81. Метод Ньютона - Канторовича. Система записана в канонической форме:

Линеаризованная система:

где W(x) - матрица частных производных векторной функции F(x⁽ⁱ⁾) в точке x⁽ⁱ⁾.

Рекуррентная формула:

На каждой итерации методом вращений решается линеаризованная система, затем по рекуррентной формуле вычисляется новое приближение.

Условие окончания счета:

A=n
I=A*A+A
DEFM A+I-17
10 FOR B=1 TO A:PRINT B;:INPUT I(I+B):NEXT B
20 INPUT "e",H
30 FOR B=1 TO A:I(I+B)=I(I+B)+H:GOSUB 190:FOR C=1 TO A
40 I(A*C+B)=I(C):NEXT C:I(I+B)=I(I+B)-H:NEXT B:GOSUB 190
50 FOR B=1 TO A:FOR C=1 TO A:I(A*B+C)=(I(A*B+C)-I(B))/H:NEXT C
60 NEXT B:FOR B=1 TO A-1:FOR C=B+1 TO A
70 F=I(B*A+B):G=I(C*A+B):IF F=0;IF G=0 THEN 120
80 E=SQR (F*F+G*G):D=G/E:E=F/E
90 FOR F=1 TO A:G=E*I(B*A+F)+D*I(C*A+F)
100 I(C*A+F)=E*I(C*A+F)-D*I(B*A+F):I(B*A+F)=G:NEXT F
110 G=E*I(B)+D*I(C):I(C)=E*I(C)-D*I(B):I(B)=G
120 NEXT C:NEXT B
130 $=CHR 1:FOR B=A TO 1 STEP -1:D=0:IF B=A THEN 150
140 FOR C=B+1 TO A:D=D+I(I-A+C)*I(B*A+C):NEXT C
150 I(I-A+B)=(I(B)-D)/I(B*A+B):I(I+B)=I(I+B)-I(I-A+B)
160 IF ABS I(I-A+B)>H;$=""
170 NEXT B:IF $="" THEN 30
180 FOR B=1 TO A:PRINT B;I(I+B):NEXT B:GOTO 20
190 J=...f₁{J(I),K(I),...,I(I+A)}...
200 K=...f₂{J(I),K(I),...,I(I+A)}...
... I(A)=...f_n{J(I),K(I),...,I(I+A)}...:RETURN

Размер: 566+, I(n*(n+2))

Пример:

В программе:

5.3. Поиск экстремальных точек

Поиск экстремальных точек функции - операция, применяемая при оптимизации решений. Соответствующие прикладные программы получили чрезвычайно широкое распространение. Ресурсы "МК 85" не позволяют реализовывать оптимизационные методы в приемлемом для использования виде. Исключение составляют лишь некоторые алгоритмы непрерывной безусловной оптимизации.

Одномерная задача представлена двумя высокоэффективными методами, первый из которых основан на процедуре деления отрезка, второй - на замене функции интерполяционным многочленом.

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

5.3.1 Одномерная задача

Программа 82. Метод золотого сечения. Процесс сходится к локальному минимуму. Для определения максимума следует изменить знак алгебраического отношения в строке 30. В качестве исходных данных указывается интервал поиска [a,b]. Деление отрезка в пропорции золотого сечения позволяет обойтись минимальным количеством вычислений F(x). Требуемое число итераций:

Положить: x_b=a+d(b-a); x_a=b-d(b-a);
Пока x_b-x_a> e, повторять:
если F(x_b)>F(x_a), тогда b=x_b; x_b=x_a; x_a=b-d(b-a);
иначе a=x_a; x_a=x_b; x_b=a+d(b-a).

10 INPUT "a",A,"b",B:G=SQR 5/2-.5:H=0:GOSUB 50:H=1:GOSUB 50
20 IF C=D;PRINT "x";C,"F";E:GOTO 10
30 H=0:IF E>F;H=1
40 A(H)=D(-H):D(-H)=C(H):F(-H)=E(H):GOSUB 50:GOTO 20
50 C(H)=A(H)+G*(B(-H)-A(H)):E(H)=...f{C(H)}...:RETURN

Размер: 162+, H

Пример: F(x)=|x-1|; a=0; b=2.

В программе: 50 ... :E(H)=ABS (C(H)-1):RETURN

Программа 83. Метод квадратичной интерполяции:

где s - шаг опорных точек.

Если F(x) - гладкая, процесс сходится; при F"(x₀)>0 - к минимуму, при F"(x₀)<0 - к максимуму.

10 INPUT "s",A,"x",B
20 E=0:C=B:GOSUB 50:E=C:C=B-A:GOSUB 50:D=C:C=B+A:GOSUB 50
30 C=(D-C)/(D+C)/2:B=B+C*A:IF ABS C>1E^-9 THEN 20
40 PRINT "x";B,"F";E:GOTO 10
50 C=...f{C}-E...:RETURN

Размер: 123+, E

Пример: F(x)=exp(-x²); s=1; x₀=1/2.

В программе: 50 C=EXP(-C*C)-E:RETURN

5.3.2 Многомерная задача

Программа 84. Метод покоординатного спуска. Требует задания начального приближения r₀, начального шага поиска s и допустимой линейной погрешности e. По каждой из n координат проводится грубый одномерный поиск минимума методом поразрядного приближения с шагом s и s=-s/3, затем шаг уменьшается в той же пропорции и процедура повторяется. Условие завершения - s<e.

10 A=n:FOR D=1 TO A:PRINT D;:INPUT G(D):NEXT D:GOSUB 80
20 INPUT "s",C,"e",E
30 FOR D=1 TO A:B=C
40 G(D)=G(D)+C:G=F:GOSUB 80:IF F<G THEN 40
50 IF B=C;C=-C/3:GOTO 40
60 NEXT D:IF ABS C>E THEN 30
70 FOR D=1 TO A:PRINT D;G(D):NEXT D:PRINT "F";F:GOTO 10
80 F=...f{H,I,...,G(A)}...:RETURN

Размер: 159+, G(n)

Пример: F(x,y,z)= e^x+y+z/x/y²/z³
x₀=0,5; y₀=0,5; z₀=0,5; s=0,5; e= 10^-4.

В программе:

Программа 85. Метод наискорейшего спуска. Требует задания начального приближения r₀, начального шага поиска s и допустимой линейной погрешности e. В текущей точке r_i численно определяется направление, противоположное местному градиенту минимизируемой функции. Вдоль этого направления проводится точный одномерный поиск экстремума поразрядным приближением с кратностью шага 3. Найденная точка принимается за новое приближение r_i+1, для которого процедура повторяется с начальным шагом s=|r_i+1 - r_i|. Условие окончания счета - s<e.

10 A=n:FOR B=1 TO A:PRINT B;:INPUT H(B):NEXT B
20 INPUT "s",G,"e",E
30 IF G<E;FOR B=1 TO A:PRINT B;H(B):NEXT B:PRINT "F";F:GOTO 20
40 C=G:H=0:GOSUB 90:D=F:FOR B=1 TO A:H(B)=H(B)-E:GOSUB 90
50 H(B+A)=(F-D)/E:H(B)=H(B)+E:H=H+H(B+A)-2:NEXT B:H=SQR H:G=0
60 FOR B=1 TO A:H(B)=H(B)+C/H*H(B+A):NEXT B:G=G+C
70 GOSUB 90:IF F≥D;C=-C/3:IF ABS C<E/5 THEN 30
80 D=F:GOTO 60
90 F=...f{I,J,...,H(A)}...:RETURN

Пример: 258+, H(2*n)

Пример: F(x,y)=x²+y²; x₀=1; y₀=1; s=1,5; e =10^-5.

В программе:

5.4. Численное интегрирование

Методы численного интегрирования имеют дело с определенным интегралом вида:

Несобственные интегралы заменяются определенным интегралом путем выбора соответствующих пределов a, b и деления интервала интегрирования на части. Возможны следующие варианты:

1. Если функция содержит разрыв второго рода (f(x)=г), то выделяется e-окрестность разрыва, в которой интеграл вычисляется многократно, с поэтапным уменьшением шага интегрирования, до тех пор, пока результат не сойдется (функция должна быть интегрируемой). Затем с обычным шагом вычисляются интегралы вокруг окрестности разрыва. Суммарный интеграл:

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

где s - шаг интегрирования, не сойдется. f(x) должна быть абсолютно интегрируемой - быстро затухать при стремлении x к бесконечности. Шаг интегрирования увеличивают по ходу вычислений, используя правило: если (1-I_i-1/I_i)<e, тогда s=s+s. Здесь e - порог, задаваемый в зависимости от требуемой точности.

5.4.1 Квадратурные формулы Ньютона - Котеса

Интегрирование ведется по равноотстоящим отсчетам. Представлены формулы 2-го и 6-го порядков. Интеграторы других порядков могут быть получены путем незначительной модификации программы 87.

Интегрируемая кривая записывается в тексте программы. В процессе счета задаются пределы интегрирования [a,b] и число интервалов разбиения n, влияющее на точность.

Программа 86. Формула Симпсона.

Формула точна для многочлена третьей степени.

10 INPUT "a",A,"b",B,"n",C:C=(B-A)/C/2:D=0:GOSUB 40:E=F
20 IF A≤B;GOSUB 40:D=2-D:E=E+(2+D)*F:GOTO 20
30 PRINT (E-F)*C/3:GOTO 10
40 F=...f{A}...:A=A+C:RETURN

Размер: 109+, F

Пример: f(x)=x³; a=0; b=1; n=2.

В программе: 40 F=A-3: ...

Программа 87. Формула Ньютона - Котеса.

Формула точна для многочлена седьмой степени и применяется к громоздким зависимостям.

10 INPUT "a",A,"b",B,"n",C:G=272:H=27:I=216:J=82
20 C=(B-A)/C/6:GOSUB 50:E=41*F
30 IF A>B;PRINT (E-41*F)*C/140:GOTO 10
40 FOR D=-2 TO 3:GOSUB 50:E=E+G(ABS D)*F:NEXT D:GOTO 30
50 F=...f{A}...:A=A+C:RETURN

Размер: 142+, J

Пример: f(x)=x⁷; a=0; b=2; n=2.

В программе: 50 F=A-7: ...

5.4.2 Квадратурные формулы Чебышева и Гаусса

Интегрирование ведется по неравноотстоящим отсчетам. Общая формула:

где k - число отсчетов и порядок интегратора;
A_j - удельный вес отсчета;
f(t_j) - отсчет интегрируемой функции f(t) в точке -1 ≤ t_j ≤ 1.

В формулах Чебышева все A_j=2/k.

Программа 88. Формула Чебышева 3-го порядка.

10 G=SQR 2/4:INPUT "a",A,"b",B,"n",C:C=(B-A)/C:F=0
20 FOR A=A+C/2 TO B STEP C:FOR D=-1 TO 1:E=A+SGN D*F(ABS D)*C
30 E=...f{E}...
40 F=F+E:NEXT D:NEXT A:PRINT F*C/3:GOTO 10

Размер: 105+, G

Пример: f(x)=x³; a=0; b=1; n=2.

В программе: 30 F=F+E-3:NEXT D:NEXT A:PRINT F*C/3:GOTO 10

Программа 89. Формула Чебышева 4-го порядка.

10 INPUT "a",A,"b",B,"n",C:C=(B-A)/C:F=0
20 FOR A=A+C/2 TO B STEP C:D=C
30 E=.0937962:GOSUB 50:E=.3972724:GOSUB 50:IF D=C;D=-D:GOTO 30
40 NEXT A:PRINT F*C/4:GOTO 10
50 E=A+D*E:F=F+...f{E}...:RETURN

Размер: 134+, H

Пример: f(x)=x⁴; a=0; b=1; n=2.

В программе: 50 E=A+D*E:F=F+E-4:RETURN

Программа 90. Формула Гаусса 3-го порядка.

10 G=SQR .15:H=8:I=5:INPUT "a",A,"b",B,"n",C
20 C=(B-A)/C:F=0:FOR A=A+C/2 TO B STEP C
30 FOR D=-1 TO 1:E=A+SGN D*F(ABS D)*C
40 E=...f{E}...
50 F=F+H(ABS D)*E:NEXT D:NEXT A:PRINT F*C/18:GOTO 10

Размер: 128+, I

Пример: f(x)=x⁵; a=0; b=1; n=2.

В программе: 40 E=E-5: ...

Программа 91. Сложная квадратура Гаусса для двумерного интеграла. Вычисляется интеграл по прямоугольной области:

Интервал [a₁,b₁] делится на n₁ отрезков, [a₂,b₂] - на n₂ отрезков. В результате получается n₁n₂ элементарных прямоугольников, к каждому из которых применяется квадратура, включающая четыре точки с координатами (¦Ц(1/3), ¦Ц(1/3)) и равными весами.

10 INPUT "a1",A,"b1",B,"n1",C,"a2",D,"b2",E,"n2",F:G=.5/SQR 3
20 H=0:C=(B-A)/C:F=(E-D)/F:FOR A=A+C/2 TO B STEP C
30 FOR L=D+F/2 TO E STEP F:FOR I=-1.5 TO 1.5
40 J=A+SGN I*G*C:K=L+SGN (ABS I-1)*G*F
50 J=...f{J,K}...
60 H=H+J:NEXT I:NEXT L:NEXT A:PRINT H*C*F/4:GOTO 10

Размер: 186+, L

Пример: f(x)=x₁+x₂-x₁x₂; a₁=0; b₁=2; n₁=4; a₂=0; b₂=3; n₂=5.

В программе: