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

sg = sglist[k]
pit.plot(xl, yl, sg, label=sl) pit .xlabelCx’) pit.grid(True) pit.legend(loc=0) pit.show()
Расчетные данные, приведенные на рис. 9.1, демонстрируют недостатки по­линомиальной аппроксимации при больших п.
Упражнение 9.2 Напишите программу для интерполирования данных на основе естественного интерполяционного кубического сплайна. С помотаю этой программы интерполируйте данные в равноотстоящих узлах для функции Руте f(x) = (1 + 25а;²)^-1 на интервале [—1,1] при п — 4,6,10.
Сначала получим расчетные формулы для коэффициентов естественного ку­бического сплайна.
Введем обозначения
= х_{ - i, i = 1,2,..., n.
Для кубического сплайна 5₃(ж) с учетом представления (9.6) получим следу­ющую систему уравнений:
аг = /(х_{), г = 0,1,..., п — 1, (9.8)
а* + bihi+i + 7^i+i + "q"^+i ⁼ i ⁼ 0,1, • • •,ⁿ — 1, (9.9)
bi +Cih_i+1 + ^/t-₊₁ = b_i+u ?; = 0,l,...,n-2, (9.10)
Ci + dihi+i = Cj+i, i = 0,1,... ,n - 2, (9.11)
Co = 0, c_n__x + dn-xK = 0. (9.12)
Формально доопределим Cn — 0, тогда из (9.11) и второго условия (9.12) получим
di = , г = 0,1,... ,п- 1, (9.13)
h_i+1
а вместо (9.12) будем иметь
с₀ = 0, с_п = 0. (9.14)
Подстановка (9.8), (9.13) в (9.9) дает следующее представление для коэффи' циентов bi:
bi = _ ^l±i(_Ci+1 + 2_Ci), (9.15)
Ы+1 6
г = 1,2,..., n - 1.
С учетом (9.13), (9.15) соотношения (9.10) приводят к уравнению
Ci-ihi + 2Ci(hi + hi+1) -f- Cf-j-i/ij+i —
=
(9.1G)

q( _ f(^xi) -
\ h_i+i Ы ) ’
i — 1,2,..., n - 1.
Тем самым приходим к линейной системе уравнений (9.14), (9.16) с трехдиа­гональной матрицей с диагональным преобладанием. Решение этой системы всегда существует и единственно. Другие коэффициенты сплайна определя­ются в соответствии с (9.8), (9.13), (9.15).
В модуле spline функция coef spline () предназначена для вычисления коэф­фициентов интерполяционного кубического сплайна, а функция spline() — для вычисления значения кубического сплайна в заданной точке. Для нахо­ждения коэффициентов сплайна используется функция solveLU3() (решение системы уравнений с трехдиагональной матрицей) из модуля 1иЗ.
Модуль spline
import numpy as np from lu3 import solveLU3 def spline(x, y, c, xO):
и и и
Evaluates cubic spline at xO.
и и и
def find(x, xO):
ч и и и
Find the segment spanning xO
и и и
iL = 0
iR = len(x) - 1 while iR - iL > 1:
, i = (iL + iR) / 2
if xO < x [i]: iR = i else:
iL = i return iL i = find(x, xO) h = x[i+l] - x[i]
yO = ((x[i+l] - x0)**3/h - (x[i+1] - x0)*h) * c[i] /6. \ + ((xO - x[i])**3/h - (xO - x[i])*h) * c[i+l] /6. \
+ (y[i]*(x[i+l] - xO) + y[i+l]*(x0 - x[i])) / h return yO
^def coefspline(x, y):
к и и
Computes the coefficients of cubic spline.
и и и
n = len(x) - 1 a = np.ones((n+1), ’float’) b = np.zeros((n+1), ’float’) c = np.zeros((n+1), ’float’) f = np.zeros((n+1), ’float’) for i in ranged, n):
a[i] = 2. *(x[i+l] - x[i-l]) b[i] = x[i+l] - x[i] c[i] = x[i] - x[i-l]
f[i] = 6.*(y[i+l] - у[i]) / (x[i+l] - x[i]) \
- 6.*(y[i] - y[i-l]) / (x[i] - x[i-l]) return solveLU3(a, b, c, f)
Решение задачи интерполяции сплайнами для функции Рунгс дается следу­ющей программой.


Рис. 9.2 Сплайи-интерполяция функции Рунге f(x) = (1 + 25ж²) ¹ ехег -9 2 ру



imp^ort numpy as np import matplotlib.pyplot as pit from spline import spline, coefspline def f(x):
return 1. / (1 + 25*x**2)
_a = -1. b = 1- m = 200
xl = np.linspace(a, b, m) yl = f(xl) pit.plot(xl, yl) nList = [4, 6, 10] sglist = , ’:’]
for к in range(len(nList)): n = nList[k]
x = np.linspace(a, b, n+1) у = f(x)
pit.scatter(x, у, тагкег=^,о^>) c = coefspline(x, y) for i in range(m):
yl [i] = spline(x, y, c, xl [i]) si = ^,n=^> + str(n) sg = sglist [k]
pit.plot(xl, yl, sg, label=sl) pit.xlabel(’x’) pit.grid(True) pit.legend(loc=0) pit. show ()
Преимущества интерполяции сплайнами иллюстрируются рис. 9.2.
* 9.4 Задачи
Задача 9.1 Напишите программу для пахоэюдеиия значения интерполиру­ющего полинома в точке х па основе рекуррентных соотношений (алгоритм Невилля):
Ри(х) =Уг, 0 < г < п,
_Pij(x) = 1* ⁺ ^ ~ ^x)P'-iA^x\ о < i < j < n.
Xi Xj
С помощью этой программы найдите значения в точках х = 1.5, 2.5 при интерполировании функции /(х) = (arctan(l+:г²))^-1 при использовании рав­ноотстоящих узлов на интервале [—3,3] при п = 4, 6,10.
Задача 9.2 Напишите программу для интерполирования данных на оспо. ее кусочно-линейного восполнения (интерполяционного линейного сплайна) С помощью этой программы интерполируйте данные в равноотстоящие узлах для функции f(x) = (arctan(l + х²))~¹ на интервале [—3,3] при п - 4,6,10.
Задача 9.3 Напишите программу для приблиэюепия сеточной функции y_i i = 0,1,...,п полиномом р_т(х) = с₀ + с_хх + ... + с_тх^т, т < п методов наименьших квадратов. Работоспособность программы проиллюстрируйте па примере аппроксимации
Уг = 1 -cosfe),
Xi = ih, h = 1/n, г = 0,1, ..:,n для n — 10 и m = 1,2,3.
Задача 9.4 Напишите программу для приблиоюения сеточной функции у_и г = 0,1,...,п функцией £(х) = ае^Ьх методом наименьших квадратов. С по­мощью этой программы аппроксимируйте сет,очную функцию
Vi = 1 - COS(Xi),
Xi=ih, h=l/n, г = 0,
для n = 10.

