Формула Симпсона. Алгоритм, блок-схема. Геометрическая иллюстрация. Оценка погрешности.

 

b

n

Найдем коэффициенты формулы ò f (x) = (b - a)å yi H i,

a

i =0

где H i

=

1 n

(-1) n - i t (t -1)...(t - n)

dt

, i=0,1,…,n при n= 2.

n ò

i! (n - i)! (t

- i)

При i = 0

2 t (t -1)(t -2)

æ t 3

3 t 2

ö

æ

ö

H

=

dt

=

(t 2

- 3 t + 2) dt =

ç

-

+ t ÷

=

- 6

+ 4

=

×

=

ò

ò

ç

÷

× 2 × t

ç

÷

4 3 6

è

ø

4 è 3

ø

При i = 1

2 (-1) t (t -1)(t - 2)

æ t 3

ö

1 æ 8

ö

H

=

dt

= -

(t

- 2 t) dt = -

ç

- t

÷

= -

- 4

=

×

=

ò

ò

ç

÷

1×1× (t -1)

ç

÷

2 3 3

è

ø

2 è 3

ø

При i = 2

2 t (t -1)(t -2)

æ t 3

t 2

ö

æ

ö

H

=

dt

=

(t 2

- t) dt =

ç

-

÷

=

- 2

=

×

=

ò

ò

ç

÷

1× 2

× (t - 2)

ç

3 2

÷

4 3 6

è

ø

4 è 3

ø

b

n

Формула ò f (x) = (b - a)å yi H i на отрезке [ x0, x2 ] примет вид:

a

i =0

x 2

n

æ 1

ö

ò f (x) dx»(x 2- x 0)×å H i yi =2 h ç

y 0+

y 1

+

y 2

÷

x 0

i =0

è 6

ø

x 2		2 h æ							ö
ò	f (x) dx»		+ 2 y 1 +
	ç		y 0		y		÷
x 0		è					ø

При n=2m применив формулу к каждой паре частичных отрезков [ x_2i-2, x_2i ] (i=1,2,…,m)

получим формулу Симпсона:

b

2 h æ y + y

ö

ò f (x) dx»

ç

2 m

+ 2 y 1 + y 2

+... + 2 y 2 m -1 ÷ =

a

è

ø

=

2 h æ y

- y

2 m

+ (2 y

+ y

) + (2 y

+ y

) +... + (2 y

+ y

ö

ç

2 m -1

) ÷

è

2 m

ø

b

b

Рассмотрим погрешность: ò f (x) dx =ò Ln (x) dx + Rn (x)

a

a

На отрезке [ a, b ] R =

h 4(b - a)

f IV (x), x Î[ a, b ]

n

или

Rn

£ M

b - a

h 4

, где M = max

f IV (x)

x Î[ a, b ]

При

вычислении

по

методу

повторного

счета можно использовать формулу:

R 2 n

£

I n

- I 2 n

. Если при вычислении интеграла требуемая точность не достигнута (т.е.

 

I _n - I _{2 n} >15 e),предусматривается повторный счет с шагом,уменьшенным вдвое.

 

Программа вычисления по формуле Симпсона методом повторного счета:

program lab4_2;

var n: integer;

S,a,b,e,h,x,I_n,I_n2,M: real;

function f(x: real):real;

 

begin {записать, функцию в виде f:=[математическое выражение]} f:=sin(x); end; begin

 

write('Введите концы отрезка интегрирования: '); readln(a,b);

write('Введите погрешность e: '); readln(e);

I_n:=0;

n:=4;

Repeat

 

h:=abs(b-a)/n;

s:=(f(a)-f(b))/2;

x:=a+h;

repeat s:=s+2*f(x)+f(x+h);

x:=x+2*h;

until x>=b;

 

I_n2:=2*h*s/3;

n:=n*2;

M:=abs(I_n - I_n2);

I_n:=I_n2;

Until M<=15*e;

writeln('Интеграл I=',I_n2:12:7);

readln;

 

end.

 

Решение обыкновенных дифференциальных уравнений. Численное интегрирование уравнений порядка выше, чем первый. Решение систем дифференциальных уравнений.

