Приближение сплайнами

Цель лекции: Описать стандартные методы построения кубических сплайнов. Реализовать на языке C# процедуру интерполирования функций с помощью кубических сплайнов.

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

Наиболее успешными средствами интерполяции числовых функций являются методы, основанные на сплайн-интерполяции. Термин "сплайн" происходит от названия гибкой рейки, с помощью которой проводят гладкие кривые через заданные точки. В заданных точках сплайн закреплялся физически. В математике под сплайном понимают интерполирующую функцию, производные которой могут иметь разрывы в узловых точках. Классическим вариантом сплайна является кубические сплайны. Эти сплайны иногда называют сплайнами Шёнберга по имени математика, который ввел эти сплайны на современном уровне. Кубические сплайны представляют собой дважды непрерывно дифференцируемые функции, которые на отрезка между соседними узлами являются кубическими многочленами.

Пусть также дано разбиение отрезка $$[a,b]$$, и значения функции в узлах$$0=x_0<x_1<\dots<x_N=b,$$ $$f(x_n)=f_n,\quad n=0,\dots,N.$$ Кубическим сплайном называется такая функция $$S_N(x)$$, $$x\in[a,b]$$, что для этой функции выполнены следующие условия:

$$S_N(x_n)=f_n$$, $$n=0,1,\dots,N$$.

$$S_N(x)$$, $$S'_N(x)$$, $$S''_N(x)$$ являются непрерывными функциями.

На каждом отрезке $$[x_{n-1},x_n]$$, $$n=1,\dots,N$$ функция $$S_N$$ представима в виде$$S_N(x)=a^n_3x^3+a^n_2x^2+a^n_1x+a^n_0,\quad x\in[x_{n-1},x_n], n=1,\dots,N.$$

Как мы увидим в дальнейшем - для однозначного задания сплайна нам необходимо дополнительно добавить краевые условия на искомый сплайн. Мы будем рассматривать наиболее характерные условия:

$$S'_N(x_0)=m_0$$, $$S'_N(x_N)=m_N$$.

$$S''_N(x_0)=M_0$$, $$S''_N(x_N)=M_N$$.

Перейдем к вопросу построения сплайна. Разумеется для того, чтобы задать сплайн необходимо вычислить набор коэффициентов $$a^N_i$$, $$i=0,1,2,3$$ ; $$n=1,2,\dots,N$$. Однако этот путь является неэффективным. Дело в том, что эти коэффициенты не являются независимыми - на них наложены условия, чтобы обеспечить непрерывность и непрерывную дифференцируемость первых и вторых производных. Мы будем использовать для задания сплайнов другой подход.

Введем обозначение$$S'_N(x_n)=m_n,\quad n=0,1,\dots,N.$$ Запишем сплайн $$S_N$$ в следующей форме$$S_N(x)=f_n(1-t)^2(1+2t)+f_{n+1}t^2(3-2t)+m_nh_nt(1-t)^2-m_{n+1}h_nt^2(1-t),$$ где $$t=\frac{(x-x_n)}{h_n}$$, $$h_n=x_{n+1}-x_n$$.

Условия непрерывности второй производной сплайна определяют систему линейных уравнений$$\lambda_nm_{n-1}+2m_n+\mu_nm_{n+1}= 3\left(\mu_n\frac{f_{n+1}-f_n}{h_n}+\lambda_n\frac{f_n-f_{n-1}}{h_{n-1}}\right),$$ где$$\mu_n=\frac{h_{n-1}}{h_{n-1}+h_n},\quad \lambda_n=1-\mu_n.$$ К этим условиям следует добавить краевые условия. В итоге для мы имеем следующую систему уравнений$$\begin{array}{c} 2m_0+\mu_0m_1=d_0 \\ \\ \lambda_nm_{n-1}+2m_n+\mu_nm_{n+1}=d_n,\ n=1,\dots,N-1 \\ \\ \lambda_Nm_{N-1}+2m_N=d_N \\ \end{array}$$ где$$d_n=3\left(\mu_n\frac{f_{n+1}-f_n}{h_n}+\lambda_n\frac{f_n-f_{n-1}}{h_{n-1}}\right), \ n=1,\dots,N-1.$$ Для условий первого типа $$\mu_0=\lambda_N=0,\quad d_0=2m_0,\quad d_N=2m_N$$, а для условий второго типа$$\mu_0=\lambda_N=1,\quad c_0=3\frac{f_1-f_0}{h_0}-\frac{h_0}{2}M_0,\ c_N=3\frac{f_N-f_{N_1}}{h_{N-1}}+\frac{h_{N-1}}{2}M_N,$$

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