Численное интегрирование
З
Основные обозначения


f(x) — подынтегральная функция д(х) — весовая функция


х₀ < х_х < ... < х_п — узлы квадратурной формулы Cj, г = 0,1,... ,п — коэффициенты квадратурной формулы L_n(x) — интерполяционный многочлен n-го порядка

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

1. 
   Задачи приближенного вычисления интегралов
   

Рассматривается задача приближенного вычисления определенных интегра­лов
( 10.1)
па некотором классе функций f(x) с заданной весовой функцией р(.т).
С этой целью подынтегральная функция задается в отдельных точках х, отрезка |а,6], i = 0,1,...,п. Под квадратурной формулой понимается при­ближенное равенство
( 10.2)

в
котором Ci, г = 0,... , п — коэффициенты квадратурной формулы. Через
определим погрешность квадратурной формулы.
Минимизация погрешности (увеличение точности) квадратурной формулы на выбранном классе функций может достигается прежде всего за счет вы­бора коэффициентов квадратурной формулы а также за счет выбора узлов интегрирования.

2. 
   Алгоритмы приближенного вычисления интегралов
   

Рассматриваются простейшие квадратурные формулы прямоугольников, трапеций и Симпсона для приближенного вычисления определенных инте­гралов. Строятся квадратурные формулы интерполяционного типа с задан­ными узлами квадратурной формулы, обсуждаются вопросы оптимизации квадратурные формулы за счет выбора узлов.
Классические квадратурные формулы составного типа
Б
удем рассматривать задачу вычисления интеграла
т
.е. весовая функция g(x) = 1. Представим интеграл в виде суммы по частич­ным отрезкам:
Приведем некоторые простейшие квадратурные формулы для приближенно­го вычисления интеграла на частичном отрезке Х{].
Формула
( 10.3)
называется формулой прямоугольников на частичном отрезке где
Xi-1/2 = (Zi-1 + Xi)/2.
Для случая равномерной сетки
uj = {ж | х = Xi = а + ih, г = 0,1,..., гг, nh = b — а}
суммирование (10.1) по г от 1 до N дает составную формулу прямоугольников
[ f(x)dx « Y, f{xi-i_/2)h.
Ja ,_i
П