Дифференциальное уравнение 1-го порядка, разрешенное относительно производной, имеет вид:

 

y ¢ = f (x, y)

 

Решением

 

 

Дифференциального

 

 

Уравнения

 

(1)

 

(1) называется функция

 

 

y(x),

 

подстановка которой в уравнение обращает его в тождество: y ¢(x) = f (x, y (x)).

 

График решения y=y(x) называется интегральной кривой.

 

Задача Коши для дифференциального уравнения(1)состоит в том,чтобы найтирешение дифференциального уравнения (1), удовлетворяющее начальному условию y(x₀)=y₀ (2).Пару чисел (x₀,y₀) называют начальными данными.

 

Решение задачи Коши называется частным решением дифференциальногоуравнения (1)при условии(2).

 

Геометрически задача Коши означает, что требуется найти интегральную кривую y=y(x),проходящую через заданную точку (x₀,y₀).

Теорема о существовании и единственности решения задачи Коши.

 

Пусть функция f(x,y) – правая часть уравнения y ¢ = f (x, y) - непрерывна вместе со

 

своей частной производной по переменной y	¶ f (x, y)	в некоторой области D на
dy

плоскости. Тогда при любых начальных данных (x₀,y₀)ÎD задача Коши имеет единственное решение y=y(x).

 

При выполнении условий теоремы через точку (x₀,y₀) на плоскости проходит единственная интегральная кривая.

 

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

 

· аналитические методы (решение в виде аналитического выражения);

· графические методы (решение в виде графика);

· численные методы (решение в виде таблицы).

 

Численное решение задачи Коши состоит в том, чтобы получить искомое решение y(x) в виде таблицы его приближенных значений аргумента x на некотором отрезке [ a, b ]:

x0=a, x1, x2, …, xm=b

(3)

Точки (3) называют узловыми, множество этих точек называют сеткой наотрезке [ a, b ].

 

Как правило, используют равномерную сетку с шагом h: h = ^b _m ^{- a} x_i=x₀+ih (i=0, 1, …, m)

Приближенные значения численного решения задачи Коши в узловых точках обозначим y_i.

y_i» y(x_i),где(i=0, 1, …, m)

Начальное условие выполняется точно: y₀ = y(x₀).

 

Величина погрешности численного решения задачи Коши на сетке отрезка [ a, b ] оценивается величиной d = max{ yi - y (xi)},

1£ i £ m

 

т.е. расстоянием между векторами приближенного решения (y₀, y₁, …,y_m) и точного решения (y(x₀), y(x₁), …,y(x_m)) на сетке по m -норме.

{ y’=f(x, y); (1)

y(x₀)=y₀;

 

Решить обыкновенное диф. ур. – найти функ y=y(x), которая при подстановке даёт тождество.

 

Существование и единственность решения ур. (1) обеспечиваются теоремой. Теорема Пикара: если функ. f определена и непрерывна в некоторой области G,

 

определяемой неравенствами:

|x–x₀|≤a, |y–y₀|≤b,

и удовлетворяет в этой области условию Липшица по y:

|f(x, y₁)–f(x, y₂)|≤M|y₁–y₂|,

 

то на некотором отрезке |x–x₀|≤h, где h – положительное число, существует, и при том только одно, решение y=y(x) уравнения (1), удовлетворяющее начальному условию

y₀=y(x₀).

M=max|f’_y(x, y)| –константа Липшица.

Для нормальных систем ОДУ используются те же методы,только в векторной форме:

y’=f(x, y),где y=(y₁(x),...,y_n(x)); f(x, y)=(f₁,...,f_n);

Уравнения более высокого порядка принято сводить к системе:

 

Пусть:

y(n)=F(x, y(n-1),..., y);

 

Обозначим:

z₁=y;

z₂=y’;

………

z_n=y^(n-1);

Вместо ур. запишем систему с переменной z:

|z’₁=z₂;

|z’₂=z₃;

{………

|z’_n-1=z_n;

|z’_n=F(x, z_n,…, z₁);

нормальная сист. диф. ур.

}

M1(x1,y1) равен x2=x1+h