Использование методов Симпсона и Гаусса-Лобатто

В системе MATLAB методы численного интегрирования более высокого порядка точности реализуются функциями quad (метод Симпсона) и quad1 (метод Гаусса-Лобатто). Следует отметить, что обе указанные функции реализуют адаптивные алгоритмы, т.е. отсутствует необходимость контроля со стороны пользователя достигнутой точности вычисления интеграла посредством сопоставления результатов, полученных при использовании разных шагов интегрирования (эту процедуру адаптивные функции реализуют самостоятельно). Заметим, что функция quad1 имеет более высокий порядок точности по сравнению с функцией quad, что позволяет выполнять интегрирование гладких функций с достаточно большими шагами интегрирования при сохранении точности (т.е., по сути, при меньшем объёме вычислительной работы). Функции quad и quad1 имеют достаточно длинные списки формальных параметров (как, впрочем, и многие другие функции в системе MATLAB), передаваемых в функцию при её вызове. Минимальный список формальных параметров данных функций включает в себя три параметра: дескриптор подынтегральной функции, нижний и верхний пределы интегрирования:

q=quad (fun, a, b)

q1=quad1 (fun, a, b)

где fun — указатель на подынтегральную функцию; a и b — верхний и нижний пределы интегрирования.

Четвёртый (дополнительный) параметр списка, задающий требуемую точность результата вычислений, является необязательным (по умолчанию точность принимается равной 10^-6).

Вычислим теперь рассмотренный выше интеграл, используя функции quad и quad1. Имеем:

>> f=inline ('1-x+x.^2-x.^3', 'x');

>> TolX=1e-5;

>> [I, cnt]=quad (f, a, b, TolX);

-12.749999999999996

cnt=

>> [I,cnt]=quad1 (f, a, b, TolX)

-12.750000000000002

cnt=

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

>> f=inline ('sin(x)', 'x');

>> TolX=1e-5;

>> [I, cnt]=quad (f, -pi, pi, TolX);

2.220446049250313e-16

cnt=

>> [I,cnt]=quad1 (f, -pi, pi, TolX)

2.015274060089494e-17

cnt=

т.е. в данном случае уже функция quad1 более эффективна (установленная точность была достигнута за 18 интеграций против 21 для функции quad).

Выполнение работы на ЭВМ.

Для выполнения создается М-файл. Ниже приведен текст М-файла.

function Lab_3_5

a=input('Введите a=');

b=input('Введите b=');

eps=input('Введите eps=');

nmax=input('Введите nmax=');

fprintf('\n Значения интеграла S по разным методам:');

% Вычисление интегралов с использованием функции MATLAB

f=@(x)p(x);

squad=quad(f,a,b);

fprintf('\n - Симпсона (функция quad), S=%10.3f',squad);

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

[spr,npr]=INT_Pr(a,b,eps,nmax);

fprintf('\n - прямоугольников, S=%10.3f',spr);

fprintf(' (n=%4d)',npr);

[str,ntr]=INT_Tr(a,b,eps,nmax);

fprintf('\n - трапеций, S=%10.3f',str);

fprintf(' (n=%4d)',ntr);

[ssimps,nsimps]=INT_Simps(a,b,eps,nmax);

fprintf('\n - Симпсона, S=%10.3f',ssimps);

fprintf(' (n=%4d)',nsimps);

end

function [s,nf]=INT_Pr(a,b,eps,nmax)

n=4;

s1=pr(a,b,n);

n=2*n;

s2=pr(a,b,n);

while(abs(s2-s1)>eps&n<nmax)

s1=s2;

n=2*n;

s2=pr(a,b,n);

end

s=s2;

nf=n;

end

function [s,nf]=INT_Tr(a,b,eps,nmax)

n=4;

s1=tr(a,b,n);

n=2*n;

s2=tr(a,b,n);

while(abs(s2-s1)>eps&n<nmax)

s1=s2;

n=2*n;

s2=tr(a,b,n);

end

s=s2;

nf=n;

end

function [s,nf]=INT_Simps(a,b,eps,nmax)

n=4;

s1=simps(a,b,n);

n=2*n;

s2=simps(a,b,n);

while(abs(s2-s1)>eps&n<nmax)

s1=s2;

n=2*n;

s2=simps(a,b,n);

end

s=s2;

nf=n;

end

function y=pr(a,b,n)

h=(b-a)/n;

x=a+h/2;

s=0;

for i=1:n

s=s+p(x);

x=x+h;

end

y=h*s;

end

function y=tr(a,b,n)

h=(b-a)/n;

x=a;

s=(p(a)+p(b))/2;

for i=2:n

x=x+h;

s=s+p(x);

end

y=h*s;

end

function y=simps(a,b,n)

h=(b-a)/n;

x=a+h;

s=p(a)+p(b);

z=1;

for i=2:n

s=s+(3+z)*p(x);

z=-z;

x=x+h;

end

y=h/3*s;

end

function y=p(x)

y=1-x+x.^2-x.^3;

end

Результаты расчета в командном окне при a=0, b=3, eps=0.001,

nmax=1000:

Значения интеграла S по разным методам:

- Симпсона (функця quad), S= -12.750

- прямоугольников, S= -12.750 (n= 256)

- трапеций, S= -12.750 (n= 256)

- Симпсона, S= -12.750 (n= 8)

Ниже приведены пояснения к приведенному тексту программы.

1. В системе MATLAB методы численного интегрирования высокого порядка точности реализуются, в частности, функцией quad, в основе которой лежит метод Симпсона. В данную функцию заложен адаптивный алгоритмы, т.е. отсутствует необходимость контроля со стороны пользователя достигнутой точности вычисления интеграла посредством сопоставления результатов, полученных при использовании разных шагов интегрирования (эту процедуру адаптивная функция реализует самостоятельно). Функция quad имеет достаточно длинный список формальных параметров, передаваемых в функцию при ее вызове. Минимальный список формальных параметров данных функций включает в себя три параметра: дескриптор подынтегральной функции, нижний и верхний пределы интегрирования:

q=quad(fun,a,b), где fun – указатель на подынтегральную функцию (чтобы создать указатель на функцию перед именем функции следует добавить символ @); a и b – верхний и нижний пределы интегрирования. Четвертый (дополнительный) параметр списка, задающий требуемую точность результата вычислений, является необязательным (по умолчанию точность принимается равной ). .

2. В приведенной программе реализованы все три рассматриваемые в настоящем курсе формулы численного интегрирования: в функциях pr и INT_Pr – вычисление интеграла по формуле прямоугольников; в функциях tr и INT_Tr – вычисление интеграла по формуле трапеций; в функциях simps и INT_Simps – вычисление интеграла по формуле Симпсона.

3. Переменная (в программе – nmax) в представленной программе (а также в последующих программах, относящихся к данной

лабораторной работе) определяет количество отрезков, на которые делится интервал интегрирования для вычисления интеграла с заданной точностью.

При этом значение n последовательно удваивается до тех пор, пока изменение приближенного значения интеграла не станет меньше заданного малого числа (в программе – eps), значение которого вводится. При расчете было принято: nmax=1000, eps=0.001.