Гауссов интеграл ошибок матлаб

Описание

f = @(z) 1./(2.*z-1);
contour_segments = [1+1i 0+1i 0-1i 1-1i];
q = quadgk(f,1,1,'Waypoints',contour_segments)

fun = @(x) sin(2e4*pi*x);
[Q,errbnd] = quadgk(fun,0,pi)

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is  5.7e-01. The integral may not exist, or it may be difficult to approximate numerically. Increase MaxIntervalCount to 1272 to enable QUADGK to continue for another iteration.

[Q,errbnd] = quadgk(fun,0,pi,'MaxIntervalCount',1e5)

errbnd <= max(AbsTol,RelTol*abs(q))

errbnd <= max(AbsTol,RelTol*abs(q))

errbnd <= max(AbsTol,RelTol*abs(q))

Представленный в R2007b

Write a Matlab function that applies the Gauss three point rule to N sub-intervals of $[a, b].$ The input parameters should be the name of the function being integrated, $a, b,$ and $N$.

Attempt:

For $int^b_a f(x) dx$ the Gaussian 3 point quadrature is given by:

$$frac{b-a}{18} left( 5f left( frac{a+b}{2}- sqrt{frac{3}{5}} frac{b-a}{2} right) +8f left( frac{a+b}{2} right) + 5 f left( frac{a+b}{2}+ sqrt{frac{3}{5}} frac{b-a}{2} right) right).$$

Here is my Matlab code that uses this equation to approximate $int^b_a f(x) dx$:

function t = gauss3(func,a,b)
j=(((a+b)/2)-sqrt(3/5)*((b-a)/2));
k=((a+b)/2);
l=(((a+b)/2)+sqrt(3/5)*((b-a)/2));
t=((b-a)/18)*(5*feval(func,j) + 8*feval(func,k) + 5*feval(func,l))

function y=f(x)
y=exp(-x^2);
> t = gauss3(@f,0,1)

So how can I modify this code to take into account the sub-intervals $N$? The function should start as function t = gauss3(func,a,b,N).

P.S. In the code above I am trying to approximate the integral $int^1_0 e^{-x^2} dx,$ using the Gaussian three point rule with error $<10^{-10}.$

Edit:

function y = gaussq(func,a,b,N);

h = (b-a)/N;
y = 0;
for i = 1:N,
    y = y+gauss3(func,a+(i-1)*h,a+i*h);
end
% y0 = 100;
% for N = logspace(0,6,7),
%     y1 = gaussq(@f,a,b,N);
%     err = abs((y1-y0)/y0);
%     if err < 1.0e-10,
%         break;
%     end
%     y0 = y1;
% end

Where:

function y=f(x)
y=exp(-x^2);

Write a Matlab function that applies the Gauss three point rule to N sub-intervals of $[a, b].$ The input parameters should be the name of the function being integrated, $a, b,$ and $N$.

Attempt:

For $int^b_a f(x) dx$ the Gaussian 3 point quadrature is given by:

$$frac{b-a}{18} left( 5f left( frac{a+b}{2}- sqrt{frac{3}{5}} frac{b-a}{2} right) +8f left( frac{a+b}{2} right) + 5 f left( frac{a+b}{2}+ sqrt{frac{3}{5}} frac{b-a}{2} right) right).$$

Here is my Matlab code that uses this equation to approximate $int^b_a f(x) dx$:

function t = gauss3(func,a,b)
j=(((a+b)/2)-sqrt(3/5)*((b-a)/2));
k=((a+b)/2);
l=(((a+b)/2)+sqrt(3/5)*((b-a)/2));
t=((b-a)/18)*(5*feval(func,j) + 8*feval(func,k) + 5*feval(func,l))

function y=f(x)
y=exp(-x^2);
> t = gauss3(@f,0,1)

So how can I modify this code to take into account the sub-intervals $N$? The function should start as function t = gauss3(func,a,b,N).

P.S. In the code above I am trying to approximate the integral $int^1_0 e^{-x^2} dx,$ using the Gaussian three point rule with error $<10^{-10}.$

Edit:

function y = gaussq(func,a,b,N);

h = (b-a)/N;
y = 0;
for i = 1:N,
    y = y+gauss3(func,a+(i-1)*h,a+i*h);
end
% y0 = 100;
% for N = logspace(0,6,7),
%     y1 = gaussq(@f,a,b,N);
%     err = abs((y1-y0)/y0);
%     if err < 1.0e-10,
%         break;
%     end
%     y0 = y1;
% end

Where:

function y=f(x)
y=exp(-x^2);

Системы
MatLabSciLab обладают большими возможностями
программирования и визуализации
результатов научных исследований при
решении задач вычислительной математики.
Численные методы в этих системах
предлагаются уже в виде готового
инструментария — набора функций, их
реализующих. Задачей данной работы
является ознакомление с наиболее часто
встречающимися на практике численными
методами, которые используются при
решении нелинейных уравнений, при
аппроксимации и интегрировании функций.
Приближѐнное вычисление определѐнных
интегралов. С вычислением определѐнных
интегралов самого общего вида:

приходится
сталкиваться при решении многих задач
математики, физики, химии и т.д. Как
правило, на практике получить первообразную
аналитически и затем воспользоваться
формулой Ньютона-Лейбница удаѐтся
крайне редко. В этих случаях интеграл
вычисляется приближѐнно с помощью
одного из методов численного интегрирования.
В методе трапеций используется линейная
интерполяция, т.е. подынтегральная
функция заменяется кусочно-линейной,
а еѐ график представляется в виде
ломаной, соединяющей отдельные точки
yi
= f(xi),
где xi
– координаты узлов сетки, введѐнной
на интервале интегрирования. При этом,
очевидно, на каждом из интервалов от
xi
до xi+1
искомые площади вычисляются как площади
обычных трапеций, после чего суммируются.
В итоге формула трапеций при постоянном
шаге сетки h
= xi+1
— xi
= (b
— a)/n
(n
– число интервалов разбиения сетки)
принимает вид:

Эта
схема интегрирования чрезвычайно проста
для реализации и имеет 2-ой порядок
точности относительно шага сетки. Тем
не менее, для экономии времени пользователя
в MatLab
еѐ реализовали в виде функции trapz(x,
y),
а в SciLab
– inttrap(x,
y).
Здесь x
и y
– массивы одного и того же размера 1хn,
содержащие, соответственно, координаты
узлов выбранной сетки и значения функции
в этих узлах. Не многим сложнее
оказывается формула Симпсона 4-го
порядка точности, получаемая при кусочной
интерполяции подынтегрального выражения
с помощью парабол:

В
MatLab формула Симпсона реализуется
функцией quad(fun, a, b [,tol]), где fun –
подынтегральная функция, a, b – концы
интервала интегрирования, tol – требуемая
точность (необязательный параметр).
Подынтегральная функция может быть
предварительно определена в m-файле,
тогда еѐ имя следует указать с дескриптором
@, или еѐ можно указать непосредственно
при вызове quad в виде текстовой строки,
заключѐнной в одинарные кавычки.
Например, вызов quad(@myfun, 0, 1, 1e-3) вычислит
интеграл от myfun(x) на интервале [0, 1] с
точностью 0.001. Файл функции myfun может
быть примерно следующим: function m = myfun(x)

m
= x.*x;

return

Аналогично,
производится интегрирование при
непосредственном задании функции,
например: quad(‘x.*x’,
0, 1, 1e-3).
Примечание: не забывайте, что операции
умножения, деления и возведения в степень
необходимо указывать в форме с «.», т.к.
они производятся над массивами.

Интегрирование
методом трапеций.

Приведенные
ниже функции выполняют численное
интегрирование хорошо известным методом
трапеций и методом трапеций с накоплением.

•
trapz(Y)
возвращает определенный интеграл,
используя интегрирование методом
трапеций с единичным шагом между
отсчетами. Если Y
– вектор, то trapz(Y)
возвращает интеграл элементов вектора
Y,
если Y
– матрица, то trapz(Y)
возвращает векторстроку, содержащую
интегралы каждого столбца этой матрицы.

•
trapz(X,Y)
возвращает интеграл от функции Y
по переменной X,
используя метод трапеций (пределы
интегрирования в этом случае задаются
начальным и конечным элементами вектора
X).

•
trapz(…,dim)
возвращает интеграл по строкам или по
столбцам для входной матрицы, в зависимости
от значения переменной dim.

Примеры:

>>
y=[1,2,3,4]

y
= 1 2 3 4

>>
trapz(y)

ans
= 7.5000

>>
X=0:pi/70:pi/2; Y=cos(X); Z = trapz(Y)

Z
= 22.2780

•
cumtrapz(Y)
возвращает численное значение
определенного интеграла для функции,
заданной ординатами в векторе или
матрице Y
с шагом интегрирования, равным единице
(интегрирование методом трапеций с
накоплением). В случае когда шаг отличен
от единицы, но постоянен, вычисленный
интеграл достаточно умножить на величину
шага. Для векторов эта функция возвращает
вектор, содержащий результат интегрирования
с накоплением элементов вектора Y.
Для матриц – возвращает матрицу того
же размера, что и Y,
содержащую результаты интегрирования
с накоплением для каждого столбца
матрицы Y.
cumtrapz(X,Y)
выполняет интегрирование с накоплением
от Y
по переменной X,
используя метод трапеций. X
и Y
должны быть векторами одной и той же
длины, или X
должен быть векторомстолбцом, а Y
– матрицей. cumtrapz(…,
dim)
выполняет интегрирование с накоплением
элементов по размерности, точно
определенной скаляром dim.
Длина вектора X
должна быть равна size(Y,dim).

Интегрирование
методом квадратур

