Mathcad - ЛР5

Сравним метод центральных прямоугольников с методом Симпсона при равном числе узлов:

n := 2 , 4 .. 20

E_Centr := 0

integraln :=

IntСentr(f , a, b, n)

E_Centræ n

ö :=

true − integraln

ç

−1÷

è 2

ø

max_P(E_Centr , 9) = 0.076141573677683

min(E_Centr) = 5.9341029956772 × 10− 4

max_P(E_Simpson , 9) = 0.627482713901151

min(E_Simpson) = 1.24684471785041 × 10− 5

n2 := 1 .. 9

0.015

0.01

E_Centrn2

E_Simpsonn2

0.005

0

0

5

10

n2

Вывод: В ходе проведенного вычиcлительного эксперимента метод Симпсона показал более высокую скорость сходимости.

Задача 5.3. Убедиться, что квадратурные формулы Гаусса с N (N=1,2,3,4) узлами точны для всех многочленов степени m ≤ 2N −1.

Порядок решения задачи:

1)Составить программу (можно назвать, к примеру, IntGauss( f , a,b, N ) ),

вычисляющую определенный интеграл произвольной функции на произвольном отрезке [a,b] на основе квадратурной формулы Гаусса с заданным количеством узлов N (N=1,2,3,4).

2)	Составить несколько многочленов	Pm ( x)		= c0 + c1 x + c 2 x 2 + ... + c m x m для
	различных степеней m = 0..7 .
3)	Аналитически (по формуле Ньютона-Лейбница) найти интегралы для всех
	многочленов Im = òb Pm (x)dx = (c0 x + c1	x2	+ c2	x3	+ ... + cm	xm+1	) |ba .
		2		3
	a					m +1

4)Найти значение тех же интегралов на основе квадратурной формулы Гаусса

IGm,N = IntGauss(Pm , a,b, N ) , m = 0..7 , N=1,2,3,4.

5)		Для каждого фиксированного N (N=1,2,3,4) вывести на один чертеж значенияIm и
		IGm,N .			Убедиться, что c точностью до погрешности округления Im = IGm,N																		для
		любых m ≤ 2N −1.
6)		Сделать выводы.
Решение задачи:
æ	2	1	5		0.3478548451	ö	æ0		-0.5773502692	-0.7745966692				-0.8611363116					ö
ç			9			÷
							ç	0	0.5773502692	0				-0.3399810436					÷
ç			8			÷	ç												÷
ç	0	1			0.6521451549	÷	t := ç	0	0	0.7745966692				0.3399810436					÷
			9
A := ç						÷	ç		0	0				0.8611363116					÷
ç	0	0	5		0.6521451549	÷	è0												ø
			9
ç						÷
ç						÷
è0		0	0		0.3478548451 ø		IntGauss(f , a, b, N) :=				h ¬	b - a
											2
											S ¬ 0
											for i Î 0 .. N - 1
															æ a + b					b - a			ö
											S ¬ S + Ai , N−1×f				è			+				×ti, N−1
																2			2
																							ø
											S ¬ h×S
											S

Составим 8 многочленов с произвольными коэффициентами: a := 0		b := 1
	x := a, a + 0.01 .. b
P0	(x) := 5×x0
P1	(x) := 11 + x

P2(x) := 5 + 9×x - x2

P3(x) := 8 + 2×x + 9×x2 - 7×x3

P4(x) := x + 4×x2 + 5×x3 + x4

P5(x) := 31 - 2×x + 7×x2 - x3 + 12×x4 + 5×x5

P6(x) := 32×x3 + 6×x5 - 3×x6

P7(x) := 10 + 8×x + 4×x2 + 5×x3 + 7×x6 + 8×x7

Вычислим значение представленных многочленов аналитически:

IP0 := 5×b - 5×a

æ

2

ö

æ

2

ö

IP1

:= ç

11×b +

b

÷

-

ç11×a +

a

÷

è

2 ø

è

2 ø

æ

2

3

ö

æ

2

3

ö

IP2

:= ç

5×b + 9×

b

-

b

÷

- ç5×a + 9×

a

-

a

÷

è

2

3 ø

è

2

3 ø

æ

2

3

4

ö

æ

2

3

4

ö

IP3

:= ç

8×b + 2×

b

+ 9×

b

- 7×

b

÷

-

ç

8×a + 2×

a

+ 9×

a

- 7×

a

÷

2

è

2

3

4 ø

è

3

4 ø

æ

2

3

4

5 ö æ

2

3

4

5

ö

IP4

:= ç

b

+ 4×

b

+ 5×

b

+

b

÷

-

ç

a

+ 4×

a

+ 5×

a

+

a

÷

è 2

3

4

5 ø

è 2

3

4

5 ø

æ

2

3

4

5

6

ö

æ

2

3

IP5

:= ç

31×b - 2×

b

+ 7×

b

-

b

+ 12×

b

+ 5×

