Кубический сплайн

Пусть некоторая функция f(x) задана на отрезке $[a,b]$, разбитом на части $[x_{i-1},x_{i}]$, $a=x_{0}<x_{1}<...<x_{N}=b$. Кубическим сплайном называется функция S(x), которая:

• на каждом отрезке $[x_{i-1},x_{i}]$ является полиномом третьей степени;
• имеет непрерывную вторую производную на всём отрезке $[a,b]$;
• в точках $x_{i}$ выполняется равенство $S(x_{i})=f(x_{i})$;
• $S''(a)=S''(b)=0$.

По построению сплайн S(x) интерполирует функцию f в точках $x_{i}$.

Теорема: Существует ровно одна функция S(x), удовлетворяющая этим условиям.

Построение

Обозначим: $h_{i}=x_{i}-x_{i-1}$

На каждом отрезке $[x_{i},x_{i+1}]$ функция $S(x)$ есть полином третьей степени $S_{i}(x)$, коэффициенты которого надо определить. Запишем для удобства $S_{i}(x)$ в виде:

$S_{i}(x)=a_{i}+b_{i}(x-x_{i})+{c_{i} \over 2}(x-x_{i})^{2}+{d_{i} \over 6}(x-x_{i})^{3}$

тогда

$S_{i}\left(x_{i}\right)=a_{i},\quad S'_{i}(x_{i})=b_{i},\quad S''_{i}(x_{i})=c_{i}$

Условия непрерывности всех производных до второго порядка включительно записываются в виде
$S_{i}\left(x_{i-1}\right)=S_{i-1}(x_{i-1})$
$S'_{i}\left(x_{i-1}\right)=S'_{i-1}(x_{i-1})$
$S''_{i}\left(x_{i-1}\right)=S''_{i-1}(x_{i-1})$

а условия интерполяции в виде 

$S_{i}\left(x_{i-1}\right)=f(x_{i-1})$

Отсюда получаем формулы для вычисления коэффициентов сплайна:

$a_{i}=f\left(x_{i}\right)$

$h_{i}c_{i-1}+2(h_{i}+h_{i+1})c_{i}+h_{i+1}c_{i+1}=6\left({{f_{i+1}-f_{i}} \over {h_{i+1}}}-{{f_{i}-f_{i-1}} \over {h_{i}}}\right)$

$d_{i}={{c_{i}-c_{i-1}} \over {h_{i}}}$

$b_{i}=-{1 \over 2}h_{i}c_{i}-{1 \over 6}h_{i}^{2}d_{i}+{{f_{i}-f_{i-1}} \over {h_{i}}}$

Если учесть, что $c_{0}=c_{n}=0$, то вычисление с можно провести с помощью метода прогонки для трехдиагональной матрицы.

Реализация на языке Haskell

<source lang="hs"> // Интерполирование функций естественными кубическими сплайнами // Разработчик: Назар Андриенко Email: nuzikprogrammer@gmail.com // Автор - хуй!

Data Structures ^ Algorithms

> module Main where

> import Array > import Char

The Sigma ^ Bracketing datatypes. Sigma is the type of the elements that can be multiplied together; Bracketing is the type of bracketings of lists of elements;

> data Sigma = A | B | C > deriving (Show, Read, Eq, Ord, Enum, Ix)

> data Bracketing = S Sigma | Br Bracketing Bracketing

For convenience, we make Bracketing an instance of the Show ^ Read typeclasses. You don't need to know how this works unless you're interested, ^ you won't be asked to modify this part of the numFactorsogram.

