Функция Power предназначена для возведения значения Xx в целую степень P. Функция Deg_Rad предназначена для перевода значения угла Gr из градусной

меры в радианы.

4.4 Примеры использования предопределенных процедур.

Пример №1 Отделение корней путем анализа графика функции.

Использование некоторых подпрограмм модуля Librery1 рассмотрим на примере нахождения корня уравнения:

f(x) =k1 x −cos(k2 x) =0

(3.26)

где k1 и k2 - коэффициенты функции.

Сначала произведем отделение корней с использованием графического представления функции на экране монитора. Для этого разработана программа Sample_1.

В этой программе вслед за заголовочным оператором помещается ключевое слово uses после которого располагается список имен модулей к процедурам которых будут обращения в теле программы. В данном случае это модули Librery1 и Crt. Следует обратить также внимание на использование в данной программе определения типа переменных Vec, принятого в модуле Library1. При определении функции F используются глобальные переменные K1 и K2 отражающие значения коэффициентов уравнения. Поэтому соответствующие формальные и фактические параметры, отражающие значения коэффициентов уравнения носят одинаковые имена. Такой прием обеспечивает определенные удобства, но в некоторых случаях опасен (когда программа громоздка, а пользователь недостаточно компетентен). Все формальные параметры функции F – параметры - значения, поэтому ключевое слово var перед списком формальных параметров не ставится.

В программе после ввода исходных данных производится (в цикле) расчет значений аргумента и функции. Затем осуществляется обращение к процедуре Grafic. Обратите внимание на то, что при обращении к процедуре Grafic формальные параметры, отражающие маркеры и точки излома линии графика совпадают. Следовательно, на графике маркеры будут расположены в точках излома. К выбору количества точек на графике следует подходить осторожно. Большое количество точек (а следовательно и маркеров) затруднит рассмотрение графика, малое количество плохо отразит особенности функции. Разумным представляется назначить N = 20 – 50.

Ниже приводится текст программы Sample_1.pas, протокол ввода данных и график функции выводимый на экран монитора.

program Sample_1; (* Рисование графика функции *) uses Library1, Crt;

var

N,I : Integer; X,Y : Vec;

K1,K2 : Real; Xmin,Xmax,Dx : Real; function F(XX:Real):Real; begin

F:=K1*XX-cos(K2*XX);

100