b

÷

- ç

31×a - 2×

a

+ 7×

a

4

è

2

3

5

6 ø

è

2

3

æ

4

6

7

ö

æ

4

6

7

ö

IP6

:= ç

32×

b

+ 6×

b

- 3×

b

÷

-

ç

32×

a

+ 6×

a

- 3×

a

÷

4

4

è

6

7 ø

è

6

7 ø

æ

2

3

4

7

8

ö

æ

2

3

IP7

:= ç

10×b + 8×

b

+ 4×

b

+ 5×

b

+ 7×

b

+ 8×

b

÷ -

ç

10×a + 8×

a

+ 4×

a

7

è

2

3

4

8 ø

è

2

3

	4	5		6		ö
-	a	+ 12×	a	+ 5×	a	÷
	4
		5		6 ø

4		7		8		ö
+ 5×	a	+ 7×	a	+ 8×	a	÷
4		7		8 ø

Произведем расчет значений интегралов численным методом Гаусса при различном количестве узлов N для каждого из представленных выше многочленов степеней m:

N := 1 .. 4

m := 0

IGm , N−1 := IntGauss(P0 , a, b , N)

m := 1

IGm , N−1 := IntGauss(P1 , a, b , N)

	m := 2
	IGm , N−1 := IntGauss(P2 , a, b , N)
	m := 3
	IGm , N−1 := IntGauss(P3 , a, b , N)		Значения интегралов, полученные методом Гаусса для
	m := 4
			различных значений количества узлов находятся в
	IGm , N−1 := IntGauss(P4 , a, b , N)		строках нижепредставленной таблицы.
	m := 5		æ	5	5	5	5
	IGm , N−1 := IntGauss(P5 , a, b , N)
			ç	11.5	11.5	11.5	11.5
	m := 6		ç
			ç	9.25	9.16666666666367	9.16666666667559	9.1666666666699
	IG	:= IntGauss(P6 , a, b , N)	IG = ç	10.375	10.2499999999955	10.2500000000134	10.250000000004
	m , N−1		ç	2.1875	3.27777777781721	3.28333333321462	3.2833333332897
	m := 7		ç
			ç	32.53125	35.1805555556569	35.3166666663355	35.316666666548
	IGm , N−1 := IntGauss(P7 , a, b , N)
			ç
			è4.140625		8.55555555571503	8.57249999952149	8.5714285712524

Произведем сравнение точности численного метода Гаусса при фиксированном количестве узлов со значениями, вычисленными аналитически.

m := 0 .. 7

N := 1

40

IPm

20

IGm , N−1

0
	0	5										10
										m

N := 2

40
IPm
20
IGm , N−1
0	0	5	10
N := 3		m
40

IPm

20

IGm , N−1

0
	0	5										10
										m

N := 4

40

IPm

20

IGm , N−1

0
	0	5										10
										m

Произведем сравнение точности вычислений численного метода Гаусса для различного количества узлов:

m := 0 .. 7 N := 1 .. 4

Em , N−1 := IPm − IGm , N−1

	æ	0	0		0			0
	ç	0	0		0			0
	ç
	ç		2.99493763122882 ´	− 12		´	− 12	3.29158922340866 ´	− 1
	ç	0.083333333333334		10	8.92619311798626		10		10
	ç	0.125	4.49240644684323 ´	− 12	1.33884014985597	´	− 11	4.9382720135327 ´	− 12
	ç			10			10		10
E =			5.55555551612397 ´ 10− 3			´ 10− 10		4.36197744591027 ´ 10− 1
	ç	1.09583333333333			1.18712151220279
	ç					´ 10− 10		1.18234311230481 ´ 10− 1
	ç 2.78541666666667		0.136111111009789		3.31148442000995
	ç				1.07142809292071 ´ 10− 3			1.7610979341498 ´ 10− 10
	ç 4.43080357142857		0.015873015713542
	ç							1.3245227137304 ´ 10− 10
	è 2.78645833333334		0.398148148045578		0.012500000393182

Известно, что квадратурная формула Гаусса дает точность вычисления интеграла для многочленов вплоть до степени m, сравнимую с точностью

округления, если m<=2*N-1, где N - количество узлов. Для проверки данного утверждения выберем наибольшую возможную степень многочлена для каждого из N:

N := 1 .. 4

mN := 2×N - 1

E(mN), N−1 =

0

4.49240644684323·10 -12

3.31148442000995·10 -10

1.3245227137304·10 -10

Вывод: В ходе проведенного вычиcлительного эксперимента квадратурная формула Гаусса выдавала результат, совпадающий с аналитически вычисленным

значением интеграла с точностью до погрешности округления для всех многочленов степени до m включительно, если выполнялось следующее неравенство: m<=2*N-1, где N - количество узлов.

Задача 5.4. Вычислить значение интеграла òb f (x)dx из задачи 5.1, используя квадратурную

a

