Ирина Эланс

Автор который поможет с любыми образовательными и учебными заданиями

Приближенное вычисление интеграла

Содержание

Введение 2

1. Различные методы вычисления определенных интегралов 3

1.1. Метод Симпсона для интегрирования функций F(x) по

заданному промежутку и его реализация на языке Pascal 4

1.2. Метод Симпсона для интегрирования функции от двух

переменных F(x,y) по прямоугольной двумерной области и его

реализация на языке Pascal 5

1.3. Метод Ромберга и его реализация на языке Pascal 7

1.4. Метод Гаусса и его реализация на языке Pascal 10

Заключение 16

Литература 17

Введение

Система программирования Турбо Паскаль представляет собой единство двух в известной степени самостоятельных начал: компилятора с языка программирования Паскаль (язык назван в честь выдающегося французского математика и философа Блеза Паскаля (1623-1662)) и некоторой инструментальной программной оболочки, способствующей повышению эффективности создания программ.

Паскаль – гибкий и развитый в отношении типов данных язык. Привлекательны его рекурсивные возможности, а также поддержка технологии объектно-ориентированного программирования.

Изучение программирования на языке Паскаль может дать хороший старт в огромный и увлекательный мир программирования. Обучение языку программирования проходит намного более эффективно с изучением примеров.

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

Различные методы вычисления определенных интегралов.

Приближенное вычисление интеграла,

I =

Основано на его замене конечной суммой:

I_n =

k F(x_k),

где w_k – числовые коэффициенты, а x_k – точки отрезка [x₀,x₁]. Приближенное равенство

I ≈ I_n называется квадратурной формулой, точки x_k – узлами квадратурной формулы, а числа w_k – коэффициентами квадратурной формулы. Разные методы приближенного интегрирования отличаются выбором узлов коэффициентов. От этого выбора зависит погрешность квадратурной формулы.

R_n =‌

‌.

В модуле integral реализовано несколько методов численного интегрирования как для простых (одномерных), так и для кратных (многомерных) интегралов.

В функции simpson реализован стандартный метод Симпсона для интегрирования функции F(x) по заданному промежутку, когда число разбиений интервала выбирается заранее. Функция double_simpson является прямым обобщением метода Симпсона на случай интегрирования функции от двух переменных F(x,y) по прямоугольной двумерной области.

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

В функции romberg запрограммирован еще один метод адаптивного интегрирования – метод Ромберга, в настоящее время, вероятно, один из наиболее популярных. Имеются также функция gauss – одномерная версия метода интегрирования Гаусса. Интерфейсная секция модуля integral приведена в листинге 1.1.

Листинг 1.1. Интерфейсная секция модуля integral.

Unit integral;

Interface

Const

Max_dim =10;

Max_deg=96;

Type

Real_fun=function(x:real):real;

Real_fun2=function(x,y:real):real;

Real_vec=array[1..max_dim+1] of real;

Index=array[1..max_dim+1] of word;

Vec_fun=function(j:word; x:real_vec):real;

Var

no_evaluations, highest_level:word;

function simpson(F:real_fun; x0,x1:real; div_no:word):real;

function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;

function adaptive_simpson(F:real_fun;x0,x1,eps,eta:real):real;

function romberg(f:real_fun; x0,x1,eps,eta:real; min,max:word):real;

function gauss3(F:real_fun;x0,x1:real; n:word):real;

procedure compute_gauss_coeffs(deg:word);

function gauss(Freal_fun:x0,x1:real; deg:word):real;

1.1. Метод Симпсона для интегрирования функций F(x) по заданному промежутку и его реализация на языке Pascal.

Перейдем к секции реализации. Она начинается описанием функции simpson. Стоит сказать несколько слов о выборе узлов и коэффициентов квадратурной формулы Симпсона. Идея трехточечного метода Симпсона заключается в следующем.

Пусть x_m – это средняя точка интервала [x₀, x₁] и пусть Q(x) – единственный полином второй степени, который интерполирует (приближает) подынтегральную функцию F(x) по точкам x₀, x_m и x₁. Искомый интеграл аппроксимируется интегралом от функции Q(x):

I≈

Это оценка точна, если F(x) является полиномом степени 3.

В функции simpson интервал интегрирования делится на div_no равных частей, а трехточечная формула Симпсона применяется к каждому такому интервалу. Параметрами функции simpson (листинг1.2) являются, по порядку, подынтегральная функция, нижняя и верхняя границы интервала интегрирования и количество подынтервалов.

