Интерполяционные сплайны на прямоугольных сетках

Введение

математика интерполяционный программа многочлен

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

Теория приближений активно используется при построении численных алгоритмов, а также при сжатии информации.

В вычислительной математике существенную роль играет интерполяция функций. Она представляет собой один из простых и старейших способов восстановления неизвестной функции на некотором множестве по набору известных значений в отдельных точках этого множества. Для нахождения многих функций, оказывается, эффективно приблизить их алгебраическими многочленами. Наиболее известная интерполяционная формула была открыта Лагранжем в 1975 г. В Scopus опубликовано немало работ связанных с полиномом Лагранжа (см. [1], [2], [3]).

1. Построение интерполяционного многочлена из пространства

Дана функция f(x,y), необходимо найти интерполяционный многочлен и найти точность приближения.

Общая формула интерполяционного многочлена Лагранжа p(x,y) ϵ, имеет вид

 ϵ R.

На квадрате Q[0;1]² строятся узлы. D = dim =(m+1)(n+1) - минимально число узлов интерполяции, гарантирующее единственное решение задачи.

Выбор узлов выглядит так,

Рис. 1 Выбор узлов

Ось X - на отрезке [0;1] делится на m равных отрезков,

Ось Y - на отрезке [1;0] делится на n равных отрезков, а точки пересечения и будут узлами интерполяции.