end;
begin
ClrScr;	Коэффициенты функции:');
Writeln('
Write('	k1=');	Readln(K1);
Write('	k2=');	Readln(K2);
Write('	Количество точек на оси X - N='); Readln(N);
Write('	xmin=');	Readln(Xmin);
Write('	xmax=');	Readln(Xmax);
Dx:=(Xmax-Xmin)/(N-1);
for I:=1 to N do
begin
X[I]:=Xmin+Dx*(I-1);
Y[I]:=F(X[I]);
end;
Grafic(N,X,Y,N,X,Y);
Readln;
end.
Протокол ввода данных при выполнении программы Sample_1:
Коэффициенты функции:
k1=	0.200
k2=	2.000
Количество точек на оси X - N=50
xmin=	0.000
xmax= 10.000

Рис. 3.4 Изображение на экране монитора графика выводимого программой Sample_1

101

Как видно из рис. 3.4 уравнение (3.26) в положительной части оси x имеет три корня x1 ≈ 0.7, x2 ≈ 2.7, x3 ≈ 3.5.

Пример № 2. Уточнение значения корня методом Ньютона.

В качестве примера демонстрирующего уточнение значения корня уравнения (3.26) методом Ньютона рассмотрим программу Sample_2.pas. В этой программе производится обращение к процедуре Newton определенной в модуле Librery1. Будем искать уточненное значение только первого корня так, как уточнение остальных корней производится аналогично.

Обратим внимание на определение функции Fdivdf, которая отражает отношение функции (3.26) к ее производной. В нашем случае:

	f	=	( k1 x − cos( k 2	x )	(3.27)
	f '		( k1 + k 2 sin( k 2	x )

Поскольку в приведенной программе используется обращение к функции (Fdivdf), определенной пользователем в самой программе, как к параметру, то необходимо задать директиву компилятору {$F+}, определяющую модель дальнего вызова процедур и функций.

program Sample_2; (* Поиск корня уравнения методом Ньютона *) uses Library1,Crt;

var

K1,K2,X,E : Real; {$F+}

function Fdivdf(X:real):Real; begin

Fdivdf:=(K1*X-cos(K2*X))/(K1+K2*sin(K2*X)); end;

begin

ClrScr;

Writeln('Значения коэффициентов функции:');

Write('	k1=');	Readln(K1);
Write('	k2=');	Readln(K2);

Write('Начальное приближение X0='); Readln(X); Write('Абсолютный итерационный допуск e='); Readln(E); Newton(X,E,Fdivdf);

Writeln('Значение корня X=',X:8:3); Readln;

end.

Протокол обмена программы Sample_2:

Значения коэффициентов функции: k1=0.2

k2=2.0

Начальное приближение X0=0.700 Абсолютный итерационный допуск e=0.005

Значение корня X=

0.714

102

Пример № 3. Уточнение значения корня методом половинного деления.

В программе Sample_3.pas уточнение значения корня производится с использованием метода половинного деления (процедура Dix определенная в модуле Library1).

program Sample_3; uses Library1,Crt; var

Kod : Integer; A,B,K1,K2,X,E : Real;

{$F+}

function F(X:Real):Real; begin

F:=K1*X-cos(K2*X); end;

begin

ClrScr;

Writeln('Границы интервала:');

Write('	a=');	Readln(A);
Write('	b=');	Readln(B);

Writeln('Коэффициенты функции:');

Write('	K1=');	Readln(K1);
Write('	K2=');	Readln(K2);

Write('Абсолютный итерационный допуск e='); Readln(E); Dix(A,B,E,X,F,Kod);

if Kod=0 then Writeln('Значение корня X=',X:8:3)

else Writeln('Значение корня вне заданного интервала!'); Readln;

end.

Протокол обмена программы Sample_3:

Границы интервала: a= 0.00

b= 1.00

Коэффициенты функции: K1= 0.20

K2= 2.00

Абсолютный итерационный допуск e=0.00500 Значение корня X= 0.715

Пример № 4. Уточнение значения корня методом хорд.

В программе Sample_4.pas уточнение значения корня производится с использованием метода хорд (процедура Chord определенная в модуле Library1).

program Sample_4; (* Нахождение корня методом хорд *) uses Library1,Crt;

var

Kod : Integer;

103

A,B,K1,K2,X,E : Real; {$F+}

function F(X:Real):Real; begin

F:=K1*X-cos(K2*X); end;

begin

ClrScr;

Writeln('Границы интервала:');

Write('	a=');	Readln(A);
Write('	b=');	Readln(B);

Writeln('Коэффициенты функции:');

Write('	K1=');	Readln(K1);
Write('	K2=');	Readln(K2);

Write('Абсолютный итерационный допуск e='); Readln(E); Chord(A,B,E,X,F,Kod);

if Kod=0 then Writeln('Значение корня X=',X:8:3) else begin

Writeln('*** Нет корня в заданном интервале'); Writeln('или на интервале несколько корней! ***');

end; Readln;

end.

Протокол обмена программы Sample_4:

Границы	интервала:
a=	0.00
b=	1.00

Коэффициенты функции:

K1=	0.20
K2=	2.00

Абсолютный итерационный допуск e=0.00500 Значение корня X= 0.714

Пример № 5. Решение СЛАУ методом Гаусса с выбором главного элемента.
Рассмотрим СЛАУ :
2 x1 + 1 x2 + 1 x3 = 7
1 x1 + 3 x2 + 1 x3 = 10	(3.28)

1 x1 + 2 x2 + 3 x3 = 14

Для нахождения решения рассмотренной СЛАУ будем использовать процедуру Gauss, определенную в модуле Library1. Для удобства ввода значений коэффициентов матрицы и вектора правых частей используем процедуру VvodA, определенную в том же модуле. Проверку решения осуществим путем сравнения правых частей уравнений (заданных и полученных путем подстановки решения в уравнения СЛАУ).

program Sample_5; uses Library1,Crt; var

104

I,J,N :		Integer;
A :		Mat;
X	:	Vec;
S,P	:	Real;

begin

ClrScr;

VvodA(N,A);

Gauss(N,A,X,S);

if S=0.0 then begin

Write('Матрица вырожденная'); Readln;

Exit;

end else begin

Writeln('Решение СЛАУ методом Гаусса.'); for I:=1 to N do

Writeln('	x[',I:2,']=',X[I]:10:3);
end;	Проверка	***');
Writeln('***
for I:=1 to N	do begin

P:=0.;

for J:=1 to N do P:=P+A[I,J]*X[J];

Writeln('i=',I:2,'			задано b=',A[I,N+1]:8:3,
end;	'	получено b=',P:8:3);
Readln;
end.
Протокол обмена программы Sample_5:
Ввод данных для решения СЛАУ
Порядок СЛАУ		N= 3		a[i,j]
Ввод значений коэффициентов матрицы
и элементов		вектора правой части		b[i)]
-------------------------------------------------
i - индекс		строки,	j - индекс столбца;
-------------------------------------------------
a[ 1, 1]=	2.000
a[ 1, 2]=	1.000
a[ 1, 3]=	1.000	1]= 7.000
a[ 2, 1]=	b[
	1.000
a[ 2, 2]=	3.000
a[ 2, 3]=	1.000	2]= 10.000
a[ 3, 1]=	b[
	1.000
a[ 3, 2]=	2.000
a[ 3, 3]=	3.000	3]= 14.000
	b[
Решение СЛАУ методом Гаусса.
x[ 1]=	1.000
x[ 2]=	2.000

105

x[ 3]=

3.000

***Проверка ***

i= 1		задано b=		7.000	получено b=		7.000
i=	2	задано	b=	10.000	получено	b=	10.000
i=	3	задано	b=	14.000	получено	b=	14.000

Пример № 6. Решение СЛАУ методом Зейделя.

Для решения СЛАУ (3.28) итерационным методом Зейделя воспользуемся процедурой Zeidel, определенной в модуле Library1. Для проверки решения путем сравнения правых частей уравнений (заданных и полученных путем подстановки решения в уравнения СЛАУ) воспользуемся процедурой Testing, определенной в том же модуле...

program Sample_6; uses Library1,Crt; var

I,J,N,M : Integer;

A : Mat;

X : Vec;

E,P : Real; begin

ClrScr;

Write('Порядок матрицы N='); Readln(N); for I:=1 to N do begin

for J:=1 to N do begin

Write('a[',I:2,',',J:2,']='); Readln(A[I,J]); end;

Write('	b[',I:2,']='); Readln(A[I,N+1]);
end;

Write('Абсолютный итерационный допуск e='); Readln(E); Zeidel(N,A,X,E,M);

Writeln('Решение СЛАУ методом Зейделя.'); Writeln('Число итераций =',M:4);

for I:=1 to N do Writeln(' x[',I:2,']=',X[I]:10:3); Testing(N,A,X);

Readln; end.

Протокол обмена программы Sample_6:

Порядок матрицы N=3

a[ 1, 1]=	2.
a[ 1, 2]=	1.
a[ 1, 3]=	1.
b[ 1]=	7.
a[ 2, 1]=	1.
a[ 2, 2]=	3.
a[ 2, 3]=	1.
b[ 2]=	10.
a[ 3, 1]=	1.
a[ 3, 2]=	2.
a[ 3, 3]=	3.

106

	b[ 3]= 14.
Абсолютный итерационный допуск e=0.005
Решение СЛАУ методом Зейделя.
Число итераций = 6
x[ 1]=				1.001
x[ 2]=				1.999
x[ 3]=				3.000	значений вектора правых частей СЛАУ
Проверка по невязке
i=1		заданное b=			7.000
i=2	полученное b=				7.001
		заданное b=			10.000
i=3	полученное b=				9.999
		заданное b=			14.000
	полученное b=				14.000
Пример № 7. Аппроксимация табулированной функции полиномом.
	Рассмотрим табулированную функцию:
i		1		2		3		4	5	6	7	8
xi		0		0.15		0.27		0.43	0.52	0.65	0.81	1.0
yi		0		0.0.31		0.52		0.63	0.8	0.85	0.98	1.0

Необходимо данный набор значений отразить (аппроксимировать) полиномом второго порядка. Для аппроксимации воспользуемся процедурами Gram и Gauss, определенными в модуле Library1. Процедура Gram формирует матрицу Грамма, а процедура Gauss предназначена для решения СЛАУ методом Гаусса с выбором главного элемента. Кроме вывода полученных значений коэффициентов аппроксимирующего полинома в программе процедурой Grafic обеспечивается вывод результатов в виде графика. Причем заданные значения табулированной функции отражаются на графике маркерами, а результат аппроксимации - ломаной.

program Sample_7; uses Library1,Crt; var

I,J,K,L,N,M : Integer; Xmin,Xmax,Dx,Dy,S : Real;

X,Y,Xg,Yg : Vec; A : Mat; Koef : Vec;

begin

ClrScr;

Write(' Количество точек табулированной функции n='); Readln(N);

for I:=1 to N do begin

Write(' x[',I:2,']='); Read(X[I]);

Write(' y[',I:2,']='); Readln(Y[I]);

if (I=1) or (Xmin>X[i]) then Xmin:=X[i]; if (I=1) or (Xmax<X[i]) then Xmax:=X[i];

107

end;

Write(' Степень аппроксимирующего полинома k='); Readln(K);

M:=K+1;

Gram(N,A,X,Y,M,K);

Gauss(M,A,Koef,S); if S=0 then begin

Write('Матрица вырожденная!'); Readln;

Exit;

end else begin

Writeln('Коэффициенты аппроксимирующего полинома.'); Writeln(' y = a0 + a1*x + a2*x^2 +...+ ak*x^k');

for I:=1 to M do begin

Writeln(' a[',(i-1):1,']=', Koef[I]:10:4);

end; end;

Writeln('Результаты аппроксимации табулированной функции полиномом.');

for I:=1 to N do begin

S:=0.0;

for J:=1 to M do S:=S+Koef[J]*Power(X[I],(J-1)); Dy:=S-Y[I];

Writeln('X=',X[I]:8:3,'; Yтаб=',Y[I]:8:3, '; Yaпр=', S:8:3,'; dY=',Dy:8:5);

end; Readln;

(* Подготовка данных для построения графика *) Dx:=(xmax-xmin)/100.;

for I:=1 to 101 do begin

S:=0.0; xg[i]:=xmin+dx*(i-1);

for J:=1 to M do S:=S+Koef[J]*Power(Xg[I],(J-1)); Yg[i]:=S;

end; Grafic(N,X,Y,101,Xg,Yg);

end.

Протокол обмена программы Sample_7:

Количество точек табулированной функции n=8 x[ 1]=0.

y[ 1]=0.

x[ 2]=0.15

y[ 2]=0.31

x[ 3]=0.27

y[ 3]=0.52

x[ 4]=0.43

y[ 4]=0.63

x[ 5]=0.52

y[ 5]=0.80

108

	x[ 6]=0.65
		y[ 6]=0.85
	x[ 7]=0.81
		y[ 7]=0.98
	x[ 8]=1.00
		y[ 8]=1.00
	Степень аппроксимирующего полинома k=2
Коэффициенты аппроксимирующего полинома.
y = a0 + a1*x + a2*x^2 +...+ ak*x^k
	a[0]=	0.0210
	a[1]=	1.9579
	a[2]=	-0.9799
Результаты аппроксимации табулированной функции полиномом.
X=	0.000;	Yтаб=	0.000;	Yaпр=	0.021;	dY= 0.02098
X=	0.150;	Yтаб=	0.310;	Yaпр=	0.293;	dY=-0.01738
X=	0.270;	Yтаб=	0.520;	Yaпр=	0.478;	dY=-0.04183
X=	0.430;	Yтаб=	0.630;	Yaпр=	0.682;	dY= 0.05168
X=	0.520;	Yтаб=	0.800;	Yaпр=	0.774;	dY=-0.02589
X=	0.650;	Yтаб=	0.850;	Yaпр=	0.880;	dY= 0.02958
X=	0.810;	Yтаб=	0.980;	Yaпр=	0.964;	dY=-0.01607
X=	1.000;	Yтаб=	1.000;	Yaпр=	0.999;	dY=-0.00107

Графическое представление результатов аппроксимации на экране монитора приводится на рис.3.5.

Рис.3.5 Результат аппроксимации табулированной функции полиномом.

109

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

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

Пусть решению подлежит система алгебраических уравнений:

2 x12 + x 22 + 0.5 x 32 = 6.5;
x13	+ 3 x 2 + 0.5 x 33 = 7.5;	(3.29)
x12 + 0.5 x 2 + 3 x 32 = 5.
Для данной системы известно точное решение x1=1, x2=2, x3=1.
Систему нелинейных уравнений (3.29) можно привести к виду:
a11 x1 + a12 x12 + a13 x3 = a14 ;
a21 x1	+ a22 x12 + a23 x3 = a24 ;	(3.30)
a31 x1	+ a32 x12 + a33 x3 = a34

где	a11 = 2 x1;	a12 = x2; a13 = 0.5 x3;		a14 = 6.5;
	a21 = x12;	a22 = 3 ;	a23 = 0.5 x32;	a24 = 7.5;	(3.31)
	a31 = x1;	a32 = 0.5;	a33=3 x3;	a34 =5 .
Решение системы (3.30) будем вести итерационным методом. В начале зададимся

значениями x1, x2, x3 в первом приближении. Затем на каждой итерации будем рассчитывать значения коэффициентов согласно (3.31) и решать систему (3.30) методом Гаусса с главным элементом. Итерационный процесс завершим тогда, когда для каждого искомого значения xi выполнится условие:

xim − xim−1

≤ ε,

(3.32)

где ε – абсолютный итерационный допуск, m – номер итерации.

В некоторых случаях итерационный процесс может сходиться очень медленно или вообще быть расходящимся. Для обеспечения устойчивости и быстрой сходимости

итерационного процесса часто используют релаксационную формулу (3.33):
	xim = ω xi + (1- ω) xim-1,	(3.33)
где xi	– решение, полученное на m-ой текущей итерации методом Гаусса,
xim-1 – решение, фиксированное на предыдущей итерации,
xim	- решение, принятое на текущей итерации,

ω - параметр релаксации; его значения должны находиться в диапазоне 0 ..1. Таким образом, в качестве результата на текущей итерации принимается некоторое значение лежащее между значением принятым на предыдущей итерации и решением, полученным методом Гаусса на текущей итерации. Чем меньше значение параметра релаксации, тем принятое решение ближе к результату предыдущей итерации. Если параметр релаксации равен 1, то влияние предыдущего приближения не учитывается. Изменяя значение ω можно в определенной мере управлять сходимостью итерационно-

го процесса решения системы уравнений вида (3.30).

Оговоренный подход реализован в программе Sample_7.pas. Константами в программе задаются: порядок матрицы N=3 и значение параметра релаксации W=0.5.

В программе определяется процедура Mm, предназначенная для пересчета значений коэффициентов матрицы. Для учета сходимости итерационного процесса вводится переменная Ig. В начале каждого прохода цикла Repeat – Until (итерации) значение Ig обнуляется. В процессе счета, если условие сходимости (3.32) не выполняется, к значению Ig прибавляется единица. Таким образом, только в том случае если для всех иско-

110

мых переменных выполнится условие (3.32), значение Ig к концу цикла останется равным нулю и можно будет завершить итерационный процесс.

Массив Xx отражает значения, принимаемые в качестве решения. Массив X отражает решение, полученное методом Гаусса.

В протоколе обмена приведены результаты итерационного решения.

program Sample_8; uses Library1,Crt; const N:Integer=3;

W:Real=0.5;

var

Iter,Ig,I,J : Integer; A : Mat;

XX,X : Vec;

E,S,P : Real;

procedure M(Var A:Mat; Var X:Vec); begin

A[1,1]:=2.*X[1]*X[1]; A[1,2]:=X[2]; A[1,3]:=0.5*X[3]; A[1,4]:=6.5; A[2,1]:=X[1]*X[1]; A[2,2]:=3.0; A[2,3]:=0.5*X[3]*X[3]; A[2,4]:=7.5; A[3,1]:=X[1]; A[3,2]:=0.5; A[3,3]:=3.0*X[3]; A[3,4]:=5;

end;

begin

ClrScr;

Writeln('Ввод начальных значений искомых величин:'); Write('X[1]='); Readln(Xx[1]);

Write('X[2]='); Readln(Xx[2]);

Write('X[3]='); Readln(Xx[3]); Write('Итерационный допуск E='); Readln(E); Iter:=0;

repeat

Iter:=Iter+1;

Ig:=0;

M(A,Xx);

Gauss(N,A,X,S);

for I:=1 to N do begin

if abs(Xx[I]-X[I])>E then Ig:=Ig+1; Xx[I]:=W*X[I]+(1.-W)*Xx[I];

end; until Ig=0;

Writeln('* Решение * Количество итераций =',Iter:3); for I:=1 to N do begin

111

Writeln(' x[',I:2,']=',X[I]:10:3); end;

Readln; end.

Протокол обмена программы Sample_8:

Ввод начальных значений искомых величин: X[1]=0.5

X[2]=0.5

X[3]=0.5

Итерационный допуск E=0.005

* Решение * Количество итераций = 21 x[1]= 1.003

x[2]= 1.996 x[3]= 0.999

112