Листинг 1.2. Функция simpson модуля integral

Implementation

Uses crt;

Var

zero, weight:array [1..max_deg] of real;

Function simpson(F:real_fun; x0,x1:real; div_no:word):real;

Var

x, dx, sum:real;

j:word;

begin

dx:=(x1-x0)/(2.0*div_no);

sum:=F(x0)+F(x1);

x:=x0;

for j:=1 to 2*div_no-1 do

begin

x:=x+dx;

if Odd (j) then

sum:=sum+4.0*F(x)

else

sum:=sum+2.0*F(x);

end;

simpson:=dx*sum/3.0;

end;

1.2. Метод Симпсона для интегрирования функции от двух переменных F(x,y) по прямоугольной двумерной области и его реализация на языке Pascal.

Функция double_simpson (листинг 1.3.) является, по существу, прямым обобщением одномерного метода Симпсона на случай вычисления двойного интеграла по прямоугольной области.

Листинг 1.3. Функции double_simpson и simple_simpson модуля integral

Function double_simpson(F:real_fun2; x0,x1,y0,y1:real; x_div,y_div:word):real;

var

dx,dy,x,sum:real;

i:word;

function simple_simpson(x:real):real;

var

y,sum:real;

j,v:word;

begin

sum:=F(x,y0)+F(x,y1);

y:=y0;

for j:=1 to 2*y_div-1 do

begin

y:=y+dy;

if Odd(j) then

sum:=sum+4.0*F(x,y)

else

sum:=sum+2.0*F(x,y);

end;

simple_simpson:=sum;

end;{simple_simpson}

begin{doudle_simpson}

dx:=(x1-x0)/(2.0*x_div);

dy:=(y1-y0)/(2.0*y_div);

x:=x0;

sum:=simple_simpson(x0)+simple_simpon(x1);

for i:=1 to 2*x_div-1 do

begin

x:=x+dx;

if Odd(i) then

sum:=sum+4.0*simple_simpson(x)

else

sum:=sum+2.0*simple_simpson(x);

end;

double_simpson:=dx+dy*sum/9.0;

end;{double_simpson}

Недостатком рассмотренных функций интегрирования является то, что они не дают возможности явно задать точность вычисления интеграла. Точность связана с количеством точек разбиения, но ее значение в этих функциях не определяется с адаптивным выбором шага разбиения. Такой функцией является adaptive_simpson. Параметры eps и eta задают соответственно абсолютную и относительную погрешности. Их роль поясняется следующим неравенством:

Функция adaptive_simpson (листинг1.4) использует рекурсивную процедуру simpson3point, которая вычисляет значение интеграла по интервалу [x₀, x₀+δx], где x₀ – не обязательно исходная левая граничная точка.

Если трехточечный метод Симпсона не дает достаточную точность на данном интервале, этот интервал делится на три равные части, и метод вновь применяется к каждой из полученных частей. В результате получим 7 точек разбиения, но вычислять функцию F(x) придется только в четырех из них, поскольку значения в других трех точках уже известны.

При адаптивном разбиении имеется одна тонкость. При переходе к подынтервалам, составляющим одну треть от исходного, чтобы получить новые абсолютную и относительную погрешности, надо поделить eps и eta на .

Листинг 1.4. Функция adaptive_simpson модуля integral

Function adaptive_simpson(F:real_fun:x0,x1,eps,eta:real):real;

const

max_level=35;

var

k,nest_level:word;

integral_abs:real;

function simpson3poin(x0,delta_x, estimate, integral_abs,

eps,eta,left,middle,right:real):real;

var

dx3,sum,eps3,eta3,factor,left_integ,

middle_integ, right_integ,F1,F2,F4,F5:real;

begin

Inc(nest_level);

dx3:=delta_x/3.0;

F1:=F(x0+0.5*dx3);

F2:=F(x0+dx3);

F:=F(x0+2.0*dx3);

F5:=F(x0+2.5*dx3);

Inc(no_evaluations,4);

factor:=dx3/6.0;

left_integ:=factor*(left+4.0*F1+F2);

middle_integ:=factor*(F2+4.0*middle+F4);

right_integ:=factor*(F4+4.0*F5+right);

sum:=left_integ+middle_integ+right_integ;

