Практикум по программированию Обработка числовых данных

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

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

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

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

5.1.Вычисление определенного интеграла

ñзаданной точностью

Задача вычисления определенного интеграла формулируетс я сле-

b

дующим образом: вычислить I = ò f (x)dx c точностью e при изве-

a

стных значениях пределов интегрирования a, b, известной точности e и заданной подынтегральной функции f(x). При вычислении интеграла с точностью e будут использоваться изложенные ранее методы прямоугольников, трапеций и парабол для определен ия значения интеграла при заданном числе разбиений интерва ла интегрирования. Для оценки погрешности вычисления интегра ла на практике используют правило Рунге [4]. Суть правила состоит в том, что выполняют вычисление интеграла с двумя разными ш а- гами изменения переменной x, а затем сравнивают результаты и получают оценку точности. Наиболее часто используемое пр авило связано с вычислением интеграла дважды: с шагом Dx и шагом Dx/2.

Для методов прямоугольников и трапеций погрешность R x/2 вычисления интеграла с шагом Dx/2 оценивается следующей формулой:

ãäå I x/2 — значение интеграла, вычисленное с шагом Dx/2; I x — значение интеграла, вычисленное с шагом Dx.

Для метода Симпсона погрешность оценивается в соответствии со следующим выражением:

Шаг интегрирования для методов прямоугольников и трапе-

ций пропорционален ε, поэтому в [4] рекомендовано начальное число разбиений выбирать согласно следующему выраже -

В программе вычисления интеграла с точностью e во внутреннем цикле находят значение определенного интеграла п ри

заданном (фиксированном) числе разбиений интервала инте грирования. Во внешнем цикле производится сравнение значений и нтегралов, вычисленных при числе шагов, равных n и 2n соответственно. Если требуемая точность не достигнута, то число р азбиений удваивается, а в качестве предыдущего значения интегр ала берут текущее и вычисление интеграла выполняется при нов ом числе разбиений.

В дальнейшем в программах при вычислении значений аргумента используют формулу арифметической прогрессии X = a + + (i − 1) x, а не формулу накопления суммы X = X + x, так как в первом случае меньше погрешность вычислений [9].

Пример программы вычисления определенного интеграла с заданной точностью методом трапеций:

program integraltrap;

{Вычисление определенного интеграла с заданной точностью методом трапеций. Подынтегральная функция - X*Exp(X). Первообразная функция (x-1)*Exp(X).} {$APPTYPE CONSOLE}

uses SysUtils, Math;

function Rus(S:String):String; var I:Byte;

begin Result:='';

for I:=1 to Length(S) do case S[I] of

'À'..'ï': Result:=Result+Chr(Ord(S[I])-64); 'ð'..'ÿ': Result:=Result+Chr(Ord(S[I])-16); '¨': Result:=Result+Chr(240);

'¸': Result:=Result+Chr(241); else

Result:=Result+S[I];

end;

end;

var

A,B,Eps,I1,I2,Itoch,X,Dx,S1:Real;

I,N,M,K:Integer; begin

WriteLn(Rus('Введите пределы интегрирования '), Rus('и точность'));

ReadLn(A,B,Eps);

//вычисление начального количества разбиений N:=Trunc((B-A)/Sqrt(Eps))+1;

//полусумма значений функции на нижнем и // верхнем пределах интегрирования S1:=(A*Exp(A) + B*Exp(B))/2.0;

//начальное значение интеграла, вычисленное //при начальном количестве разбиений Dx:=(B-A)/N;

I2:=0;

for I:=1 to N-1 do begin

X:=A+I*Dx; // текущее значение аргумента //вычисление суммы значений функции //в узлах интегрирования

I2:=I2+X*Exp(X);

end;

I2:=(S1+I2)*Dx;

//вычисление количества позиций при выводе //числа с заданной точностью M:=Trunc(-Log10(Eps))+2; K:=Trunc(Log10(Abs(I2)))+M+3; Writeln(Rus('Начальное количество разбиений=')

