Реализация на языке C# Платформа .NET

⇐ ПредыдущаяСтр 2 из 2

// Интерполирование функций естественными кубическими сплайнами using System; class CubicSpline{ SplineTuple[] splines; // Сплайн // Структура, описывающая сплайн на каждом сегменте сетки private struct SplineTuple { public double a, b, c, d, x; } // Построение сплайна // x - узлы сетки, должны быть упорядочены по возрастанию, кратные узлы запрещены // y - значения функции в узлах сетки // n - количество узлов сетки public void BuildSpline(double[] x, double[] y, int n) { // Инициализация массива сплайнов splines = new SplineTuple[n]; for (int i = 0; i < n; ++i) { splines[i].x = x[i]; splines[i].a = y[i]; } splines[0].c = splines[n - 1].c = 0.0; // Решение СЛАУ относительно коэффициентов сплайнов c[i] методом прогонки для трехдиагональных матриц // Вычисление прогоночных коэффициентов - прямой ход метода прогонки double[] alpha = new double[n - 1]; double[] beta = new double[n - 1]; alpha[0] = beta[0] = 0.0; for (int i = 1; i < n - 1; ++i) { double hi = x[i] - x[i - 1]; double hi1 = x[i + 1] - x[i]; double A = hi; double C = 2.0 * (hi + hi1); double B = hi1; double F = 6.0 * ((y[i + 1] - y[i]) / hi1 - (y[i] - y[i - 1]) / hi); double z = (A * alpha[i - 1] + C); alpha[i] = -B / z; beta[i] = (F - A * beta[i - 1]) / z; } // Нахождение решения - обратный ход метода прогонки for (int i = n - 2; i > 0; --i) { splines[i].c = alpha[i] * splines[i + 1].c + beta[i]; } // По известным коэффициентам c[i] находим значения b[i] и d[i] for (int i = n - 1; i > 0; --i) { double hi = x[i] - x[i - 1]; splines[i].d = (splines[i].c - splines[i - 1].c) / hi; splines[i].b = hi * (2.0 * splines[i].c + splines[i - 1].c) / 6.0 + (y[i] - y[i - 1]) / hi; } } // Вычисление значения интерполированной функции в произвольной точке public double Interpolate(double x) { if (splines == null) { return double.NaN; // Если сплайны ещё не построены - возвращаем NaN } int n = splines.Length; SplineTuple s; if (x <= splines[0].x) // Если x меньше точки сетки x[0] - пользуемся первым эл-тов массива { s = splines[0]; } else if (x >= splines[n - 1].x) // Если x больше точки сетки x[n - 1] - пользуемся последним эл-том массива { s = splines[n - 1]; } else // Иначе x лежит между граничными точками сетки - производим бинарный поиск нужного эл-та массива { int i = 0; int j = n - 1; while (i + 1 < j) { int k = i + (j - i) / 2; if (x <= splines[k].x) { j = k; } else { i = k; } } s = splines[j]; } double dx = x - s.x; // Вычисляем значение сплайна в заданной точке по схеме Горнера (в принципе, "умный" компилятор применил бы схему Горнера сам, но ведь не все так умны, как кажутся) return s.a + (s.b + (s.c / 2.0 + s.d * dx / 6.0) * dx) * dx; }}

Реализация на Python 2.7

from collections import defaultdictfrom matplotlib import mlabimport bisect, pylab, math def drange(start, stop, step): while start < stop: yield start; start += step; class Dot: def __init__(self, x, y): self.x, self.y = [x, y] class Tuple: a, b, c, d, x = [0., 0., 0., 0., 0.] def buildSpline(dots): for i in range(len(dots)): splines[i].x, splines[i].a = dots[i].x, dots[i].y alpha, beta = [defaultdict(lambda: 0.), defaultdict(lambda: 0.)] for i in range(1, len(dots)-1): C = 4. * in_step F = 6. * ((dots[i + 1].y - dots[i].y) / in_step - (dots[i].y - dots[i - 1].y) / in_step) z = (in_step* alpha[i - 1] + C) alpha[i] = -in_step / z beta[i] = (F - in_step* beta[i - 1]) / z for i in reversed(range(1, len(dots) - 1)): splines[i].c = alpha[i] * splines[i+1].c + beta[i] for i in reversed(range(1, len(dots))): hi = dots[i].x - dots[i-1].x splines[i].d = (splines[i].c - splines[i-1].c) / hi splines[i].b = hi * (2.0 * splines[i].c + splines[i - 1].c) / 6.0 + (dots[i].y - dots[i-1].y) / hi def calc(x): distribution = sorted([t[1].x for t in splines.items()]) indx = bisect.bisect_left(distribution, x) if indx == len(distribution): return 0 dx = x - splines[indx].x return splines[indx].a + splines[indx].b * dx + splines[indx].c * dx**2 / 2. + splines[indx].d * dx**3 / 6.#============================================ in_func = lambda x: math.sin(x)in_min_x = 0in_max_x = 25 in_step = 2.5 out_min_x = 0out_max_x = 25out_step = 0.1 #============================================ #build modelsplines = defaultdict(lambda: Tuple())buildSpline(map(lambda x: Dot(x, in_func(x)), [x for x in drange(in_min_x, in_max_x+1, in_step)])) #print resultfor x in drange(out_min_x, out_max_x, out_step): print str(x) + ';' + str(calc(x)) #build graphicsxlist = mlab.frange (out_min_x, out_max_x, out_step)pylab.plot(xlist, [calc(x) for x in xlist])pylab.plot(xlist, [in_func(x) for x in xlist])pylab.show()

Квадратичный сплайн

Определение сплайна:

Сплайном Sm(x) называется функция определенная на [a,b] , p раз непрерывно дифференцируемая, и такая, что на каждом промежутке [xk-1,xk](k=1, 2, 3, ... ,n) это многочлен степени m
Разность d=m−p между степенью сплайна и показателем гладкости называется дефектом сплайна.

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

На отрезке [a, b] заданы (n+1) точка - узлы интерполяции.

x0,x1,x2, . . . ,xn

и известны значения функции y=f(x) в этих точках. Задана табличная функция:

f(x0)=y0, f(x1)=y1, f(x2)=y2, . . ., f(xn)=yn

Необходимо найти приближенное значение функции в точках, отличных от заданных.

3.1 Квадратичный сплайн дефекта 1:

Квадратичным сплайном дефекта 1 интерполирующем на отрезке [a,b] данную функцию f(x) называется функция S2(x) , такая что:

1) S2(xk)=fk - условие интерполяции.

2) S'2(k)(xk)=S'2(k+1)(xk) - равенство первых производных в узлах (гладкая стыковка)(k=1, 2, 3, ... ,n-2)

3) S'2(x0)=CилиS'2(xn)=C - краевое условие, С - произвольное число

Сформируем систему для 1-го условия:

a1x20+b1x0+c1=f0a1x21+b1x1+c1=f1a2x21+b2x1+c2=f1a2x22+b2x2+c2=f2. . . . .. . . . .anx2n-1+bnxn-1+cn=fn-1anx2n+bnxn+cn=fn

Для 2-го условия:

2a1x1+b1=2a2x1+b22a2x2+b2=2a3x2+b3. . . . . .2an-1xn-1+bn-1=2anxn-1+bn

Для 3-го условия: 2a1x0=C или 2anxn=C.

Объединив полученные системы в одну и решив её получим матрицу коефицентов ai,bi,ci,i=(1,2,. . .,n)

Формирование систем

Сформируем системы.

Для 1-го условия:

a1x20+b1x0+c1=f0a1x21+b1x1+c1=f1a2x21+b2x1+c2=f1a2x22+b2x2+c2=f2. . . . .. . . . .anx2n-1+bnxn-1+cn=fn-1anx2n+bnxn+cn=fn

Для 2-го условия:

2a1x1+b1=2a2x1+b22a2x2+b2=2a3x2+b3. . . . . .2an-1xn-1+bn-1=2anxn-1+bn

Для 3-го условия:

2a1x0=C или 2anxn=C

Объединив полученные системы в одну и решив её получим матрицу коефицентов

ai,bi,ci,i=(1,2,. . .,n)

Заключение

В вычислительной математике существенную роль играет интерполяция функций, т.е. построение по заданной функции другой (как правило, более простой), значения которой совпадают со значениями заданной функции в некотором числе точек. Причем интерполяция имеет как практическое, так и теоретическое значение. На практике часто возникает задача о восстановлении непрерывной функции по ее табличным значениям, например, полученным в ходе некоторого эксперимента. Для вычисления многих функций, оказывается, эффективно приблизить их полиномами или дробно-рациональными функциями. Теория интерполирования используется при построении и исследовании квадратурных формул для численного интегрирования, для получения методов решения дифференциальных и интегральных уравнений. Основным недостатком полиномиальной интерполяции является то, что она неустойчива на одной из самых удобных и часто используемых сеток - сетке с равноудаленными узлами. Если позволяет задача, эту проблему можно решить за счет выбора сетки с Чебышевскими узлами. Если же мы не можем свободно выбирать узлы интерполяции или нам просто нужен алгоритм, не слишком требовательный к выбору узлов, то рациональная интерполяция может оказаться подходящей альтернативой полиномиальной интерполяции.

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

Список использованной литературы

1.Роджерс Д.,Адамс Дж. Математические основы машинной графики. — М.: Мир, 2001.

2.Завьялов Ю. С., Леус В. А., Скороспелов В. А. Сплайны в инженерной геометрии. — М.: Машиностроение, 1985.

3.Лившиц Е.Д. Непрерывные E-выборки для приближения полиномиальными и рациональными сплайнами : Дис. … канд. физ.-мат. наук : 01.01.01 Москва.

4.Алберг Дж., Нильсон Э., Уолш Дж. — Теория сплайнов и ее приложения.

5.Винниченко Л. Ф. Экспоненциальные гистосплайны: предпосылки введения // Publishing house Education and Science s.r.o., конференция «Европейская наука XXI века», 2009.

6.Корнейчук, Н. П., Бабенко, В. Ф., Лигун, А. А. Экстремальные свойства полиномов и сплайнов / отв. ред. А. И. Степанец; ред. С. Д. Кошис, О. Д. Мельник, АН Украины, Ин-т математики. — К.: Наукова думка 1992.

7.Вершинин В. В., Завьялов Ю. С, Павлов Н. Н. Экстремальные свойства сплайнов и задача сглаживания. — Новосибирск: Наука, 1988, УДК 519.651

8.Роженко Александр Иосифович. Теория и алгоритмы вариационной сплайн-аппроксимации : Дис. … д-ра физ.-мат. наук : 01.01.07 : Новосибирск, 2003 231 c. РГБ ОД, 71:05-1/136

9.Шикин Е. В., Плис Л. И. Кривые и поверхности на экране компьютера. Руководство по сплайнам для пользователей. — М.: ДИАЛОГ-МИФИ, 1996.

10.Хакимов Б.В. Моделирование корреляционных зависимостей сплайнами на примерах в геологии и экологии. — С-Пб.: НЕВА, 2003.

Дата добавления: 2018-02-15; просмотров: 2302; Мы поможем в написании вашей работы!

Поделиться с друзьями:

⇐ Предыдущая 12

Мы поможем в написании ваших работ!