формулу Гаусса с одним, двумя, тремя, четырьмя узлами. Определить абсолютную погрешность результата. Построить график зависимости погрешности от числа узлов. Сделать выводы, сравнивая с результатами, полученными в задачах 5.1 и 5.2.

Решение задачи:

Задана функция

f(x) :=		1		на отрезке [2;2]	a := −2	b := 2
		9 − x2

true := asin

æ 2

ö

− asin

æ

−2 ö

= (1.45945531245393)

è 3

ø

è

3 ø

N := 1 .. 4

integral := 0

integralN−1 :=

IntGauss(f , a, b, N)

E_GaussN−1 :=

true − integralN−1

æ

0.126121979120599

ö

æ1.33333333333333

ö

ç

÷

ç

0.014825075420188

÷

ç

÷

ç1.44463023703374

÷

E_Gauss = ç

÷

integral =

ç1.45759141942188

÷

ç1.86389303204892 × 10− 3

÷

ç

− 4

÷

ç

÷

è1.45921307015357

ø

è 2.4224230036074 × 10

ø

E_GaussN−1

0.2

0.1

0

1

2

3

4

N

0

0

E_LeftN−1 =

0

0.329399069545899

E_SimpsonN−1 =

0

0.622003164434844

1

0.10163854521265

1

3.7183483503247·10 -3

2

0.048514093223115

2

1.0160368771579·10 -3

3

0.028198397565906

3

3.7843268755644·10 -4

		0			0
E_RightN−1 =	0	0.329399069545899	E_CentrN−1 =	0	0.045241750080838
	1	0.10163854521265		1	0.013531549751618
	2	0.048514093223115		2	6.31471233809422·10 -3
	3	0.028198397565906		3	3.62258247620195·10 -3

			æ	0.126121979120599	ö
		0
			ç	0.014825075420188	÷
	0	3.90710783354556	ç		÷
E_TrapeciaN−1 =	1	0.10163854521265	E_Gauss = ç	− 3	÷
			ç	1.86389303204892 × 10	÷
	2	0.048514093223115
			ç	− 4	÷
	3	0.028198397565906
			è	2.4224230036074 × 10	ø

Вывод: В ходе проведенного вычиcлительного эксперимента квадратурная формула

Гаусса оказалась ощутимо точнее формулы Симпсона и всех остальных методов при одном и том же количестве узлов.

Задача 5.5.	Функция y = f (x) задана на отрезке [a,b] . Требуется вычислить значение
интеграла òb	f (x)dx с заданной точностью ε =10−6 при помощи двух методов:
a

А) трапеций; В) Симпсона.

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

Теоретический материал:

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

A)

Для трапеций нужно решить неравенство:

M

2 (b - a)3

£ ε Þ n ³ ceil(

(b - a)3

M 2

12n2

12ε

B)

А для Симпсона:

M

4 (b - a)5

£ ε откуда число узлов n ³ ceil(4

(b - a)5 M 4

)

180n4

180ε

Решение задачи:

ε := 10− 6

a := 3

b := 6

f(x) :=

atan(x)

x − 2

A)

		d2
	df2(x) :=			f(x)	absdf2(x) :=	df2(x)
		2
		dx			1

M2 := absdf2(3) n := 0 .. 9

df2(z) 0.5

M2 = 0.77678432929884

0
	3	4		5							6
				z

	é					ù
			(b - a)	3
	n := ceilê				×M2ú			n = 1323
			12×e
	ë					û
							b - a
	Trapeсia(f , a, b, n) :=					h ¬
								n
						for	i Î 0 .. n
						xi ¬ a + h×i
						S ¬		f(x0) + f(xn)
							2
						for	i Î 1 .. n - 1
						S ¬ S + f(xi)
						S ¬ h×S
						S

integral_T := Trapeсia(f , a, b, n) integral_T = 2.67094

B)

d2

df4(x) :=

f(x)

absdf4(x) :=

df2(x)

2

dx

1

M4 := absdf2(3)

n := 0 .. 9

M4 = 0.77678432929884

absdf4(z) 0.5

0

3

4

5

6

z

é4

ù

(b - a)5

ê

ú

Количество узлов должно быть четным

n := ceilê

×M4ú

n = 33

180×e

ë

û

n := 34

Увеличение количества узлов не ухудшает точности

b - a

Simpson(f , a , b, n) :=

h ¬

n

for

i Î 0 .. n

xi ¬ a + h×i

S ¬ f(x0) + f(xn)

for

i Î 1 ..

æ nö

è 2 ø

S ¬ S + 4×f(x2×i-1)

for

i Î 1 ..

n - 1

2

S ¬ S + 2×f(x2×i)

S ¬

h

×S

3

S

integral_S :=

Simpson(f , a, b, n)

integral_S = 2.67094

Вычислим теперь значение интеграла встроенными средствами MathCAD:

ób

integral_M := ô f(x) dx integral_M = 2.67094

õa