,N:4,Rus(' Интеграл='),I2:K:M); //цикл вычисления интеграла с точностью repeat

//Переменной I1 присваивается текущее //приближенное значение интеграла I1:=I2;

//Переменной I2 присваивается полусумма //значений функции на концах интервала //интегрирования

I2:=S1;

//Удвоение количества разбиений N:=N*2;

//Шаг изменения переменной интегрирования Dx:=(B-A)/N;

//Внутренний цикл для вычисления суммы //значений функции при фиксированном //количестве шагов

for I:=1 to N-1 do begin

//Вычисление текущего значения аргумента X:=A+I*Dx;

//Накопление суммы значений функции //в узлах интегрирования I2:=I2+X*Exp(X);

end;

//Текущее значение интеграла I2:=I2*Dx;

//Анализ точности вычисления until Abs(I2-I1)/3<Eps;

WriteLn(Rus('Значение интеграла ='),I2:K:M ,Rus(' при количестве разбиений, равном ') , N);

//Точное значение интеграла Itoch:=(B-1)*Exp(B)-(A-1)*Exp(A); WriteLn(Rus('Точное значение интеграла равно ')

,Itoch:K:M);

ReadLn;

end.

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

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

ReadLn(A,B,Eps); N:=Trunc((B-A)/Sqrt(Eps))+1; S1:=(A*Exp(A) + B*Exp(B))/2.0; S:=0;

Dx:=(B-A)/N;

//Цикл вычисления суммы значений функции //в узлах интегрирования при начальном //разбиении интервала интегрирования. for I:=1 to N-1 do

begin X:=A+I*Dx;

S:=S+ X*Exp(X); end; I2:=Dx*(S+S1);

S2:=S; repeat

I1:=I2;

N:=N*2; Dx:=(B-A)/N;

//Координата первого нового узла интегрирования. X1:=A+Dx;

S:=0;

//Цикл вычисления суммы значений функции //в новых узлах интегрирования.

for I:=1 to N div 2 do begin

X:=X1+2*(I-1)*Dx;

S:=S+ X*Exp(X); end;

//Сумма значений функции в узлах интегрирования //при новом числе разбиений.

S2:=S+S2;

I2:=(S1+S2)*Dx; until Abs(I2-I1)/3<Eps;

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

Xi = Xíà÷ + x i/6 + x(i − 1), i = 1, n; Xj = Xíà÷ + 5 x j/6 + x( j − 1), j = 1, n. В приведенных выражениях Xíà÷ — нижний предел интег-

рирования, x — шаг, n — количество разбиений на предыдущей итерации.

При вычислении интеграла по методу Симпсона можно непосредственно использовать суммы значений функции, вычи с- лявшиеся на предыдущей итерации. При этом надо иметь в вид у, что узлы с нечетными номерами получают на очередной итера - ции четные индексы, и сумма значений функции для этих узло в должна браться с коэффициентом 2 (на предыдущей итерации с коэффициентом 4). Узлы с четными номерами и при новой итерации останутся четными, поэтому сумма значений в этих уз лах берется с тем же коэффициентом 2.

5.2.Задания для самостоятельной работы

Âтаблице приведены подынтегральные функции. Вычислить значение интеграла каждым из трех методов (прямоугольник ов, трапеций, парабол) с заданной точностью. Пределы интегрир о- вания брать из интервала, установленного заданием. Сравнить скорость сходимости используемых методов. Обеспечить в п рограмме ввод точности, пределов интегрирования, коэффицие н- тов, присутствующих в функциях, вывод с поясняющим тексто м вычисленных значений интеграла каждым методом, количест во разбиений при использовании каждого метода.

5.3. Вычисление наибольшего (наименьшего) значения функции с заданной точностью на заданном интервале

Необходимость организации вложенного цикла возникает т акже и при решении задачи нахождения наибольшего (наименьшего) значения функции на заданном интервале с заданной т оч- ностью. Рассмотрим сначала задачу определения с заданной точ- ностью аргумента, при котором функция достигает своего экстремального значения.

