31. Некоторые вычислительные алгоритмы

В этом разделе приводятся некоторые алгоритмы приближенных вычисле­ний: метод простой итерации и метод Ньютона решения алгебраических уравне­ний; метод Гаусса и итерационные методы решения линейных систем; метод наименьших квадратов для аппроксимации таблично заданной функции; квадра­турные формулы трапеций и Симпсона для приближенного вычисления опреде­ленных интегралов; метод Эйлера и метод Рунге-Кутта численного решения за­дачи Коши. Мы опустим обоснование этих численных методов, а приведем лишь общие схемы алгоритмов и примеры программ.

Приближенное решение алгебраических уравнений

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

Метод простой итерации определен для уравнения вида х=Ф(х), любое уравнение F(x)=0, очевидно, можно преобразовать к такому виду (однако не вся­кое такое преобразование будет подходящим для метода простой итерации). Ал­горитм метода простой итерации крайне легок и заключается в следующем: вы­бирается некоторое начальное приближение X0, затем вычисляются X1 X2 и так далее по формуле

X+1 = Ф (X) ; j=0,1,... ^ до тех пор, пока не выполнится условие / X+1- X / < е , где е - некоторое наперед заданное маленькое число, которое называется погрешностью, или точностью. Во многих случаях удается найти корень с абсолютной компьютерной точно­стью, т.е. для некоторого конечного j выполняется условие X+1= X ; при этом не нужно задавать никакой точности, однако автор не гарантирует, что это справед­ливо для всех уравнений. Метод простой итерации сходится, т.е. за конечное число итераций дает значение корня с заданной точностью, при удачном выборе вида функции Ф(х) и нулевого приближения. Запишем программу, реализующую метод простой итерации:

TYPE FuncType = FUNCTION (x:Real):Real; FUNCTION Fi(x:Real):Real; FAR;

{ Решаем уравнение xA5 - xA2 - ln(2+xA2 )=0, для метода простой итерации

преобразуем его к виду x = (xA2+ln(2+xA2))A0.2 }

BEGIN Fi:=Exp(0.2*Ln(Sqr(x)+Ln(2+Sqr(x)))); END;

FUNCTION SimpleIteration(x0:Real; f:FuncType):Real;

{ Ищем корень с абсолютной точностью }

VAR x1 : Real;

BEGIN x1:=x0;

REPEAT x0:=x1; x1:=f(x0); UNTIL x1=x0; SimpleIteration:=x1;

END;

CONST X0 = 1.0;

VAR Root : Real;

BEGIN Root:=SimpleIteration(X0,Fi);

WRITELN('Корень уравнения x=Ф(x) равен ',Root); WRm^N^ этой точке x-Ф^) равно ',Root-Fi(Root));

Метод Ньютона применяется для уравнения F(x)=0 и записывается в виде:

XJ+1 = Xj - 1,2,... F'(XJ)

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

TYPE FuncType = FUNCTION (x:Real):Real;

FUNCTION F(x:Real):Real; FAR;

{ Решаем уравнение       -       - ln(2+xA2 )=0 }

VAR x2 : Real;

BEGIN x2:=Sqr(x); F:=x*Sqr(x2)-x2-Ln(2+x2); END; FUNCTION dF(x:Real):Real; FAR;{ Производная функции F }

VAR x2 : Real;

BEGIN x2:=Sqr(x); dF:=5*Sqr(x2)-2*x*(1+1/(2+x2)); END; FUNCTION Newton(x0:Real; f,df:FuncType):Real; { Ищем корень с абсолютной точностью }

VAR x1 : Real;

BEGIN x1:=x0;

REPEAT x0:=x1; x1:=x0-f(x0)/df(x0); UNTIL x1=x0;

Newton:=x1; END;

CONST X0 = 1.0;

VAR Root : Real;

BEGIN Root:=Newton(X0,F,dF);

WRITELN('Корень уравнения F(x)=0 равен ',Root); WRITELN('В этой точке F(x) равна ',F(Root));

END.

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

Следующая группа вычислительных алгоритмов применяется для решения систем линейных алгебраических уравнений A*X=B. Для не слишком больших матриц применяют точный метод исключения Гаусса с выбором главного эле­мента. Система решается в два этапа: прямой ход приводит матрицу A к тре­угольному виду:

aij = aij ~ akj ' dik;bi = bi ~ bk ' dik;dik =

i = k + 1,k + 2,K ,n;j = k,k + 1,K ,n;k = 1,K ,n.

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

=   max   {alk}. i = kK n

Обратный ход дает искомый вектор X :