Приведем реализацию нашего класса

$$\begin{verbatim} class TSpline { double[] Xn; int N; TRFunc f; TProgon Progon; double[] m; public TSpline(double[] Xn, TRFunc f) { this.N = Xn.Count(); this.Xn = Xn; this.f = f; } double h(int n) { return Xn[n+1] - Xn[n]; } double mu(int n) { return h(n - 1) / (h(n - 1) + h(n)); } double la(int n) { return 1 - mu(n); } \end{verbatim}$$ $$\begin{verbatim} public void Build(double m0, double mN) { Progon = new TProgon(); Progon.a = new double[N]; Progon.b = new double[N]; Progon.c = new double[N]; Progon.d = new double[N]; int n; for (n = 1; n < N-1; n++) { Progon.a[n] = 2.0; Progon.b[n] = la(n); Progon.c[n] = mu(n); Progon.d[n] = 3.0 * (mu(n) * (f.CalcY(Xn[n + 1]) - f.CalcY(Xn[n])) / h(n) + la(n) * (f.CalcY(Xn[n]) - f.CalcY(Xn[n - 1])) / h(n - 1)); } Progon.a[0] = 2; Progon.a[N-1] = 2; Progon.b[0] = 0; Progon.b[N-1] = 0; Progon.c[0] = 0; Progon.c[N-1] = 0; Progon.d[0] = 2 * m0; Progon.d[N-1] = 2 * mN; Progon.CalcZ(); } \end{verbatim}$$ $$\begin{verbatim} public double CalcSpline(double x) { int n; for (n = 0; n < N-2; n++) { if ((x >= Xn[n]) (x < Xn[n + 1])) { break; } } double t = (x - Xn[n]) / h(n); return f.CalcY(Xn[n]) * (1 - t) * (1 - t) * (1 + 2 * t) + f.CalcY(Xn[n + 1]) * t * t * (3 - 2 * t) + Progon.z[n] * h(n) * t * (1 - t) * (1 - t) - Progon.z[n + 1] * h(n) * t * t * (1 - t); } } \end{verbatim}$$

Теперь испытаем наш сплайн на функции $$\sin x$$. Мы будем строить сплайн всего по шести точкам.

$$\begin{verbatim} TSinFunc F = new TSinFunc(0, 2 * Math.PI); int n; int N = 5; double[] Xn = new double[N + 1]; double dx = (F.Get_b() - F.Get_a()) / (double)N; for (n = 0; n <= N; n++) { Xn[n] = F.Get_a() + (double)n * dx; } TSpline Spline = new TSpline(Xn, F); Spline.Build(1, 1); \end{verbatim}$$ $$\begin{verbatim} StreamWriter Fout = File.CreateText("spline.txt"); double x = F.Get_a(); dx = (F.Get_b() - F.Get_a()) / (10.0 * (double)N); while (x <= F.Get_b()) { Fout.WriteLine("{0}\t{1}\t{2}\t{3}", x, F.CalcY(x), Spline.CalcSpline(x), F.CalcY(x) - Spline.CalcSpline(x)); x += dx; } Fout.Close(); \end{verbatim}$$

На рисунке 15.1 мы приводим погрешность нашей интерполяции. Наибольшая погрешность не превосходит $$0.01$$, что оказывается очень неплохим результатом.

(рис 15.1) Погрешность кубического сплайна

Ключевые термины

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

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

Краткие итоги: Рассмотрены методы интерполяции функций на основе кубических сплайнов. Приведена объектно-ориентированная реализация построение кубического сплайна.

Практическое занятие "Аппроксимация функций"

Цель занятия

Практическое проведение интерполяции для конкретных функций и изучение свойств интерполяционных функций.

Практическая задача

