Вычисление интеграла.

← →
ezorcist (2006-12-17 14:34) [0]

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

пишу:

lim:=20; sum:=0; step:=(xk-xn)/lim; x1:=xn; repeat x2:=x1+step; x01:=func(x1); x02:=func(x2); If abs(x01-x02) > eps then begin sum:=0; inc(lim,10); step:=(xk-xn)/lim; x1:=xn; x2:=0; end else begin sum:=sum+func(x1)*step; x1:=x2; end; until x2>=xk;

ранее написана сама function func;

Вопрос: почему работате так медленно и можно ли сделать быстрей? при достаточно резко возрастающей функции и точности eps=0.001 можно смело идти обедать...

заранее спасибо

← →
Desdechado © (2006-12-17 17:56) [1]

А ты прикинь, сколько раз при шаге 0.001 будет вычислена твоя FUNC.
Тем паче, что в sum:=sum+func(x1)*step она зачем-то повторно вычисляется. Да и сама FUNC неизвестная. Мож, ты там вес вселенной вычисляешь.

← →
Думкин © (2006-12-17 19:51) [2]

Можно.
Если резко возрастает, то у тебя шаг начинает подыскиваться через уменьшение. В итоге, если функция с большим значением производной то, этот шаг становится очень маленьким. Если производная порядка 1000, то шаг в тысячу раз меньше eps. Но этого мало, он до этого мало, ты к этому шагу идешь черепашьими щагами - увеличивая число разбиений не умножением, а на определенное количество. В итоге чтобы уменьшить шаг в 10 раз, тебе надо сделать 18 попыток. Чтобы в 100 раз ...И на каждом шаге считать функцию, притом от 2 до 3 раз.

← →
TUser © (2006-12-17 22:39) [3]

Тут на быстро возрастающих/убывающих участках вся функция начинает считаться с начала.

← →
Bolt (2006-12-17 23:04) [4]

Для начала АЛГОРИТМ немного поправить надо, а потом перенести на Assembler

← →
Kostafey © (2006-12-17 23:56) [5]

Я когда-то писал так. Притом за чаем сбегать не успевал:
unit UTypes; interface type Treal = Real; Preal = ^Treal; TFunc = function (x:Treal):Treal of object; TFunc2= function (x:Treal;y:Treal):Treal; TDouble2 = array[1..2] of Double; TFunc3= function (x:Treal;y:Treal;dy:Treal):TDouble2; implementation end.
unit UIntegrMaethod; interface uses UTypes, SysUtils; type TIntegerMethod=class(TObject)//Методы интегрирования: //Методом левых прямоугольников function IntegralLeftRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; //Методом правых прямоугольников function IntegralRightRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; //Метод центральных прямоугольников function IntegralRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; //Метод трапеций function IntegralTrap(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; //Метод Симпсона function IntegralSimps(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; procedure PrintN(n:integer); end; implementation uses Math, Ap, Unit1; //Интегрирование методом правых прямоугольников function TIntegerMethod.IntegralRightRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; var i : Integer; n : Integer; h : Double; s1 : Double; s2 : Double; begin Form1.Gauge1.Progress:=0; n := 1; h := b-a; s2 := h*F((a+b)/2); repeat n := 2*n; s1 := s2; h := h/2; s2 := 0; i := 1; Form1.Gauge1.Progress:=Form1.Gauge1.Progress+3; repeat s2 := s2+F(a+h+h*(i-1)); i := i+1; until not (i<=n); s2 := s2*h; until not (AbsReal(s2-s1)>3*Epsilon); Form1.Gauge1.Progress:=100; Result := s2; PrintN(n); end; //Интегрирование методом левых прямоугольников function TIntegerMethod.IntegralLeftRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; var i : Integer; n : Integer; h : Double; s1 : Double; s2 : Double; begin Form1.Gauge1.Progress:=0; n := 1; h := b-a; s2 := h*F((a+b)/2); repeat n := 2*n; s1 := s2; h := h/2; s2 := 0; i := 1; Form1.Gauge1.Progress:=Form1.Gauge1.Progress+3; repeat s2 := s2+F(a+h*(i-1)); i := i+1; until not (i<=n); s2 := s2*h; until not (AbsReal(s2-s1)>3*Epsilon); Form1.Gauge1.Progress:=100; Result := s2; PrintN(n); end; //Интегрирование методом центральных прямоугольников function TIntegerMethod.IntegralRect(F:TFunc; a : Double; b : Double; Epsilon : Double):Double; var i : Integer; n : Integer; h : Double; s1 : Double; s2 : Double; begin Form1.Gauge1.Progress:=0; n := 1; h := b-a; s2 := h*F((a+b)/2); repeat n := 2*n; s1 := s2; h := h/2; s2 := 0; i := 1; Form1.Gauge1.Progress:=Form1.Gauge1.Progress+3; repeat s2 := s2+F(a+h/2+h*(i-1)); i := i+1; until not (i<=n); s2 := s2*h; until not (AbsReal(s2-s1)>3*Epsilon); Form1.Gauge1.Progress:=100; Result := s2; PrintN(n); end; //.....прочие методы end.

> точность понимается разница ординат двух соседних точек
- не понимаю логику такого определения, а вдруг на отрезке [x<sub>k</sub>;x<sub>k+1</sub>] функция быстро возросла, но потом опять вернулась обратно? разве точность, это не δ >= |a - A|, т.е. максимальное отклонение от точного значения, т.е. "хочу знать интеграл с такой-то точностью"? Ну или хотя бы разница абсцисс соседних точек, т.е. h = (b-a)/n - шаг. Но если всё же определять точность как разницу ординат, то вместо увеличения числа отрезков:

> inc(lim,10);

предлагаю делить "плохие" отрезоки надвое, подобно
> y1:=abs(func(x1)); > repeat > x2:=x1+step; > y2:=abs(func(x2)); > If abs(y1-y2) > eps then > sum:=sum+subdivide(x1,x2,y1,y2,step) > else > sum:=sum+y1*step; > x1:=x2; > y1:=y2; > until x2>=xk; > > function subdivide(x1,x3,y1,y3,step: double): double; > var > x2,y2: double; > begin > x2 := (x1+x3)/2; > y2 := abs(func(x2)); > if abs(y1-y2) <= esp then > result:=y1*step > else > result:=subdivide(x1,x2,y1,y2,step / 2); > if abs(y2-y3) <= esp then > result:=y2*step > else > result:=subdivide(x2,x3,y2,y3,step / 2); > end
жжж

> Kostafey © (17.12.06 23:56) [5]
> i := 1; > repeat … > i := i+1; > until not (i<=n);

- почему не for i := 1 to n do… ?

1. Обычно количество отрезков не увеличивают на 10, а удваивают. Тогда правило Рунге, то есть оценка точности через разность двух последовательных шагов, применяется более обоснованно.
2. При методе прямоугольников с удвоением есть возможность не считать повторно значения функций в точках, в которых функция уже считалась на предыдущем шаге. Нужно просто тупо разделить интегральную сумму на два, а потом добавить к ней сумму по добавившимся точкам.

Вычисление интеграла. Найти похожие ветки