Метод
трапеций обеспечивает невысокую точность
при заданном числе шагов или дает слишком
большое число шагов при вычислениях с
заданной погрешностью. Приведенные
ниже функции осуществляют интегрирование
и двойное интегрирование, используя
более точную квадратурную формулу
Симпсона или метод Гаусса–Лобатто.
Квадратура – численный метод нахождения
площади под графиком функции f(x),
то есть вычисление определенного
интеграла вида. В приведенных ниже
формулах подынтегральное выражение
fun
обычно задается или в прямых апострофах,
или в форме handleфункции.
Функции quad
и quadl
используют два различных алгоритма
квадратуры для вычисления определенного
интеграла. Функция quad
выполняет интегрирование по методу
низкого порядка, используя рекурсивное
правило Симпсона. Но она может быть
более эффективной при негладких
подынтегральных функциях или при низкой
требуемой точности вычислений. Алгоритм
quad
в MATLAB
6 изменен по сравнению с предшествовавшими
версиями, точность по умолчанию по
сравнению с версиями 5.3x
повышена в 1000 раз (с 10–3 до 10–6). Новая
функция quad1
(квадратура Лобатто) использует адаптивное
правило квадратуры Гаусса–Лобатто
очень высокого порядка. Устаревшая
функция quad8
выполняла интегрирование, используя
квадратурные формулы Ньютона–Котеса
8-го порядка. Достижимая точность
интегрирования гладких функций в MATLAB
6 поэтому также значительно выше, чем в
предшествующих версиях. quad(fun,a,b)
возвращает численное значение
определенного интеграла от заданной
функции @fun
на отрезке [a
b].
Используется значительно усовершенствованный
в MATLAB
6 адаптивный метод Симпсона. quad(fun,a,b,tol)
возвращает численное значение
определенного интеграла с заданной
относительной погрешностью tol.
По умолчанию tol=1.e–6.
Можно также использовать вектор,
состоящий из двух элементов tol
=[rel_tol
abs_tol],
чтобы точно определить комбинацию
относительной и абсолютной погрешностей.
quad(fun,a,b,tol,trace)
возвращает численное значение
определенного интеграла и при значении
trace,
не равном нулю, строит график, показывающий
ход вычисления интеграла.
quad(fun,a,b,tol,trace,P1,P2,…)
возвращает численное значение
определенного интеграла по подынтегральной
функции fun,
использует дополнительные аргументы
P1,
P2,
…, которые напрямую передаются в
подынтегральную функцию: G=fun(X,P1,P2,…).

Вычисления
двойных и тройных интегралов

Для
вычисления двойных интегралов служит
следующая функция:

•
dblquad(fun,inmin,inmax,outmin,outmax)
вычисляет и возвращает значение двойного
интеграла для подынтегральной функции
fun(inner,outer),
по умолчанию используя квадратурную
функцию quad.
inner
– внутренняя переменная, изменяющаяся
от inmin
до inmax,
а outer
– внешняя переменная, изменяющаяся от
outmin
до outmax.
Первый аргумент @fun
– строка, описывающая подынтегральную
функцию. Эта функция должна быть функцией
двух переменных вида fout=fun
(inner,outer).
Функция должна брать вектор inner
и скаляр outer
и возвращать вектор fout,
который является функцией, вычисленной
в outer
и в каждом значении inner.

•
dblquad(fun,inmin,inmax,outmin,outmax,tol,trace)
передает в функцию dblquad
параметры tol
и trace.
Смотрите справку по функции quad
для получения информации о параметрах
tol
и trace.

•
dblquad(fun,inmin,inmax,outmin,outmax,tol,trace,order)
передает параметры tol
и trace
для функции quad
или quadl
в зависимости от значения строки order.
Допустимые значения для параметра order
– @quad
, @quadl
или имя любого определенного пользователем
квадратурного метода с таким же вызовом
и такими же возвращаемыми параметрами,
как у функций quad
и quadl.
(Например, при проверке старых программ
можно использовать @quad8
для большей совместимости с прежними
версиями MATLAB).
По умолчанию (без параметра order)
вызывается @quad,
поскольку подынтегральные функции
могут быть не гладкими.

Пример:
пусть mфайл
integ1.m
описывает функцию 2*y*sin(x)+x/2*cos(y),
тогда вычислить двойной интеграл от
этой функции можно следующим образом:

>>
result = dblquad(@integ1,pi,2*pi,0,2*pi)

result
= -78.9574

В
систему MATLAB
6.5 была введена новая функция для
вычисления тройных интегралов. Она
имеет четыре формы записи:

triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax)

triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)

triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method)

triplequad(fun,xmin,xmax,ymin,ymax,zmin,zmax,tol,method,p1,p2,…)

Параметры
этой функции уже были определены выше
(добавлены только пределы для новой
переменной z).
Пример:

>>
Q
= triplequad(‘x+y*z’,0,pi,0,1,-1,1,0.001)

Q
= 9.8696

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

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]

• #
• #
• #
• #
• #
• #
• #
• #
• #
• #
• #