integral_abs:=integral_abs- Abs(eastimate)+Abs(left_integ)+Abs(middle_integ)+Abs(right_integ);

if (nest_level>1) and ((nest_level=max_level) or

(Abs(sum- estimate)<=eps+eta*integral_abs)) then simpson3point:=sum

else

Begin

If nest_level>highest_level then

Inc(highest_level);

Eps3:=0.577*eps;

Eta3:=0.577*eta;

Left_integ:=simpson3point(x0,dx3,left_integ,integral_abs,eps3,eta3, left,F1,F2);

middle_integ:=simpson3point(x0+dx3,dx3,middle_integ, integral_abs,eps3,eta3,F2,middle,F4);

right_integ:=simpson3point(x0+2.0*dx3,dx3,right_integ,integral_abs,eps3,eta3,F4,F5,right);

simpson3poin:=left_integ+middle_integ+right_integ;

end;

Dec(nest_level);

End; {simpson3point}

Begin {adaptive_simpson}

nest_level:=1;

highest_level:=1;

no_evaluations:=3;

adaptive_simpson:=simpson3point(x0,x1-x0,0.0,0.0,eps,eta,F(x0),F(x0+0.5*(x1-x0)),F(x1));

end;{adaptive_simpson}

1.3. Метод Ромберга и его реализация на языке Pascal.

Интегрирование следующим методом – методом Ромберга – основано на правиле трапеций, использующем кусочно-линейное приближение для интегрируемой функции. Основной факт относительно погрешности в методе трапеций следующий.

Теорема. Пусть F(x) – гладкая функция на интервале [a,b], и этот интервал делится на т равных частей, каждая длиной h =. Пусть I(h) обозначает соответствующее приближение метода трапеций:

I(h) =

где f_i=F(a+jh) – значение интегрируемой функции в точке a+jh.

Тогда

Где a_k – некоторая константа.

Основное здесь то, что погрешность в методе трапеций может быть выражена рядом по четным степеням шага интегрирования h. Построим таблицу значений T_ik:

В нулевой строке T_0k = I((b – a)/2k), так что T₀₀,T₀₁,… являются последовательными приближениями метода трапеций для интеграла, каждое с удвоенным по сравнению с предыдущим числом интервалов. Согласно приведенной выше теореме,

