Практикум j практическое примщенше численных методов

Е


^ш(х)
(х - xjw’ixi)


f(Xi),

одынтегральная функция f(x) в (10.1) заменяется интерполяционным мно­гочленом Лагранжа
где
и(х) = П(^ж “ ^х*)> ^ш'(^х) = Т^^х)-
г=0 ^йХ
230 Численное интегрирование
Д
Ci =


= / в{*
J а


ш{х)


dx, г = 0, l,...,n.


(10.G)

ля коэффициентов квадратурной формулы (10.2) получим представление
Квадратурные формулы Гаусса
Для повышения точности квадратурной формулы можно оптимизировать выбор не только коэффициентов квадратурной формулы с*, г — 0,1,..., п, по и узлов интерполяции ж», г = 0,1,... , п. Квадратурные формулы интерполя­ционного типа (10.2), (10.6) являются точными для алгебраических полино­мов степени п. За счет выбора узлов интерполирования строятся квадратур­ные формулы наивысшей алгебраической степени точности (квадратурные формулы Гаусса), которые точны для любого алгебраического многочлена степени 2п + 1.
Потребуем, чтобы квадратурная формула (10.1) была точна для любого ал­гебраического многочлена степени т. Это означат¹, что формула точна для функций f(x) = x^OL₎ а = 0,1,..., т:
/ь ^п
g(x)x^adx = ^2 CiX“, г = 0,1,..., т. (10.7)
г=0
Для определения 2п + 2 неизвестных с*, ж», г = 0,1,... , п имеем нелинейную систему (10.7) m + 1 уравнений.
Для знакопостоянной весовой функции д{х) система уравнений (10.7) при т = 2п+1 имеет единственное решение. При этом квадратурная формула яв­ляется интерполяционной, т.е. коэффициенты вычисляются согласно (10.6), а узлы должны быть такими, чтобы многочлен
Ч*) = Р(ж - Xi)
г=0
был ортогонален с весом д(х) любому многочлену q(x) степени меньше n + 1:
/ g(x)uj(x)q(x)dx = 0.
J а

3. 
   Упражнения
   

Упражнение 10.1 Напишите программу для приблиэюсппого вычисления интеграла от функции f(x) на интервале [а, Ь] с использованием квадро' туриой формулы трапеций па основе рекуррентного уточнения при увели- чспии числа частичных отрезков в два раза. С помощью этой программы
рассчитайте значение интеграла ошибок