:=(; i=0,…,m; j=0,…,n;

 - узлы сетки, получающиеся при разбиении координатных сторон квадрата.

Пусть f ϵ C(Q),  Задача интерполяции заключается в нахождении такого многочлена p ϵ , для которого

p( = ; i=0,…,m; j=0,…,n; (1)

Решение задачи (1) дается следующим аналогом формулы Лагранжа. Положим

Тогда многочлен

принадлежит . Тогда для нахождения точности разделим отрезок [0;1] и [1;0] на 1000 (к примеру, можно делить на большее число, чем больше число, тем больше вероятность достигнуть наилучшей точности) равных по длине отрезков. Получим по x и y 1000 новых точек равноудаленных друг от друга. Тогда нахождение точности сводится к формуле

. Приближение кусочно-полиномиальными функциями

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

Введем теперь множество π) кусочно-полиномиальных функций степени (m,n), подчинённых разбиению π. Каждая функция из этого множества на прямоугольнике из разбиения π совпадает с некоторым многочленом p из пространства .

Рассмотрим вопрос о нахождении кусочно-полиномиальной функции g* наилучшего приближения. Такая функция определяется соотношением

Теоретически вопрос решается следующим образом. Нужно на каждом прямоугольнике из разбиения π построить для функции f многочлен наилучшего приближения из пространства . Но так как эта задача является весьма сложной, то в моей работе я использую интерполяционный многочлен Лагранжа для приближения данной функции (построение интерполяционного многочлена рассматривалось ранее в данной работе). Таким образом, если разбиение π заранее задано, то нетрудно построить кусочно-полиномиальную функцию, которая дает приближение почти как наилучшая.

Перейдя к более сложной ситуации, когда разбиение заранее не задано: при этом нужно найти кусочно-полиномиальную функцию, приближающую заданную функцию,  c точностью ε. Разбиение, которому будет подчинена кусочно-полиномиальная функция, зависит и от функции , и от точности ε. Такой процесс приближения называется адаптивным, так как разбиение подстраивается под характер особенностей функции. В частности, если у функции нет особенностей, то разбиение будет состоять из почти равных прямоугольников. Ниже будет рассмотрен один из таких алгоритмов.

Уточним постановку задачи. На квадрате  задана непрерывная функция : кроме того, заданы числа  и ε>0. Нужно построить такое разбиение единичного квадрата π= π(), что

3. Описание алгоритма программы и примеры её реализации

На квадрате  построим интерполяционный многочлен  из пространства  и найдем погрешность приближения. Если погрешность меньше , то  искомая функция. В этом случае разбиение π состоит из единственного квадрата . Пусть погрешность приближения (т.е. ) больше . Разделим квадрат  на четыре равных квадрата:

[

[

На каждом из этих квадратов построим интерполяционный многочлен (i,j=0,1). Вычислим погрешности:

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

Примеры работы программы:

Здесь приведен пример, когда многочлен , то есть состоит из одних констант. Многочлен здесь красного цвета, а сама функция зеленого. При уменьшении точности количество квадратов будет расти, а соответственно и число многочленов:

Далее рассмотрим пример, когда меняется не точность, а степени многочлена:

При увеличении степеней многочлена видно, что многочлен всё сильнее приближается к исходной функции.

Заключение

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

В данной работе мною было изучено построение интерполяционного многочлена Лагранжа двух переменных x и y, степени не больше m и n соответственно. Реализован алгоритм приближения интерполяционными сплайнами на прямоугольных сетках. Исходя из тестов программы, можно сказать, что при увеличении степеней многочлена и при фиксированной точности, мы достигаем этой точности за меньшее число шагов. При фиксированных степенях полинома и при уменьшении точности, мы достигаем этой точности за большее число шагов, что приводит к более долгому процессу нахождения всех полиномов, т.к. их становится больше.

Список литературы

1.   Невский М.В., Иродова И.П. Некоторые вопросы теории приближения функций. Ярославль, 1999.

2.   Брудный Ю.А. Теория приближения. Ярославль, 1981.

3.   « И.А. Шакиров, “О влиянии выбора узлов лагранжевой интерполяции на точные и приближенные значения констант Лебега”, Сиб. матем. журн., 55:6 (2014), 1404-1423».

4.      «И.А. Шакиров, “О тригонометрическом интерполяционном полиноме Лагранжа, имеющем минимальную норму как оператор из C2π в C2π”, Изв. вузов. Матем., 2010, №10, 60-68».

.        «Е.А. Волков, “О применении интерполяционного многочлена Лагранжа при решении методом сеток задачи Дирихле для уравнения Пуассона”, Ж. вычисл. матем. и матем. физ., 4:3 (1964), 466-472».

Приложение

Описание программной среды и языков программирования:

Средой программирования была выбрана Microsoft Visual Studio 2013. такой выбор был сделан в силу ее ориентации под программные продукты, работающие в ОС Windows, которая в свою очередь является самой популярной операционной системой. Языком, на котором написана программа стал C#, с подключенной компонентой Windows Forms. Это взаимодействие позволило представить программу пользователю в дружественном и интуитивно понятном интерфейсе.

Код программы:

CalculateCore.cs

using System;System.Collections.Generic;System.Linq;

Interpolation

{Point3D

{double X { get; set; }double Y { get; set; }double Z { get; set; }

CalculationCore

{int _m;int _n;int M

{{ return _m - 1; }{ _m = value + 1; }

}

int N

{{ return _n - 1; }{ _n = value + 1; }

}

double CalculateDeltaFunctionsNorm(Func<double, double, double> func1, Func<double, double, double> func2,x0, double y0, double x1, double y1)

{values1 = CalculatePointsForFunctionOnSquare(func1, x0, y0, x1, y1, 10, 10).Select(v => v.Z).ToArray();values2 = CalculatePointsForFunctionOnSquare(func2, x0, y0, x1, y1, 10, 10).Select(v => v.Z).ToArray();

values1.Select((t, i) => Math.Abs(t - values2[i])).Concat(new[] {0.0}).Max();

}

List<Point3D> CalculatePointsForFunctionOnSquare(Func<double, double, double> function, double x0, double y0,x1, double y1, int xPointCount, int yPointCount)

{result = new List<Point3D>();xStep = (x1 - x0)/xPointCount;yStep = (y1 - y0)/yPointCount;

(int i = 0; i <= xPointCount; i++)(int j = 0; j <= yPointCount; j++).Add(new Point3D()

{= x0 + i*xStep,= y0 + j*yStep,= function(x0 + i*xStep, y0 + j*yStep)

});

result;

}

Func<double, double, double> CalculatePolynomOnSquareByFunction(Func<double, double, double> function, double x0, double y0,x1, double y1)

{xStep = (x1 - x0) / _m;yStep = (y1 - y0) / _n;

lComponents = new List<Func<double, double>>();mComponents = new List<Func<double, double>>();(int i = 0; i <= _m; i++).Add(CalculateComponent(i, x0, x1, _m));

(int j = 0; j <= _n; j++).Add(CalculateComponent(j, y0, y1, _n));<double, double, double> polynom = (x, y) =>

{r = 0;

(M == 0 || N == 0)= function((x0 + x1)/2, (y0 + y1)/2);

{(int i = 0; i <= _m; i++)(int j = 0; j <= _n; j++)+= function(x0 + i*xStep, y0 + j*yStep)*lComponents[i](x)*[j](y);

}

r;

};

polynom;

}

Func<double, double> CalculateComponent(int k, double p0, double p1, int count)

{x =>

{result = 1.0;

step = (p1 - p0)/count;

for (int i = 0; i <= count; i++)(i!= k)*= (x - (p0 + i*step))/((k - i)*step);

result;

};

}

}

}

MainForm.cs

using System;System.Collections.Generic;System.Linq;System.Threading.Tasks;System.Windows.Forms;Tao.FreeGlut;Tao.OpenGl;

Interpolation

{partial class MainForm: Form

{readonly Func<double, double, double> _baseFunction = (x, y) => Math.Sin(x*10) + Math.Sin(y);double _eps = 0.001;int QuadCount;<List<Point3D>> _functionPoints;<List<Point3D>> _polynomPoints;MainForm()

{();

.InitializeContexts();

();

}

void InitGlut()

{.glutInit();.glutInitDisplayMode(Glut.GLUT_RGB | Glut.GLUT_DOUBLE | Glut.GLUT_DEPTH);

.glClearColor(255, 255, 255, 1);

.glViewport(0, 0, OglOutput.Width, OglOutput.Height);

.glMatrixMode(Gl.GL_PROJECTION);.glLoadIdentity();.gluPerspective(45, (float)OglOutput.Width / OglOutput.Height, 0.1, 200);.glMatrixMode(Gl.GL_MODELVIEW);.glLoadIdentity();

.glEnable(Gl.GL_DEPTH_TEST); .glPolygonMode(Gl.GL_FRONT_AND_BACK, Gl.GL_FILL);

}

private async void StartQuadsInterpolation()

{

_eps = double.Parse(textBoxEps.Text);calculations = new CalculationCore {M = int.Parse(textBoxM.Text), N = int.Parse(textBoxN.Text)};quadCount = 0;maxDelta;Task.Factory.StartNew(() =>

{

{++;quadSize = 1.0/Math.Pow(2, quadCount);deltas = new List<double>();

(int i = 0; i < Math.Pow(2, quadCount); i++)(int j = 0; j < Math.Pow(2, quadCount); j++)

{polynom = calculations.CalculatePolynomOnSquareByFunction(_baseFunction,*quadSize, j*quadSize,

(i + 1)*quadSize, (j + 1)*quadSize);.Add(calculations.CalculateDeltaFunctionsNorm(_baseFunction, polynom,*quadSize, j*quadSize,

(i + 1)*quadSize, (j + 1)*quadSize));

}

= deltas.Max();

} while (maxDelta > _eps);

});

.Text = "Number decompositions: " + (quadCount);= (int)Math.Pow(2, quadCount);(calculations);();

void DrawGraph()

{.glClear(Gl.GL_COLOR_BUFFER_BIT | Gl.GL_DEPTH_BUFFER_BIT);

.glLoadIdentity();

.glPushMatrix();.glTranslated(0, -0.27, -10.0/(zoomBar.Value+1));.glRotated(90, -1, 0, 0);

.glRotated(xAngleBar.Value*5, 1, 0, 0);.glRotated(yAngleBar.Value*5, 0, 1, 0);.glRotated(zAngleBar.Value*5, 0,0, 1);

();();

.glPopMatrix();.glFlush();.Invalidate();

}void DrawPolygons()

{(checkBoxFunction.Checked)

{.glColor3f(0.0f, 0.7f, 0.0f);(_functionPoints);

}

(checkBoxPolynom.Checked)

{.glColor3f(0.7f, 0.0f, 0.0f);(_polynomPoints);

}

}

void CalculatePoints(CalculationCore calculations)

{

_functionPoints = new List<List<Point3D>>();

_polynomPoints = new List<List<Point3D>>();

quadSize = 1.0/QuadCount;(int i = 0; i < QuadCount; i++)(int j = 0; j < QuadCount; j++)

{

_functionPoints.Add(calculations.CalculatePointsForFunctionOnSquare(_baseFunction,*quadSize, j*quadSize,

(i + 1)*quadSize, (j + 1)*quadSize, 10, 10));

var polynom = calculations.CalculatePolynomOnSquareByFunction(_baseFunction,*quadSize, j*quadSize,

(i + 1)*quadSize, (j + 1)*quadSize);

_polynomPoints.Add(calculations.CalculatePointsForFunctionOnSquare(polynom,*quadSize, j*quadSize,

(i + 1)*quadSize, (j + 1)*quadSize, 10, 10));

}

}

static void DrawSurface(List<List<Point3D>> points)

{(int k = 0; k < points.Count; k++)

{.glBegin(Gl.GL_TRIANGLES);

{internalCount = (int) Math.Sqrt(points[k].Count);(int i = 0; i < internalCount - 1; i++)(int j = 0; j < internalCount - 1; j++)

{.glVertex3d(points[k][i*internalCount + j].X, points[k][i*internalCount + j].Y,[k][i*internalCount + j].Z);.glVertex3d(points[k][i*internalCount + j + 1].X, points[k][i*internalCount + j + 1].Y,[k][i*internalCount + j + 1].Z);.glVertex3d(points[k][(i + 1) * internalCount + j].X, points[k][(i + 1) * internalCount + j].Y,[k][(i + 1) * internalCount + j].Z);.glVertex3d(points[k][i * internalCount + j + 1].X, points[k][i * internalCount + j + 1].Y,[k][i * internalCount + j + 1].Z);.glVertex3d(points[k][(i + 1) * internalCount + j].X, points[k][(i + 1) * internalCount + j].Y,[k][(i + 1) * internalCount + j].Z);.glVertex3d(points[k][(i + 1) * internalCount + j + 1].X, points[k][(i + 1) * internalCount + j + 1].Y,[k][(i + 1) * internalCount + j + 1].Z);

}

}.glEnd();

}

}

void DrawGrid()

{.glBegin(Gl.GL_LINES);

{.glColor3f(1.0f, 0, 0);.glVertex3d(0, 0, 0);.glVertex3d(2, 0, 0);.glColor3f(0, 1.0f, 0);.glVertex3d(0, 0, 0);.glVertex3d(0, 2, 0);

.glColor3f(0, 0, 1.0f);.glVertex3d(0, 0, 0);.glVertex3d(0, 0, 2);

Gl.glColor3f(0.5f, 0.5f, 0.5f);(int i = 0; i < QuadCount; i++)

{.glVertex3d((1.0/QuadCount)*(i + 1), 0, 0);.glVertex3d((1.0/QuadCount)*(i + 1), 1, 0);.glVertex3d(0, (1.0/QuadCount)*(i + 1), 0);.glVertex3d(1, (1.0/QuadCount)*(i + 1), 0);

}

}.glEnd();

}

void drawButton_Click(object sender, EventArgs e)

{();

}

void AngleBar_Scroll(object sender, EventArgs e)

{(_functionPoints!= null && _polynomPoints!= null)();

}

}

}

Program.cs

using System;System.Collections.Generic;System.Linq;System.Threading.Tasks;System.Windows.Forms;

Interpolation

{class Program

{

/// <summary>

/// The main entry point for the application.

/// </summary>

[STAThread]void Main()

{.EnableVisualStyles();.SetCompatibleTextRenderingDefault(false);

Application.Run(new MainForm());

}

}

}