Данная задача может быть сформулирована следующим образом. Заданы некоторая функция f (x) и некоторый интервал [a, b]. Известно, что на заданном интервале функция имеет один экстремум, известен и вид экстремума (максимум или миниму м). Требуется с заданной точностью ε найти значение аргумента, при котором достигается экстремум функции.

Решить поставленную задачу можно, используя прием программирования — вычисление максимума или минимума. При этом шаг изменения аргумента следует задать равным требу емой точности ε. Однако такой подход может потребовать большого количества вычислений. Сократить количество выполняемы х операций можно за счет использования следующего алгоритма. С на- чала вычисляют значения функции при грубом значении шага изменения аргумента. При этом очередное значение функции сравнивают с ранее вычисленным максимальным значением (д о-

пустим, что определяется максимум функции); если текущее з на-

чение превышает максимальное, т. е. f (xi) > fmax , òî fmax := f (xi), и производится вычисление значения функции в следующей точ-

ке. Если же на очередном шаге будет выполнено условие

f (xi) < fmax , то это означает, что максимум функции уже пройден, т. е. находится на интервале (xi − 2h, xi), где h — шаг изменения

аргумента. В этом случае шаг изменения аргумента уменьшае тся (обычно в два раза) и на вновь полученном интервале вычисл я- ют значения функции и максимум при новом шаге изменения аргумента. Процесс продолжается до тех пор, пока h не станет меньше или равным ε.

Таким образом, во внутреннем цикле вычисляют значения функции и максимум на заданном интервале изменения аргумента при заданном шаге. Во внешнем цикле задают новый интервал поиска максимума и новый шаг изменения аргумента. Фрагмент программы вычисления максимума функции y = 2x3 + 10x2 + 6x − 20 в интервале [a, b] имеет вид

//Вычисление максимума функции с использованием //итерационного (вложенного) цикла

N2:=0; //Количество вычислений значения функции ReadLn(A,B,Eps,H);

Xn:=A;

while H>Eps do begin

I:=0;

//Установка начального значения для максимума Ymax:=((2*Xn+10)*Xn+6)*Xn-20;

Xmax:=Xn;

//Цикл определения максимума функции //на очередном интервале изменения аргумента repeat

I:=I+1;

//Вычисление текущего значения аргумента X:=Xn+(I-1)*H;

//Вычисление значения функции в очередной точке Y:=((2*X+10)*X+6)*X-20;

N2:=N2+1;

//Определение нового максимума //и соответствующего значения аргумента if Y>Ymax then

begin Ymax:=Y; Xmax:=X;

end;

until (Y<Ymax);

//Определение левого конца нового интервала //вычисления максимума

Xn:= Xmax-H;

//Уменьшение шага изменения аргумента в два раза H:=H/2;

end;

Поиск минимума функции можно свести к поиску максимума, если функцию f (x) заменить функцией −f (x).

Рассматриваемую задачу обобщим на случай, когда функция не имеет внутри заданного интервала экстремум, а изменяет ся внутри него монотонно, причем заранее характер поведения функции не известен. Тогда следует говорить о нахождении наи - большего (наименьшего) значения функции на заданном интервале. Наибольшее (наименьшее) значение функции достиг а- ется либо в критической точке (где первая производная рав на нулю), либо на конце интервала.

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

until (Y<Ymax)or(X>=B);

//Определение левой границы нового интервала //вычисления наибольшего значения.

Xn:= Xmax-H;

//Проверка выхода нового значения аргумента //за пределы заданного интервала.

if Xn<A then Xn:=A;

Если известны аналитические выражения для первой и второй производных, то можно найти значение аргумента, при ко тором первая производная обращается в нуль (решить уравнен ие, например, методом деления отрезка пополам), и определить з нак второй производной в найденной точке. Знак второй произво д- ной позволит определить вид экстремума (максимум или мин и- мум) или его отсутствие на заданном интервале. При отсутст вии экстремума следует сделать запрос о нахождении наибольш его или наименьшего значения функции.