Эксперименты лаба5,6(2курс)

#i, j переменные циклов

#xp середины частичных отрезков

#n количество частичных отрезков

#x массив узлов

#h длины частичных отрезков n:=1; x[1]:=a;

for i from 1 while (x[i]<=b) do n:=n+1;

x[i+1]:=x[1]+sh*i; # заполнение массива узлов end do;

if x[n]<b then x[n]:=b; end if; for i from 1 to n-1 do

#вычисление длин частичных отрезков h[i]:=abs(x[i+1]-x[i]);

#нахождение середин частичных отрезков xp[i]:=h[i]/2+x[i];

end do; slog:=0;

for j from 1 to n-1 do

slog:=slog+simplify(h[j]*f(xp[j])); # формула прямоугольников end do;

slog:=evalf(slog);

return slog; # вывод результата end proc;

Проверим результат работы процедуры > prl(0,1,0.0001);

0.5222818376

> evalf(int((sin(x))-xˆ 3*ln(x),x=0..1));

0.5221976941

4.9.2.Метод трапеций

>restart;

>f:=x-> evalf(sin(x)-xˆ 3*ln(x));

# вычисление определенного интеграла методом трапеций > trp:=proc(a,b,sh)

11

#b, a верхний и нижний пределы интегрирования соответственно

#sh шаг разбиения отрезка [a,b]

local i,j,n,x,y,h,slog;

#i, j переменные циклов

#n количество частичных отрезков

#x массив узлов

#y массив значений функции в полученных узлах

#h длины частичных отрезков n:=1; x[1]:=a; y[1]:=f(x[1]);

for i from 1 while (x[i]<=b) do n:=n+1;

x[i+1]:=x[1]+sh*i; # заполнение массива узлов

#заполнение массива значений функции в полученных узлах

y[i+1]:=f(x[i+1]); end do;

if x[n]<b then x[n]:=b; y[n]:=f(x[n]); end if; for i from 1 to n-1 do

h[i]:=abs(x[i+1]-x[i]); # вычисление длин частичных отрезков end do;

slog:=0;

for j from 1 to n-1 do

slog:=slog+evalf(h[j]*(y[j]+y[j+1])); # формула трапеций end do;

slog:=evalf(1/2*slog);

return slog; # вывод результата end proc;

Проверим работу процедуры > trp(0.001,1,0.001);

0.5230383120

> evalf(int((sin(x))-xˆ 3*ln(x),x=0..1));

0.5221976941

4.9.3.Метод Симпсона

>restart;

>f:=x-> evalf(sin(x)+ln(x)*xˆ 3);

12

> Simpson:=proc(a,b,n)

#b, a верхний и нижний пределы интегрирования соответственно

#n количество частичных отрезков

local i,h,x,y,slog,slogn,slogch;

#i переменная циклов

#h шаг

#x массив узлов

#y массив значений функции в полученных узлах

#slogn сумма всех значений функции в узлах с нечетным номером

#slogch сумма всех значений функции в узлах с четным номером (кроме x[0] и x[n])

if (n mod 2 = 1) then print(’n-нечетное_число’); exit; end if; h:=abs(evalf(b)-evalf(a))/n; # определение шага, для заполнения

массива узлов

for i from 0 to n do

x[i]:=evalf(a)+i*h; # заполнение массива узлов

#заполнение массива значений функции в полученных узлах y[i]:=f(x[i]);

end do;

slog:=0; slogn:=0; slogch:=0; for i from 1 by 2 to (n-1) do

#вычисление суммы всех значений функции в узлах с нечетным номером

slogn:=slogn+y[i]; end do;

for i from 2 by 2 to (n-2) do

#вычисление суммы всех значений функции в узлах с четным номером (кроме x[0] и x[n])

slogch:=slogch+y[i]; end do;

slog:=(y[0]+4*slogn+2*slogch+y[n])*h/3; # формула Симпсона return slog; # вывод результата

end proc;

Проверим работу процедуры > Simpson(0.0001,Pi,22);

23.78877518

> int(f(x),x=0.0001..Pi);

23.78870622

13

4.9.4.Метод Гаусса

Таблица узлов и коэффициентов формулы Гаусса для n=1..6

>restart;

>n:=6:

#табличные значения узлов формулы Гаусса для n=6 > x[1]:=-0.9324695142:

> x[2]:=-0.6612093865: > x[3]:=-0.2386191861: > x[4]:=-x[3]:

> x[5]:=-x[2]: > x[6]:=-x[1]:

#табличные значения коэффициентов формулы Гаусса для n=6 > A[1]:=0.1713244924:

> A[2]:=0.3607615730: > A[3]:=0.4679139346: > A[4]:=A[3]:

> A[5]:=A[2]: > A[6]:=A[1]:

#процедура для вычисления интеграла методом Гаусса с 6 узлами > Gauss6:=proc(a,b)

#b, a верхний и нижний пределы интегрирования

local i,sl,t;

# i переменная циклов

14

# sl переменная для накопления суммы слагаемых

#t массив узлов for i from 1 to n do

t[i]:=(b+a)/2+((b-a)/2)*x[i]; # заполнение массива узлов end do;

sl:=0;

for i from 1 to n do sl:=sl+evalf(A[i]*f(t[i])); # формула Гаусса end do;

sl:=evalf((b-a)/2*sl); # вывод результата end proc;

#процедура для вычисления определенного интеграла (метод Гаусса) с заданной точностью eps

> Gauss:=proc(a,b,eps)

#b, a верхний и нижний пределы интегрирования

#eps точность интегрирования

local int1,int2,int,i,j,c,k,r,l;

#int1, int2 два последовательных приближения к точному решению интеграла

#int приближения к точному решению интеграла

#i, j переменные циклов

#c точки разбиения промежутка интегрирования

#k количество частичных отрезков

#r, l правые и левые границы частичных отрезков

#вычисление определенного интеграла (метод Гаусса n=6) int1:=Gauss6(a,b);

c:=(b-a)/2; # точка разбиения промежутка интегрирования

#вычисление определенного интеграла (метод Гаусса n=11) int2:=Gauss6(a,c)+Gauss6(c,b);

#количество частичных отрезков, на которое рабивается отрезок

[a,b]

k:=2;

for j from 1 while abs(int1-int2)>eps do int:=0;

c:=(b-a)/(2*k);

for i from 0 to 2*k-1 do

#вычисление границ каждого частичного отрезка

l:=a+c*i; r:=a+c*(i+1);

# вычисление определенного интеграла (n=6+2*k-1)

15

int:=int+Gauss6(l,r); end do;

int1:=int2;

int2:=int;

k:=k*2; end do;

return evalf(int2); end proc;

# задаем подынтегральную функцию

> f:=x-> sin(x)*xˆ 2*ln(x)/(xˆ 4*exp(x)ˆ 2);

Проверим работу процедур. > Gauss(1,8,0.0000000001);

0.008185355232

> evalf(int(f(x),x=1..8));

0.008185355234

> Gauss6(1,8);

0.008710562730

16