где h = ((b – a)/2k.

Отсюда следует, что

Поэтому положим

В общем случае строим j-ю строку таблицы Ромберга по формуле

а оценка погрешности имеет вид

где h = (b – a)/2k.

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

F₀ + 2F₁ + 2F₂ + …+ 2F_2n-1 + F_2n =

= (F₀ + 2F₂ + 2F₄ + …+ 2F_2n-2 + F_2n) + 2(F₁ + F₃ +…+F_2n-1).

Таким образом, чтобы модифицировать предыдущее приближение, необходимо вычислить сумму значений функции в новых средних точках. Это делается в цикле со счетчиком k. Метод Ромберга реализован в функции romberg (листинг1.5).

Листинг 1.5. Функция romberg модуля integral

Function romberg (F:real_fun;x0,x1,eps,eta:real; min,max:word):real;

const

abs_max=30;

var

p,dx,error,F_of_x0, F_of_x1, F_of_xk,

roundoff_error,integral_abs,tolerance,

previous_estimate,current_estimate,

mid_sum, temp_sum, mid_sum_abs:real;

table:array[0..abs_max] of real;

j,n:word;

k,r:longint;

done:Boolean;

denom:array[1..abs_max] of real;

begin

p:=1.0;

for k:=1 to abs_max do

begin

p:=4.0*p;

denom[k]:=1.0/(p-1.0);

end;

dx:=x1-x0;

F_of_x0:=F(x0);

F_of_x1:=F(x1);

current_estimate:=0.0;

previos_estimate:=0.0;

done:=False;

table[0]:=0.5*dx*(F_of_x0+F_of_x1);

integral_abs:=0.5*Abs(dx)*(Abs(F_of_x0)+Abs(F_of_x1));

n:=1;

r:=1;

repeat

dx:=0.5*dx;

mid_sum:=0.0;

mid_sum_abs:=0.0;

roundoff_error:=0.0;

for k:=1 to r do

begin

F_of_xk:=F(x0+(2*k-1)*dx);

mid_sum_abs:=mid_sum_abs+Abs(F_of_xk);

F_of_xk:=F_of_xk+roundof_eroor;

temp_sum:=mid_sum+F_of_xk;

roundof_error:=(mid_sum-temp_sum)+F_of_xk;

mid_sum:=temp_sum;

if KeyPressed then

Halt;

end;

table[n]:=0.5*table[n-1]+dx*mid_sum;

integral_abs:=0.5*integral_abs+Abs(dx)*mid_sum_abs;

for j:=n-1 downto 0 do

table[j]:=table[j+1]+denom[n-j]*(table[j+1]-table[j]);

if n>=min then

begin

tolerance:=eta*integral_abs+eps;

error:=Abs(table[0]-current_estimate)+Abs(current_estimate-previos_estimate);

done:=(error<tolerance);

end;

Inc(n);

done:=done or(n>max);

previous_estimate:=current_estimate;

current_estimate:=table[0];

r:=r+r;

until done;

romberg:=current_estimate;

end;

1.4. Метод Гаусса и его реализация на языке Pascal.

Теперь перейдем к гауссовским квадратурам – семейству правил интегрирования, основанных на неравномерном разбиении основного интервала интегрирования. Вообще, метод гаусса с n точками точен для полиномов степени 2n – 1. В функции gauss3 (листинг1.6.) основной трехточечный алгоритм Гаусса применяется к каждой из n равных частей интервала. Для интервала [-1,1] узлами квадратурной формулы являются нули полинома Лежандра третьей степени P₃ = (5x³ – 3x)/2, а коэффициенты выбираются специальным образом.

Листинг 1.6. Функция gauss3 модуля integral

Function gaus3(F:real_fun; x0,x1:real; n:word):real;

var

t,sum,x,z,dx:real;

i,k:word;

gzero,gweight:array[1..3] of real;

procedure initialize_constants;

var

s,t:real;

j:word;

begin

gzero[1]:=-sqrt(0.6);

gzero[2]:=0.0;

gzero[3]:=sqrt(0.6);

gweight[1]:=5.0/9.0;

gweight[2]:=8.0/9.0

gweight[3]:=5.0/9.0;

for j:=1 to 3 do

begin

gzero[j]:=0.5*(1.0+gzero[j]);

gweight[j]:=0.5*gweight[j]);

end;

end; {initialize_constants}

begin {gauss3}

initialize_constants;

dx:=(x1-x0)/n;

x:=x0;

sum:=0.0;

for i:=0 to n-1 do

begin

t:=0.0;

for k:=1 to 3 do

begin

z:=x+dx*gzero[k];

t:=t+gweight[k]*F(z);

end;

sum:=sum+dx*t;

x:=x+dx;

end;

gauss3:=sum;

end;{gauss3}

Дадим некоторый обзор некоторых свойств полиномов Лежандра. Рекурсивное определение полиномов Лежандра приводилось ранее в этой работе. Они образуют ортогональное (но не ортонормированное)семейство на промежутке [-1,1], то есть

, m≠n,

Величина второго интеграла определяет нормировку для этих полиномов. Имеет место также следующее представление полиномов Лежандра:

Другая явная формула:

Приведем несколько первых полиномов Лежандра:

P₀(x) = 1,

P₁(x) = x,

Очевидно, что в общем случае полиномы Лежандра нечетной степени являются нечетными функциями, а четной степени – четными функциями.

Нам требуется найти нули полинома P_n(x). Важно здесь то, что эти нули являются простыми и принадлежат открытому интервалу (-1,1). Таким образом,

-1<x₁<x₂<…<x_n<1, P_n(x_j) = 0.

Соответствующая формула гаусовского интегрировании (с остатком) имеет следующий вид:

В этой формуле

Где -1 <<1.

Веса задаются несколькими эквивалентными формулами

Процедура compute_gauss_coeffs (листинг1.7). предназначена для вычисления нулей и весов квадратурной формулы Гаусса. Подпрограмма legendre_poly вычисляет значения P_n(x), P_n-1(x) и P_n’(x). Последнее получается дифференцированием основной рекуррентной формулы для P_n(x):

Нули находятся предварительным делением интервала и применением метода секущих, после чего следует ньютоновские итерации, в которых используются значения производных. Затем применяется первая формула для весов, в которой вновь используются значения производных. Здесь zero – массив нулей полинома Лежандра n-й степени, а weight – массив соответствующих весов.

Приближенное вычисление интеграла 📙 Реферат → 🆔 216299