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

7Г kh ~2~



Здесь

_к = л², \_к
есть собственные значения непрерывной и сеточной задачи соответственно. Использование функции ехрт2() вместо ехрш() дает те же результаты, _а вот ехртЗО - ничего разумного (из-за плохой сходимости ряда Тейлора при вычислении матричной экспоненты).
Для полной поддержки вычислительных алгоритмов линейной алгебры мож­но добавить итерационные методы решения систем линейных уравнений. Со­ответствующие функции можно найти в модуле решения задач с разрежен­ными матрицами sparse.linalg (ранее они находились в linalg). Пример исполь­зования метода сопряженных градиентов для ранее рассмотренной задачи (см. рис. 3.34) дается следующей программой.
import пшпру as пр from scipy.sparse import linalg import matplotlib.pyplot as pit m = 25
h = 1. / (m+1)
A = np.mat(np.zeros((m, m), ’float’)) for i in range(l, m):
A[i-l,i] = - 1. / h**2 A[i,i] = 2. / h**2 A [i,i-1] = - 1. / h**2 A[0,0] = 2. / h**2 A[m-1,m-l] = 2. / h**2 b = np.mat(np.ones((m), ’float’)).T help(linalg.eg)
y, info = linalg.cg(A, b, tol=l.e-06)
x = np.linspace(h, m*h, m)
plt.plot(x, y, x, x*(l-x)/2, ’:’)
pit.xlabel(’$x$’)
plt.ylabel(’$y$’)
pit.show()
Модуль sparse.linalg помимо метода сопряженных градиентов (функция cgO) включает другие итерационные методы решения систем линейных уравнений:

• 
  cgs() — квадратичный метод сопряженных градиентов;
  
• 
  qmr() — метод квазиминимальных невязок;
  
• 
  gmresO — обобщенный метод минимальных невязок;
  
• 
  bicgO — метод бисопряженного градиента;
  
• 
  bicgstabO — метод бисопряженного градиента со стабилизацией.
  

Эти итерационные методы применяются при решении задач как с полными? так и с разреженными матрицами.
Интерполяция
Модуль interpolate пакета SciPy предоставляет инструменты для решения стандартных задач интерполяции и сглаживания полиномами и сплайнами сеточных функций одной переменной. Имеются также функции для интерпо­ляции двумерных данных, заданных на регулярной сетке или в произвольных -точках плоскости.
Пусть в узлах х₀, х±,..., х_т функция у = у(х) принимает заданные значения 2/oj 2/1» • • ч Ут • При полиномиальной интерполяции строиться полином р(х) = a_Qx^m + aix^m~^l + ....+ a_m, для которого р(х_к) = у_к,к = 0,1,т. Для построе­ния интерполирующего полинома можно воспользоваться функцией lagrangeO. Более устойчивые алгоритмы вычисления значений интерполирующего поли­нома реализованы в функциях barycentric_interpolate() и krogh_interpolate (). Приведем пример (рис. 3.39) построения интерполирующего полинома для функции Рунге у = (I + 25х²)^-1.


Рис. 3.39 Аппроксимация полиномами



import пшпру as пр from scipy import interpolate import matplotlib.pyplot as pit ^{x 55} np.linspace(-l, 1., 200)
pit.plot(х, 1./(1 + 25*x*x), label=’exact’)
xl = np.linspace(-l, 1., 7)
yl = l./(l + 25*xl*xl)
pi = interpolate.lagrange(xl, yl)
pit.plot(x, pl(x), label=^,m=6^,)
x2 = np.linspace(-l, 1., 11)
y2 = 1./Ц + 25*x2*x2)
yb = interpolate.barycentric_interpolate(x2, y2, x)
pit.plot(x, yb, label=’m=10’)
pit.grid(True)
pit.Iegend(loc=0)
pit. showO
В практике научных вычислений гораздо большее внимание уделяется кусоч­но-полиномиальным аппроксимациям (сплайнам). Основной функцией одно­мерной интерполяции является interpldO, первый аргумент которой есть массив узлов, второй — массив значений. Параметр kind задает тип кусочно- полиномиальной интерполяции. При kind = ’linear’ (значение по умол­чанию) используется кусочно-линейная интерполяция (линейные сплайны), kind = ’cubic’— кубические сплайны, значения kind = ’nearest’, ’zero’ связаны с кусочно-постоянными аппроксимациями. Целое значение kind со­ответствует степени интерполирующего полинома между узлами. При зада­нии в узлах интерполяции производных кусочно-полиномиальная аппрокси­мация строится с помощью функции piecewise_polynomial_interpolate().
Приведем пример использования функции interpldO при интерполяции функции Рунге (см. рис. 3.40).
import numpy as np from scipy import interpolate import matplotlib.pyplot as pit x = np.linspace(-l, 1., 200)
pit.plot(x, l./(l + 25*x*x), label=’exact’)
xl = np.linspace(-l, 1., 11) yl = l./(l + 25*xl*xl) pit.scatter(xl, yl)
fl = interpolate.interpld(xl, yl, kind=’nearest’)
pit.plot(x, fl(x), label=’nearest’)
f2 = interpolate.interpld(xl, yl, kind=3)
plt.plot(x, f2(x), label=’cubic’)
pit.grid(True)
pit.Iegend(loc=0)
pit. showO
Модуль interpolate содержит много инструментов для работы со сплайнами. В частности, можно не только вычислить значение сплайна в заданной точке, но и найти коэффициенты сплайна, вычислить производную сплайна, пай-


Рис. 3.40 Кусочно-полиномиальная аппроксимация



ти интеграл на заданном интервале, корни сплайна. В качестве первого ил­люстративного примера рассмотрим построение сглаживающего Б-сплайна (рис.3.41). Б-сплайи строиться с помощью функции splrepO, для вычис­ления значений сплайна используется splev(). Параметр s есть параметр сглаживания.
import numpy as np from scipy import interpolate import matplotlib.pyplot as pit x = np.linspace(0., 1., 200)
pit. plot (x, np. cos(3*np.pi*x)*np.exp(-x), ’label=’exact’)
xl = np.linspace(0., 1., 11)
yl = np.cos(3*np.pi*xl)*np.exp(-xl)
pit.scatter(xl, yl)
tckl = interpolate.splrep(xl, yl, s=0) si = interpolate.splev(x, tckl, der=0) plt.plot(x, si, label=’s=0’)
tck2 = interpolate.splrep(xl, yl, s=0.05) s2 = interpolate.splev(x, tck2, der=0) pit.plot(x, s2, label=^,s=0.05^,)





tck3 = interpolate.splrep(xl, yl, s=l.) s3 = interpolate.splev(x, tck3, der=0) pit.plot(x, s3, label=’s=l.’)
pit.grid(True) pit.legend(loc=0) pit. showO
Для интерполяции двумерных данных на двумерной сетке можно воси зоваться функциями bisplrepO (нахождение Б-сплайна) и bisplevO числение значений двумерного Б-сплайна). На рис. 3.42 показана исхо, сеточная функция па сетке 11x11, результат интерполяции (сетка 101 х — на рис. 3.43.
import numpy as np
from scipy import interpolate
import matplotlib.pyplot as pit
x, у = np.mgrid[0:1:llj, 0:1:11j]
z = y*(l-y)*x
pit.figure(1)
pit.pcolor(x, y, z, стар⁼⁾gray’) pit.colorbar()
0.150
0
0.225 ^ 0.200 0.175

Рис. 3.42 Сеточная функция
.125