В настоящем занятии мы разберем пример построения интерполяционного многочлена в форме Лагранжа. Мы рассмотрим простейшую функцию$$f(x)=\sin x$$ Возьмем известные значения этой функции в точках$$x_0=0,\ x_1=\frac{\pi}{6},\ x_2=\frac{\pi}{4},\ x_3=\frac{\pi}{3},\ x_4=\frac{\pi}{2}.$$ В этих точках наша функция, как известно, принимает следующие значения$$f_0=0,\ f_1=\frac{1}{2},\ f_2=\frac{1}{\sqrt{2}},\ f_3=\frac{\sqrt{3}}{2},\ f_4=1.$$ Построим интерполяционный многочлен в форме Лагранжа по этим точкам. Согласно известным формулам этот многочлен имеет следующий вид$$P_4(x)=\sum\limits_{i=0}^4p^i_4(x)f_i,$$ где$$p^0_4(x)=\frac{(x-x_1)(x-x_2)(x-x_3)(x-x_4)} {(x_0-x_1)(x_0-x_2)(x_0-x_3)(x_0-x_4)},$$ $$p^1_4(x)=\frac{(x-x_0)(x-x_2)(x-x_3)(x-x_4)} {(x_1-x_0)(x_1-x_2)(x_1-x_3)(x_1-x_4)},$$ $$p^2_4(x)=\frac{(x-x_0)(x-x_1)(x-x_3)(x-x_4)} {(x_2-x_0)(x_2-x_1)(x_2-x_3)(x_2-x_4)},$$ $$p^3_4(x)=\frac{(x-x_0)(x-x_1)(x-x_2)(x-x_4)} {(x_3-x_0)(x_3-x_1)(x_3-x_2)(x_3-x_4)},$$ $$p^4_4(x)=\frac{(x-x_0)(x-x_1)(x-x_2)(x-x_3)} {(x_4-x_0)(x_4-x_1)(x_4-x_2)(x_4-x_3)}.$$

Для использования этого многочлена, разумеется, нет необходимости раскрывать скобки. Для работы с этим многочленом мы напишем небольшой код на языке C#.

$$\begin{verbatim} double pi = Math.PI; int N = 4; double[] xn = { 0, pi / 6.0, pi / 4.0, pi / 3.0, pi / 2.0 }; double[] fn = { 0, 1.0 / 2.0, 1.0 / Math.Sqrt(2.0), Math.Sqrt(3.0) / 2.0, 1 }; double[] hn = new double[N + 1]; int i, j; for (i = 0; i <= N; i++) { hn[i] = 1.0; for (j = 0; j <= N; j++) { if (i == j) { continue; } hn[i] *= xn[i] - xn[j]; } } \end{verbatim}$$ $$\begin{verbatim} double x = pi / 5.0; double res = 0; for (i = 0; i <= N; i++) { double y; y = fn[i] / hn[i]; for (j = 0; j <= N; j++) { if (i == j) { continue; } y *= (x - xn[j]); } res += y; } Console.WriteLine("sin(pi/5) = {0}; Delta = {1}", res, Math.Abs(res - Math.Sin(x))); \end{verbatim}$$

В этой программе мы вычисляем значение нашего интерполяционного многочлена в промежуточной точке $$x=\frac{\pi}{5}$$ и сравниваем с "точным" значением $$\sin \frac{\pi}{5}$$. После запуска мы получим следующий результат.

$$sin(pi/5) = 0.587809514030551; Delta = 2.42617380783461E-05 $$

Можно считать, что это очень не плохой результат для интерполяции функции по столь малому количеству точек.

Лабораторная работа "Аппроксимация функций"

Цель занятия

Практическое проведение интерполяции для конкретных функций и изучение свойств интерполяционных функций.

Сценарий лабораторной работы

Напишите программу на языке C# для построения и вычисления интерполяционного многочлена в форме Лагранжа для произвольных функций и точек разбиения.

Проведите вычислительные эксперименты на примере тригонометрических функций. В этих экспериментах рассчитайте погрешность интерполяционной функции.

Выясните зависимость погрешности интерполяции в зависимости от использования равномерного и неравномерного разбиения отрезка.

Указания

Если на отрезке $$[a,b]$$ задана некоторая числовая функция $$f(x)$$, то разбиением отрезка $$[a,b]$$ называется конечное множество точек $$\{x_n\}_{n=0}^N$$ таких, что$$0=x_0<x_1<\dots<x_N=b.$$ Если также на ряду с разбиением отрезка нам дан набор чисел $$\{f_n\}_{n=0}^N$$, который имеет смысл значений функции в узловых точках$$f(x_n)=f_n,\quad n=0,\dots,N,$$ то возникает задача о построении интерполяционного многочлена.

Интерполяционный многочлен в форме Лагранжа может быть найден по формуле$$P_N(x)=\sum\limits_{i=0}^Np_N^i(x)f_i,$$ где$$p_N^i(x)=\frac{(x-x_0)\cdots(x-x_{i-1})(x-x_{i+1})\cdots(x-x_N)} {(x_i-x_0)\cdots(x_i-x_{i-1})(x_i-x_{i+1})\cdots(x_i-x_N)}.$$