> instance Show Bracketing where > showsPrec n (Br a b) > = ('(':) . showsPrec n a . showsPrec n b . (')':) > showsPrec n (S x) > = showsPrec n x > instance Read Bracketing where > readsPrec n xs = term ("(" ++ filter (not . isSpace) xs ++ ")") > where term ('(':xs) > = [(x, xs') | (x, ')':xs') <- term xs] ++ > [(Br a b, xs) | (a, xs' ) <- term xs, > (b, xs) <- term ('(':xs')] > term ('A':xs) = [(S A, xs)] > term ('B':xs) = [(S B, xs)] > term ('C':xs) = [(S C, xs)] > term _ = []

Some utility functions that you'll need to use to answer the questions.

> multiply A A = B > multiply A B = B > multiply A C = A > multiply B A = C > multiply B B = B > multiply B C = A > multiply C A = A > multiply C B = C > multiply C C = C

> factors :: Sigma -> [(Sigma,Sigma)] > factors A = [(A,C),(B,C),(C,A)] > factors B = [(A,A),(A,B),(B,B)] > factors C = [(B,A),(C,B),(C,C)]

> bracketings :: [Sigma] -> [Bracketing] > bracketings [] = [] > bracketings [z] = [S z] > bracketings w = [ Br ub vb | > (u,v) <- splits w, > ub <- bracketings u, vb <- bracketings v ]

> evaluate :: Bracketing -> Sigma > evaluate (S c) = c > evaluate (Br b1 b2) = multiply (evaluate b1) (evaluate b2)

> splits :: [a] -> [([a], [a])] > splits [x] = [] > splits (x:xs) = ([x], xs):[(x:us, vs) | (us, vs) <- splits xs]

Test cases

> testcases :: Sigma > testcases = [[A], [B], [C], [A,C], [B,A], > [A,B,A], [B,B,A], [A,B,B,A], [A,B,A,C,A,B], [A,B,B,C,A,C,B], > [B,B,C,A,B,C,A], [A,B,B,C,A,B,C], [B,B,B,B,B,B,B,B,B,B], > [A,B,B,B,B,B,B,B,B,B,B,B,B,C], > [A,B,B,B,B,B,B,B,B,B,B,B,B,B,B,B,B,C]]

> testcase :: Int -> [Sigma] > testcase i = if (i <= length testcases && i>0) > then testcases !! (i-1) > else error ("There is no testcase " ++ show i)

> testbbb :: Int -> [Sigma] > testbbb i = take i (repeat B)

> teststring :: String -> [Sigma] > teststring = map char2sigma > where char2sigma 'A' = A > char2sigma 'B' = B > char2sigma 'C' = C > char2sigma x = error ("Character "++ [x] ++ " not in Sigma")

Specification for the num problems

> possible_spec :: Sigma -> [Sigma] -> Bool > possible_spec z w = elem z (map evaluate (bracketings w))

> numbers_spec :: Sigma -> [Sigma] -> Integer > numbers_spec z w = toInteger$length (filter (== z) > (map evaluate (bracketings w)))

> brackets_spec :: Sigma -> [Sigma] -> Maybe Bracketing > brackets_spec z w = if possible_spec z w > then Just (snd$head ( filter ((== z).fst) > (zip (map evaluate wbs) wbs))) > else Nothing > where wbs = bracketings w

Task 1 *****************************************************************

Recursive solution to the first problem

 the base case: we've got all possible facotisations of z and a sing splitting of the string

> processFactors :: [(Sigma,Sigma)]->([Sigma],[Sigma])->Sigma->Bool > processFactors [] splitting z = False > processFactors (pair:pairs) splitting z = if ((possible_rec x u == True) && > (possible_rec y v == True)) > then True > else processFactors pairs splitting z > where u=fst (splitting) > v=snd(splitting) > x=fst (pair) > y=snd (pair)

 we process each pair consecutively:
                                              

> processPairs :: [([Sigma],[Sigma])]->Sigma->Bool > processPairs [] z = False > processPairs (splitting:splittings) z = if ((processFactors pairs splitting z) == True) then True > else processPairs splittings z > where pairs=factors(z)

 we apply ProcessPairs to a string w, where splittings represent all possible ways
 we can split the intial pair into two

> possible_rec :: Sigma -> [Sigma] -> Bool > possible_rec z [] = False > possible_rec z [w] = if (z==w) then True > else False > possible_rec z w = processPairs splittings z > where splittings = splits w

Task 2 *****************************************************************

Runtime analysis

Writing out a reccurence tree for T(n) we'll get (n-1) edges going out from the root (which symbolise the (n-1) possible ways of splitting the initial string). From each node there'll be 2 edges, which symbolise the way we split the string (1/n-1, 2/n-2,...,n-1/1). From each of those nodes we'll get 3 edges which symbolise the 3 ways we could get x or y. Notice there are 2 nodes of each number(symmetric in the tree): 2 nodes of 1, 2 nodes of 2,..., 2 nodes of (n-1). Thus our reccurence will look like this:

T(1)=3 T(n)=2*3T(1) + 2*3T(2) + ... + 2*3T(n-1)

OR

T(1)=3 T(n)=6*(T(1)+T(2)+...+T(n-1))

Notice T(1)=3;T(2)=6T(1)=6*3;T(3)=6T(1)+6T(2)=T(2)+6T(2)=7*T(2)=7*6T(1)=7*6*3;

      T(4)=..=7*7*6*3;...;T(n-1)=6*(7^n-3)*3

Thus, T(n) = 3(1+6+6*7+...+6*7^(n-3))=3 + 3*(6(7^(n-4)-1)/2)= 3 + 9*7^(n-4) - 9=

          = 9*7^(n-4) - 6 = O(7^n)

So the total running time is T(n)=O(7^n)

For testing the running-time on PC, it's better to use such case as possible_rec A [B,B,...,B], rather than the test-cases providing in the initial code as the running time ofthis test-case is quiet close to the worst one.

The results are:

T (possible_rec (testbbb 10)) = 10.36 secs T (possible_rec (testbbb 12)) = 163.68 secs

We can see that the running time of the algorithm grows in exponential time and the exponential factor doesn't go higher than 7.

Task 3 *****************************************************************

Recursive solution to the second problem

 possible_num is just a slight modification of possible_rec
 the key idea is that when we split a string w into u and v
 we have to multiply the number of ways we can factorize u to x
 and the number of ways we can factorize v to y to get the proper result

> possible_num :: Sigma -> [Sigma] -> Int > possible_num z [] = 0 > possible_num z [w] = if (z==w) then 1 > else 0 > possible_num z w = numPairs splittings z > where splittings = splits w

> numPairs :: [([Sigma],[Sigma])]->Sigma->Int > numPairs [] z = 0 > numPairs (splitting:splittings) z = numFactors pairs splitting z + numPairs splittings z > where pairs=factors(z)

> numFactors :: [(Sigma,Sigma)]->([Sigma],[Sigma])->Sigma->Int > numFactors [] splitting z = 0 > numFactors (pair:pairs) splitting z = if ((n1 > 0) && (n2 > 0)) > then (n1 * n2) + numFactors pairs splitting z > else numFactors pairs splitting z > where u=fst (splitting) > v=snd (splitting) > x=fst (pair) > y=snd (pair) > n1 = possible_num x u > n2 = possible_num y v

Task 4 *****************************************************************

Dynamic programming solution to the first problem

We are only interested in the cases when j>=i

> possible_dp :: Sigma -> [Sigma] -> Bool > possible_dp z w = table3 ! z ! 1 ! n > where n = length w > table3 = array (A,C) [(z, table z) | z <- [A .. C]] > table z = array (1,n) [(i, column z i) | i <- [1..n]] > column z i = array (1,n) [(j, poss z i j) | j <- [1..n]] > rec = array (1,n) [(i,head (drop (i-1) w)) | i<-[1..n]] > poss z i j = if (j>=i) then > if (j==i) then > if (z == rec ! i) then True > else False > else if any (True==) [((table3 ! fst(pair) ! i  ! k == True) && > (table3  ! snd(pair) ! (k+1) ! j == True)) > | pair <- pairs, k<-[i..j-1]] > then True > else False > else False > where pairs = factors z

Task 5 *****************************************************************

Dynamic programming solution to the second problem

numbers_dp is just a slight modification to possible_dp

> numbers_dp :: Sigma -> [Sigma] -> Integer > numbers_dp z w = table3 ! z ! 1 ! n > where n = length w > table3 = array (A,C) [(z, table z) | z <- [A .. C]] > table z = array (1,n) [(i, column z i) | i <- [1..n]] > column z i = array (1,n) [(j, poss z i j) | j <- [1..n]] > rec = array (1,n) [(i,head (drop (i-1) w)) | i<-[1..n]] > poss z i j = if (j>=i) then > if (j==i) then > if (z == rec ! i) then 1 > else 0 > else sum [(table3 ! fst(pair) ! i ! k) * (table3 ! snd(pair) ! (k+1) ! j) > | pair <- pairs,k<-[i..j-1]] > else 0 > where pairs = factors z

Task 6 *****************************************************************

Runtime analysis discussion

We have to process, n singleton lists, comparing each to A,B and C; than we have to process (n-1) lists of size 2 ([a1,a2],[a2,a3],...,[a(n-1),a(n)]), looking up if they parenthesise to A,B or C; keeping this process going we finally have to process 2 lists of size (n-1) and than 1 list of size n.

So our reccurence-formula will look like this: T(n)=2*3((n-1)*1+(n-2)*2+(n-3)*3+...3*(n-3)+2*(n-2)+1*(n-1))=

   =n+2n+3n+ ... +(n-1)*n - (1+4+9+...+(n-1)^2)=
   =n(1+2+3+...+(n-1)) - (n*(n-1)*(2n-1))/6=
   =n*(n*(n-1)/2) - (n*(n-1)*(2n-1))/6=
   = (2*n^3-6*n^2-n)/6=O(n^3)  

Thus the total runing time for possible_dp and numbers_dp is T(n)=O(n^3)

The actual running time of numbers_dp is smaller than the running time of possible_dp, because multipliing numbers is much more expensive than evaluating booleans.

This statement can be confirmed with following tests:

The test cases show that T(possible_rec)<T(possible_spec)<T(possible_dp):

T (possible_rec A (testbbb 12)) = 163.68 secs T (possible_spec A (testbbb 12)) = 10.83 secs T (possible_dp A (testbbb 12)) = 0.12 secs

T (possible_num A (testbbb 12)) = 1181.39 secs T (numbers_spec A (testbbb 12)) = 10.87 secs T (numbers_dp A (testbbb 12)) = 0.24 secs

Task 7 (optional) ******************************************************

Dynamic programming solution to the last problem

> brackets_dp :: Sigma -> [Sigma] -> Maybe Bracketing > brackets_dp z w = table3 ! z ! 1 ! n > where n = length w > table3 = array (A,C) [(z, table z) | z <- [A .. C]] > table z = array (1,n) [(i, column z i) | i <- [1..n]] > column z i = array (1,n) [(j, brcks z i j) | j <- [1..n]] > brcks z i j = error "'brcks' not implemented yet"