при х = 1 с различной точностью.
Пусть имеется 2^п частичных отрезков длиной h = (b-a)/2ⁿ. При п — 0 имеем
/о = (/(а) + /(Ь))£,
при п = 1 —
h = (/(а) + №)± + /(а + Л))А = + /(а + А))А.
В общем случае имеет место рекуррентное соотношение
2^п ~ 1
/»+1 ⁼ 7.I71 + ^ f{p> + (2i — l)h)h.
¹ г=1
Тем самым при увеличении числа частичных отрезков вычисления проводят­ся только в новых точках.
В модуле trap функция trap О вычисляет приближенное значение интеграла по рекуррентной квадратурной формуле трапеций.
def trap(f, а, Ъ, tol=l.e-6):
и и и
Integral of f(x) from a to b computed by trapezoidal rule
и и и
h = b - a
Iold = h * (f(a) + f(b)) / 2. m = 1
* kMax = 25
for k in ranged, kMax): x = a + h / 2. sum = 0.
for i in range(m): sum = sum + f(x) x = x + h
Inew = (Iold + h*sum) /2
if (k > 1) and (abs(Inew - Iold) < tol): break Iold = Inew m = m * 2 h = h / 2. return Inew
Решение задачи приближенного вычисления интеграла ошибок дается следу, ющей программой.
import math as mt from trap import trap def f(x):
return 2.*mt.exp(-x**2) / mt.sqrtCmt.pi) a = 0. b = 1.
for n in range(5, 11):
er = 10.**(-n)
Inew = trap(f, a, b, tol=er)
print 4ol = ’, er, ’Integral = ’, Inew
§В|1Щ|111в^
Точность вычислений контролируется параметром tol, при его уменьшении точность вычислений повышается (erf(ж) = 0.8427007929497...).
У
пражнение 10.2 Напишите программу для приблиэ1сепиого вычисления интеграла от функции f(x) па интервале [а, 6] с использованием квадра­турной формулы Гаусса с весом д(х) па основе табличного задания узлов и коэффициентов. Используйте эту программу для вычисления интеграла
при различном числе узлов и сравните приблиэюепиое значение интеграла с точным.
Узлы и коэффициенты для рассматриваемой квадратурной формулы Гаусса— Лежандра даны в табл. 10.1.
В модуле gauss функция gauss () вычисляет приближенное значение интегра­ла по этим табличным данным.
|||||||||||||||||(|||||||||^
import numpy as np def gauss(f, a, b, n):
I
и к и

ntegral of f(x) from a to b computed by Gauss-Legendre quadrature using m nodes.

Таблица 10.1 Квадратурная формула Гаусса—Лeoicaudpa

n	узлы Xi	коэффициенты c,
2	± 0.57735027	1.
3	± 0.77459667	0.5555.5556
	0.	0.88888889
4	± 0.86113631	0.34785485
	± 0.33998104	0.65214515
5	± 0.90617985	0.23692689
	± 0.53846931	0.47862867
	0.	0.56888889
6	± 0.93246951	0.17132449
	± 0.66120939	0.36076157
	± 0.23861919	0.46791393
7	± 0.94910791	0.12948497
	± 0.74153119	0.27970539
	± 0.40584515	0.38183005
	0.	0.41795918
8	± 0.96028986	0.10122854
	± 0.79666648	0.22238103
	± 0.52553241	0.31370665
	± 0.18343464	0.36268378

и и и
if n > 8 or n < 2:
print ’The number of nodes must be greater than 2 and less than 8’ return 0
x = np.zeros(n, ’float’) c = np.zeros(n, ’float’) if n == 2:
x[0] = 0.57735027 x[l] = - x[0] c[0] = 1. c[l] = c [0] if n == 3:
x[0] = 0.77459667 x[i] = - x[0] x[2] = 0. c[0] = 0.55555556 c [1] = c [0] c[2] = 0.88888889 if n == 4:
x[0] = 0.86113631 x[l] = - x[0] x[2] = 0.33998104 x[3] = - x [2]

с [0]	=	0.34785485
с[1]	=	с [0]
С [2]	=	0.65214515
сСЗ]	=	с [2]
и ==	5:
х[0]	=	0.90617985
х[1]	=	- х [0]
х[2]	=	0.53846931
х[3]	=	- х[2]
х [4]	=	0.
с [0]	=	0.23692689
с[1]	=.	с [0]
с [2]	=	0.47862867
с [3]	=	с [2]
с [4]	=	0.56888889
п ==	6:
х[0]	=	0.93246951
х[1]	=	- х[0]
х [2]	=	0.66120939
х[3]	=	- х [2]
х[4]	=	0.23861919
х[5]	=	- х [4]
с [0]	=	0.17132449
с[1]	=	с [0]
С [2]	=	0.36076157
с [3]	=	с [2]
с [4]	=	0.46791393
с [5]	=	с [4]
п ==	7:
х [0]	=	0.94910791
х[1]	=	- х[0]
х[2]	=	0.74153119
х[3]	=	- х[2]
х[4]	=	0.40584515
х[5]	=	- х [4]
х[6]	=■	0.
с [0]	=	0.12948497
с[1]	=	с [0]
с [2]	=	0.27970539
с[3]	=	с [2]
с [4]	=	0.38183005
с [5]	=	с [4]
с [6]	=	0.41795918
п == 8:
х[0]	=	0.96028986

х[1] = - х[0] х[2] = 0.79666648 х[3] = - х[2] х[4] = 0.52553241 х[5] = - х [43 х[6] = 0.18343464 х[7] = - х[6] с[0] = 0.10122854 с[1] = с [0] с[2] = 0.22238103 с [3] = с [2] с[4] = 0.31370665 с [5] = с [4] с[6] = 0.36268378 с [7] = с [6] cl = (Ь + а)/2. с2 = (Ь - а)/2. sum = 0.
for i in range(n):
sum = sum + c[i]*f(cl + c2*x[i]) return c2*sum
Решение задачи приближенного вычисления рассматриваемого интеграла да­ется следующей программой.