( n

bn 1

i = n - 1,n - 2,K ,1.

Несколько модифицированный метод Гаусса можно с успехом применять для вычисления определителя и обратной матрицы. Запишем программу, реализую­щую этот метод: CONST Nmax=40;

{Максимальный размер матрицы, вы можете выбрать другое число} TYPE VectorType = ARRAY[1..Nmax] OF Real;

MatrixType = ARRAY[1..Nmax] OF VectorType; FUNCTION Gauss( n : Byte {размер системы}; A : MatrixType {матрица системы};

B : VectorType {правые части}; VAR X : VectorType {вектор неизвестных}) : Boolean;

{ Функция будет возвращать TRUE, если систему удалось решить, и FALSE, если мат­рица системы вырождена }

VAR i,j,k,iMax : Byte; tmp,Max,d : Real; v : VectorType; BEGIN FOR k:=1 TO n-1 DO BEGIN { прямой ход } { ищем главный элемент }

Max:=Abs(A[k,k]); iMax:=k; FOR i:=k+1 TO n DO

IF Abs(A[i,k])>Max THEN BEGIN Max:=Abs(A[i,k]); iMax:=i; END; IF Max=0 THEN BEGIN { матрица вырождена } Gauss:=FALSE; Exit;

END;

IF iMax<>k THEN BEGIN { переставляем строки }

Tmp:=B[k]; B[k]:=B[iMax]; B[iMax]:=Tmp; v:=A[k]; A[k]:=A[iMax]; A[iMax]:=v; END;

FOR i:=k+1 TO n DO BEGIN {вычитаем из i-ой строки k-ю } d:=A[i,k]/A[k,k];

FOR j:=k TO n DO A[i,j]:=A[i,j]-d*A[k,j]; B[i]:=B[i]-d*B[k]; END;

END;

{ сейчас матрица системы - треугольная }

IF A[n,n]=0 THEN BEGIN { матрица вырождена } Gauss:=FALSE; Exit; END;

{ обратный ход }

X[n]:=B[n]/A[n,n];

FOR i:=n-1 DOWNTO 1 DO BEGIN

tmp:=B[i];

FOR j:=i+1 TO n DO tmp:=tmp-A[i,j]*X[j];

X[i]:=tmp/A[i,i]; END;

Gauss:=TRUE;

END;

VAR n,i,j : Byte; a : MatrixType; b,x : VectorType; BEGIN WRITE('Введите размер системы '); READ(n);

WRITELN('BBeflHTe расширенную матрицу системы'); FOR i:=1 TO n DO BEGIN

FOR j:=1 TO n DO READ(a[i,j]); READ(b[i]); END; IF NOT Gauss(n,a,b,x) THEN BEGIN

WRITELN('Матрица системы вырождена'); Halt; END; WRITELN('Рeшeниe системы и невязки :'); FOR i:=1 TO n DO BEGIN

FOR j:=1 TO n DO b[i]:=b[i]-a[i,j]*x[j];

WRITELN(x[i]:12,' ',b[i]:12); END;

END.

Хотя метод Гаусса и называют точным методом, невязки или погрешности, которые мы вывели в программе, не будут равны нулю. Но эти погрешности обу­словлены не неточностью самого метода, а исключительно погрешностью вы­числения арифметических операций.

Для очень больших систем, когда метод Гаусса становится неэффективным, применяют итерационные методы, например, метод простой итерации или ме­тод Зейделя. Вычисления по методу простой итерации начинаются с произ­вольного вектора X ={xi , x2 xn }. Итерационный процесс осуществляется по формуле:

■1        ■      1    ( n Л

j+1       j ,

X;        = X; Л--

i i

a

bi - X     • xJk ;i = L2,K n;j = 0,1,K

т.е. все неизвестные на следующей итерации вычисляются только через неиз­вестные на предыдущей итерации.

В методе Зейделя используется итерационная формула

1     ( i-1 n N

j+1        j , 1

x-    = x Л—

bi - Xaik • xk+1 " X• xk ;i = 1,2,K n;j = 0,1,K ,

aii   V        k=1 k=i J

в которой при вычислении очередного неизвестного используются последние найденные значения остальных неизвестных. Вычисления заканчиваются, когда невязки системы становятся достаточно малыми. Итерационные методы сходятся не для всякой матрицы. Достаточным условием сходимости является положи­тельная определенность матриц. Запишем программу, использующую метод простой итерации и метод Зейделя: CONST Nmax=6;

TYPE VectorType = ARRAY[1..Nmax] OF Real;

MatrixType = ARRAY[1..Nmax] OF VectorType; PROCEDURE SimpleIteration(n:Byte; A:MatrixType;B:VectorType;Epsilon:Real; VAR X:VectorType);

VAR i,k : Byte; X0 : VectorType; s,Max,Tmp : Real; BEGIN X0:=X;

REPEAT

FOR i:=1 TO n DO BEGIN s:=B[i];

FOR k:=1 TO n DO s:=s-A[i,k]*X0[k]; X[i]:=X0[i]+s/A[i,i]; END; { Вычисляем максимальную невязку }

Max:=0;

FOR i:=1 TO n DO BEGIN Tmp:=Abs(X[i]-X0[i]);

IF Max<Tmp THEN Max:=Tmp; END;

X0:=X;

UNTIL Max<Epsilon;

END;

PROCEDURE Zeidel (n:Byte; A:MatrixType; B:VectorType; Epsilon:Real; VAR X:VectorType);

VAR i,k : Byte; s,Max,Tmp : Real;

BEGIN REPEAT Max:=0;

FOR i:=1 TO n DO BEGIN s:=B[i];

FOR k:=1 TO n DO s:=s-A[i,k]*X[k]; Tmp:=X[i]; X[i]:=X[i]+s/A[i,i];

Tmp:=Abs(X[i]-Tmp); IF Max<Tmp THEN Max:=Tmp; END; UNTIL Max<Epsilon;

END;

CONST a : MatrixType = ( (110, 3, 4, -5, 16, 20), ( 1, В5, -6, 12, 11, -16), ( 24, -3,

54, -2, 0, 11),

( -В, 9, 0, 75, -9, 12), ( 0, 2, 3, -5, 9В, 34),   ( 16, 9, -

4, -11, 0, 5));

b : VectorType = (100,200,300,400,500,600); x0 : VectorType = (0,0,0,0,0,0);

CONST Eps = 1E-8;

VAR n,i,j : Byte; x : VectorType; d : Real;

BEGIN

n:=Nmax; x:=x0; SimpleIteration(n,a,b,Eps,x);

WRITELN('Рeшeниe системы методом простой итерации и невязки:'); FOR i:=1 TO n DO BEGIN d:=b[i];

FOR j:=1 TO n DO d:=d-a[i,j]*x[j];

WRITELN(x[i]:12,' ',d:12); END;

x:=x0;

Zeidel(n,a,b,Eps,x);

WRITELN('Рeшeниe системы методом Зейделя и невязки :'); FOR i:=1 TO n DO BEGIN

d:=b[i]; FOR j:=1 TO n DO d:=d-a[i,j]*x[j];

WRITELN(x[i]:12,' ',d:12); END;

END.

Аппроксимация таблично заданной функции методом наименьших квадратов

Этот алгоритм часто используется в прикладных задачах, когда необходимо представить некоторый массив экспериментальных или наблюдательных данных какой-нибудь функциональной зависимостью. При этом вид аппроксимирующей функции задается заранее, а ее числовые коэффициенты подбираются таким об­разом, чтобы функция наилучшим образом аппроксимировала имеющиеся дан­ные. Мы изучим способ аппроксимации алгебраическим полиномом. Пусть зада­ны массивы значений аргумента (xjj j=1,...,m и функции {yj} j=1,...,m и требуется аппроксимировать данную функцию алгебраическим полиномом степени n :

P(x) = a0 + a1 • x + a2 • x2 +K an • xn.

Коэффициенты полинома ищутся из условия минимума функционала

А = X (/'(х ) - У]) .

j=1

Каждая из величин P(xj)-yj называется невязкой в точке Xj, а величина А, очевид­но, есть сумма квадратов невязок. Таким образом, нашей задачей является выбор таких коэффициентов a0,a1,...,an, при которых сумма квадратов невязок достигает наименьшего значения - отсюда и название метода. Очевидно, что функционал достигает минимума при выполнении условий:

|А к | {р(^) _ у.) = 0;

|А . |xj( p(xj) _ у,) = 0;К; |А «|xn (p(xj) - yj) =°.

Введем для удобства записи обозначения

mm

Xk = X х] ; Yk = X x) ■ yj и перепишем систему уравнений в виде:

X1 ' a0 + X2 ' a1 + X3 ' a2 +K + Xn+1 ' an = Y1

KKKKKKKKKKKKKKKKKK

Это система линейных алгебраических уравнений относительно коэффициентов 0/; ее решение можно найти одним из известных методов.

{ Аппроксимация методом наименьших квадратов } CONST m=10; n=2;

TYPE VectorType = ARRAY[1..n+1] OF Real;

MatrixType = ARRAY[1..n+1] OF VectorType;

DataType   = ARRAY[1..m] OF Real;

XType       = ARRAY[0..2*n] OF Real;

YType       = ARRAY[0..n]  OF Real; { Используем функцию Gauss из предыдущего пункта } CONST x : DataType = (1,2, 3, 4, 5, 6, 7, 8, 9, 10);

y : DataType = (4,9,16,25,36,49,64,81,100,121); VAR i,j,k : Byte; XX : Xtype; YY : Ytype; A : MatrixType; B,Coeff : VectorType; p :

Real;

BEGIN { Вычислим величины X и Y } FOR k:=0 TO 2*n DO XX[k]:=0; FOR i:=1 TO m DO BEGIN p:=1; k:=0;

WHILE k<=2*n DO BEGIN XX[k]:=XX[k]+p; p:=p*x[i]; Inc(k); END;

END;

FOR k:=0 TO n DO YY[k]:=0;

FOR i:=1 TO m DO BEGIN p:=y[i]; k:=0;

WHILE k<=n DO BEGIN YY[k]:=YY[k]+p; p:=p*x[i]; Inc(k); END;

END;

{ сформируем матрицу системы }

FOR i:=1 TO n+1 DO FOR j:=1 TO n+1 DO A[i,j]:=XX[i+j-2]; FOR i:=1 TO n+1 DO B[i]:=YY[i-1]; { вычислим коэффициенты полинома }

IF Gauss(n+1,A,B,Coeff) THEN; { выведем результаты }

WRITELN('Коэффициенты полинома :');

FOR k:=1 TO n+1 DO WRITELN('a[',k-1,']=',Coeff[k]);

END.

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

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

и

j"f(x)dx

2f(a) + f(a + h) + f(a + 2h)+K + 2f(b)

где Н = (Ь-а)/п ; п - некоторое целое положительное число, называемое числом разбиений отрезка [а,Ь] .

Квадратурная формула Симпсона записывается в виде:

\ f(x)dx h ■ [f(a)

4f(a + h) + 2f(a + 2h) + 4f(a + 3h)+K +4f(b - h) + f(b)].

Число разбиений п для формулы Симпсона должно быть четным. Относи­тельная погрешность вычисленного значения оценивается по формуле:

5 = -г

'2n

J

2n

Здесь 1п и 12п - значения, вычисленные для числа разбиений п и 2п соответствен­но, г=3 для квадратурной формулы трапеций и г=12 для квадратурной формулы Гаусса. Чтобы начать вычисления, берут п=1 или п=2 и удваивают его до тех пор, пока относительная погрешность не станет меньше некоторого наперед за­данного маленького числа.

a

Непосредственное использование квадратурных формул крайне неэффектив­но, поскольку приходится многократно повторять выглисление подынтегральной функции в тех точках, где она уже посчитана. Наиболее эффективный алгоритм выглислений по формуле трапеций можно записать в виде:

jf{x}dx = J e (b - a);    = ; jn =       + 1 £f;  fk = f(a + k(b - a)).

a 2 i=l n

Для квадратурной формулы Симпсона наиболее эффективный алгоритм можно записать в виде:

jf(x)dx = Je (b - a); Jn = ^пзП^П; ф1 = fo +    фn = фn +     Sn = 2^fai-i .

a 2 j=1 n

Квадратурная формула Симпсона точнее, чем формула трапеций, и в боль­шинстве случаев выгчисление интеграла по формуле Симпсона происходит на­много быстрее.

TYPE FuncType = FUNCTION (x:Real):Real;

FUNCTION Trapezia(a,b:Real; f:FuncType; Epsilon:Real):Real;

{ Квадратурная формула трапеций } VAR i,n : LongInt; J1,J2,d,h,s : Real; BEGIN J1:=(f(a)+f(b))/2; n:=1;

REPEAT n:=n ShL 1; h:=(b-a)/n; s:=0;

FOR i:=1 TO n DIV 2 DO s:=s+f(a+(2*i-1)*h); J2:=J1/2+s/n; d:=Abs(J1/J2-1)/3; J1:=J2; UNTIL d<Epsilon; Trapezia:=J2*(b-a);

END;

FUNCTION Simpson(a,b:Real; f:FuncType; Epsilon:Real):Real;

{ Квадратурная формула Симпсона } VAR i,n : LongInt; J2,Jn,Q2,Qn,d,h,s : Real;

BEGIN Q2:=f(a)+2*f((a+b)/2)+f(b); J2:=(Q2+2*f((a+b)/2))/6; n:=2; REPEAT n:=n ShL 1; h:=(b-a)/n; s:=0;

FOR i:=1 TO n DIV 2 DO s:=s+f(a+(2*i-1)*h); s:=s*2; Qn:=Q2+s; Jn:=(Qn+s)/3/n; d:=Abs(J2/Jn-1)/12; J2:=Jn; Q2:=Qn;

UNTIL d<Epsilon;

Simpson:=Jn*(b-a);

END;

FUNCTION F(x:Real):Real; FAR; BEGIN F:=Sqr(Sqr(x)); END; BEGIN WRITELN(Trapezia(1,2,F,1E-6)); WRITELN(Simpson(1,2,F, 1E-6));

END.

Численное решение задачи Коши

Задача Коши для численного решения ставится следующим образом: найти решение дифференциального уравнения y'=f(x,y) с краевым условием y(x0)=y0 в точке x=X (или на некотором конечном множестве точек). Будем строить свойвычислительный алгоритм так: разобьем отрезок [х°,Х] на п равных частей точ­ками х0,х1,Х2,...,хп, при этом, очевидно, х0=х0 и хп=Х. Обозначим к=(Х-х°)/п. Зная х0 и у(х0) и пользуясь одним из известных численных методов, найдем у(х1); зная х1 и у(х1), найдем у(х2) и так далее, пока не будет найден у(хп), равный (с некото­рой погрешностью) искомой величине У. Обозначим это найденное значение Уп -для п разбиений отрезка. Теперь увеличим число разбиений вдвое так же, как при численном интегрировании, и найдем Т2п. Оценив относительную погреш-

ность 5 =

|, либо закончим вычисления, либо еще раз удвоим n. Суще­ствует много численных методов, позволяющих по известным значениям x и y(x) найтиy(x+h). Самый простой из них называется методом Эйлера:

y(x + h) = y(x) + f(x + y2,y + k)2)-h ;   k1 = f(x,y) • h.

Более сложен (но и намного более точен) метод Рунге-Кутта четвертого порядка:

y(x + h) = y(x) + k1 + 2k2 + 2k3 + k4 ;   k1 = f(x,y)• h;

6

k2 = f(x + y2,y + ky£).h;  k3 = f(x + y2,y + h;  k4 = f(x + h,y + k3)• h.

Запишем программу, реализующую эти два метода; обратите внимание на два последних оператора программы: примерно за одно и то же время метод Эй­лера достигает лишь погрешности 1E-3, в то время как метод Рунге-Кутты дает 1E-5.

TYPE Func2Type = FUNCTION (x,y:Real):Real;

FUNCTION Euler(x0,y0:Real; XX:Real; f:Func2Type; Epsilon:Real):Real; VAR n,i : LongInt; h,x,y,Y1,Y2,k1,d : Real; BEGIN n:=1; h:=XX-x0; Y1:=y0+f(x0+h/2,y0+f(x0,y0)*h/2)*h; REPEAT n:=n ShL 1; h:=(XX-x0)/n; y:=y0;

FOR i:=1 TO n DO BEGIN

x:=x0+(i-1)*h; k1:=f(x,y)*h; y:=y+f(x+h/2,y+k1/2)*h; END;

Y2:=y; d:=Abs(Y1/Y2-1); UNTIL d<Epsilon; Euler:=Y2;

END;

FUNCTION Runge_Kutta(x0,y0:Real; XX:Real; f:Func2Type; Epsilon:Real):Real; VAR n,i : LongInt; h,x,y,Y1,Y2,k1,k2,k3,k4,d : Real;

BEGIN n:=1; h:=XX-x0;

k1:=f(x0,y0)*h; k2:=f(x0+h/2,y0+k1/2)*h; k3:=f(x0+h/2,y0+k2/2)*h; k4:=f(XX,y0+k3)*h;

Y1:=y0+(k1+2*(k2+k3)+k4)/6;

REPEAT n:=n ShL 1; h:=(XX-x0)/n; y:=y0;

FOR i:=1 TO n DO BEGIN x:=x0+(i-1)*h; k1:=f(x,y)*h; k2:=f(x+h/2,y+k1/2)*h;

k3:=f(x+h/2,y+k2/2)*h; k4:=f(x+h,y+k3)*h; y:=y+(kl+2*(k2+k3)+k4)/6;

END;

Y2:=y; d:=Abs(Yl/Y2-l); UNTIL d<Epsilon; Runge_Kutta:=Y2;

END;

FUNCTION F(x,y:Real):Real; FAR; BEGIN F:=y/2/x+l/Sqrt(x); END;

BEGIN WRITELN;

WRITELN(Euler(2,0.980258l,2.5,F,lE-3)); WRITELN(Runge_Kutta(2,0.980258l,2.5,F,lE-5));

END.