Современные численные методы в объектно-ориентированном изложении на C#
[+]
Опубликован: 18.05.2011 | Доступ: свободный | Студентов: 968 / 105 | Оценка: 4.40 / 4.20 | Длительность: 12:30:00
Тема: Программирование
Специальности: Программист, Архитектор программного обеспечения
Теги: objective-c, алгоритмы, анализ, вычисления, графика, динамическая система, дифференциальные уравнения, игры, коэффициенты Фурье, кубический сплайн, матричная игра, машина поста, метод Эйлера, методы рунге-кутты, метрическое пространство, моделирование, программирование, процедуры, теория, трансцендентное уравнение, функциональное пространство, элементы
Лекция 16:

Приближение сплайнами
A
|
версия для печати
< Лекция 15 || Лекция 16: 123 || Лекция 17 >
Аннотация: Лекция посвящена вопросам интерполяции числовых функций с помощью сплайнов. Рассматриваются методы построения кубических сплайнов.
Ключевые слова: кубический сплайн, интерполирование функций, числовой функцией, сплайн, производные, математика, разбиение, функция, путь, коэффициенты, матрица, метод прогонки, погрешность

Цель лекции: Описать стандартные методы построения кубических сплайнов. Реализовать на языке 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], что для этой функции выполнены следующие условия:

  1. S_N(x_n)=f_n, n=0,1,\dots,N.
  2. S_N(x), S'_N(x), S''_N(x) являются непрерывными функциями.
  3. На каждом отрезке [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.

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

  1. S'_N(x_0)=m_0, S'_N(x_N)=m_N.
  2. 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), ( 15.1)
где
\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} ( 15.2)
где
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. Погрешность кубического сплайна

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

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

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

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

Дальше >>
< Лекция 15 || Лекция 16: 123 || Лекция 17 >

Вопросы и ответы

вопросов: 0