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

49. program integraltoch;

50. {вычисление определенного интеграла с заданной точностью методом трапеций}

51. { подынтегральная функция – X*Exp(X); первообразная функция (x-1)*Exp(X)}

52. {$APPTYPE CONSOLE}

53. uses

54. SysUtils,Math;

55. var

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

57. I,N,M,K:Integer;

58. begin

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

60. ReadLn(A,B,Eps);

61. //вычисление начального количества разбиений

62. N:=Trunc((B-A)/Sqrt(Eps))+1;

63. //начальное значение интеграла, взятое при n=1

64. I2:=(A*Exp(A)+B*Exp(B))*(B-A)/2.0;

65. //полусумма значений функции на нижнем и

66. // верхнем пределах интегрирования

67. S1:=(A*Exp(A) + B*Exp(B))/2.0;

68. //начало внешнего цикла вычисления интеграла

69. repeat

70. //предыдущему значению интеграла присваивается текущее

71. I1:=I2;

72. //начальному значению суммы присваивается полусумма

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

74. I2:=S1;

75. //шаг изменения переменной интегрирования

76. Dx:=(B-A)/N;

77. //внутренний цикл для вычисления суммы значений

78. //функции при фиксированном количестве шагов

79. for I:=1 to N-1 do

80. begin

81. X:=A+I*Dx; // текущее значение аргумента

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

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

84. end;

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

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

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

88. //вычисление количества позиций при выводе

89. //числа с заданной точностью

90. M:=Trunc(-Log10(Eps))+2;

91. K:=Trunc(Log10(I2))+M+3;

92. WriteLn('Значение интеграла =',I2:K:M

93. ,' при количестве разбиений, равном ‘, N div 2);

94. Itoch:=(B-1)*Exp(B)-(A-1)*Exp(A);//точное значение интеграла

95. WriteLn('Точное значение интеграла равно ',Itoch:K:M);

96. ReadLn;

97. end.

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

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

98. ReadLn(A,B,Eps);

99. N:=Trunc((B-A)/Sqrt(Eps))+1;

100. S1:=(A*Exp(A) + B*Exp(B))/2.0;

101. S:=0;

102. Dx:=(B-A)/N;

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

104. //вания при начальном разбиении интервала интегрирования

105. for I:=1 to N-1 do

106. begin

107. X:=A+I*Dx;

108. S:=S+ X*Exp(X);

109. end;

110. I2:=Dx*(S+S1);

111. S2:=S;

112. repeat

113. I1:=I2;

114. N:=N*2;

115. Dx:=(B-A)/N;

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

117. S:=0;

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

119. // интегрирования

120. for I:=1 to N div 2 do

121. begin

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

123. S:=S+ X*Exp(X);

124. end;

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

126. //числе разбиений

127. S2:=S+S2;

128. I2:=(S1+S2)*Dx;

129. until Abs(I2-I1)/3<Eps;

Таким же образом можно сократить количество выполняемых операций при вычислении интеграла по методу левых или правых прямоугольников. Однако на практике при использовании метода прямоугольников значение функции вычисляют не на конце интервала (левом или правом), а в его середине. В этом случае при удвоении числа разбиений ранее вычисленные значения функции (и их сумму) использовать не удается. Для использования результатов ранее выполнявшихся вычислений при применении метода прямоугольников следует количество разбиений увеличивать в три раза и вычислять значения функции в точках с абсциссами X_i = X_нач +hi/6+h(i-1), i=1,n; X_j = X_нач +5hj/6+h(j-1), j=1,n. В приведенных выражениях X_нач – нижний предел интегрирования, h – шаг интегрирования и n – количество разбиений на предыдущей итерации.

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