В этом разделе мы разрабатываем алгоритмы для вычисления значений базисных функций и их производных. Пусть `U={u_0,…,u_m}` будет векторным узлом как в уравнении (2.13), и предположим, что мы заинтересованы в базисных функциях степени `p`. Кроме того, предположим, u фиксирован, и `u∈[u_i,u_{i+1})`. Мы разрабатываем пять алгоритмов, которые вычисляют:
Мы представляем два алгоритма, которые вычисляют `p+1` функций до двух, вычислияющих только один, потому что они наиболее важны и на самом деле несколько проще.
Из Св.2.2 и в предположении, что `u∈[u_i,u_{i+1})`, следует, что мы можем сосредоточить наше внимание на функциях `N_{i-p,p},…,N_{i,p}` и их производных; Все остальные функции тождественно равны нулю, и это расточительно, чтобы фактически вычислить их. Таким образом, первым шагом в оценке является определение продолжительности узла, в котором лежит `u`. Можно использовать либо линейный или двоичный поиск векторного узла; мы представляем здесь двоичный поиск. Так как мы используем интервалы вида и `u∈[u_i,u_{i+1})`, тонкая проблема в оценке базисных функций частным случаем `u=u_m`. Лучше всего с этим справиться на самом низком уровне, установив индекс интервала `n(=m-p-1)`. Таким образом, в этом случае `u∈(u_{m-p-1},` `u_{m-p}]`. FindSpan это функция возвращающая целое число, которая означает индекс диапазона.
Int FindSpan(n,p,u,U)
{ /* Определение индекса узлового промежутка */
/* Вход: n,p,u,U */
/* Выход: индекс узлового промежутка */
If (u == U[n + 1]) return(n); /* Особый случай */
low = p; high = n + 1; /* Выполняем двоичный поиск */
mid = (low + high)/2;
while (u < U[mid] || u >= U[mid + 1])
{
If (u < U[mid]) high = mid;
else low = mid;
mid = (low + high)/2;
}
return(mid);
}
Теперь мы укажем второй алгоритм. Предполагая, что `u` в `i`том интервале, вычисление ненулевых функций, приводит к перевернутой треугольной схеме
| `N_{i-p,p}` | ||
| `N_{i-1,1}` | ||
| `N_{i,0}` | `ldots` | `vdots` |
| `N_{i,1}` | ||
| `N_{i,p}` |
Примеры
| Пример2.3 | Пусть `p=2`, `U={0,0,1,2,3,4,4,5,5,}`, и `u=5⁄2` (смотри рисунок
2.6). Тогда `i=4`, так `u∈[u_4,u_5)`. Таким образом, мы вычисляем
|
Это будет ясно читателю, который выполнял замены в этом примере, что есть много избыточных вычислений присущие уравнению (2.5). Например, выписывая функции второй степени в общих чертах, мы имеем
`N_{i-2,2}(u)={u-u_{i-2}}/{u_i-u_{i-2}}N_{i-2,1}(u)+{u_{i+1}-u}/{u_{i+1}-u_{i-1}}N_{i-1,1}(u)`(2.14)
`N_{i-1,2}(u)={u-u_{i-1}}/{u_{i+1}-u_{i-1}}N_{i-1,1}(u)+{u_{i+2}-u}/{u_{i+2}-u_i}N_{i,1}(u)`(2.15)
`N_{i,2}(u)={u-u_i}/{u_{i+2}-u_i}N_{i,1}(u)+{u_{i+3}-u}/{u_{i+3}-u_{i+1}}N_{i+1,1}(u)`(2.16)
Обратите внимание, что`(N_{i-1,1}(u))/{u_{i+1}-u_{i-1}}`
которое появляется во втором члене уравнения (2.14) появляется в первом члене уравнения (2.15); аналогичное утверждение справедливо и для второго члена уравнения (2.15) и первый член уравнения (2.16).`"left"[j]=u-u_{i+1-j}` `"right"[j]=u_{i+j}-u`
Уравнения (2.14) - (2.16), то`N_{i-2,2}(u)={"left"[3]}/{"right"[0]+"left"[3]}N_{i-2,1}(u)+{"right"[1]}/{"right"[1]+"left"[2]}N_{i-1,1}(u)`
`N_{i-1,2}(u)={"left"[2]}/{"right"[1]+"left"[2]}N_{i-1,1}(u)+{"right"[2]}/{"right"[2]+"left"[1]}N_{i,1}(u)`
`N_{i,2}(u)={"left"[1]}/{"right"[2]+"left"[1]}N_{i,1}(u)+{"right"[3]}/{"right"[3]+"left"[0]}N_{i+1,1}(u)`
Основываясь на этих наблюдениях, Алгоритм A2.2 вычисляет все отличные от нуля базисных функций и сохраняет их в массиве `N[0],…,N[p]`.
BasisFuns(i,u,p,U,N)
{ /* Вычислить неисчезающее базисную функцию */
/* Вход: i,u,p,U */
/* Выход: N */
N[0] = 1.0;
for (j=1; j<=p; j++)
{
left[j] = u - U[i + 1 - j];
right[j] = U[i + j] - u;
saved = 0.0;
for (r=0; r<j; r++)
{
temp = N[r] / (right[r + 1] + left[j - r]);
N[r] = saved + right[r + 1]*temp;
saved = left[j - r]*temp;
}
N[j] = saved;
}
}
Заметим, что Алгоритм A2.2 не только эффективен, но также гарантирует, что не будет деления на
ноль, что может произойти с непосредственным применением уравнения (2.5).
Теперь к третьему алгоритму; в частности, мы хотим, чтобы вычислить все `N_{r,p}^{(k)}(u)`, для `i-p≤r≤i` и `0≤k≤n`, где `n≤p`. Инспекция уравнения (2.10) показывает, что основные ингредиенты:
Рассматриваемый как двумерный массив размерности `(p+1)×(p+1)`, базисные функции вписываются в верхний треугольник (в том числе по диагонали), и различные узлы находятся к низу треугольника, т.е.
| `N_{i,0}(u)` | `N_{i-1,1}(u)` | `N_{i-2,2}(u)` |
| `u_{i+1}-u_i` | `N_{i,1}(u)` | `N_{i-1,2}(u)` |
| `u_{i+1}-u_{i-1}` | `u_{i+2}-u_i` | `N_{i,2}(u)` |
Примеры
| Пример2.4 | Пусть `p=2`, U`={0,0,0,1,2,3,4,4,5,5,5}`, и `u=5⁄2`. Тогда `u∈[u_4,u_5)`, и массив становится
`a_{1,0}=1/{u_6-u_4}=1/2` `a_{1,1}=-1/{u_7-u_5}=-1` `a_{2,0}=a_{1,0}/{u_5-u_4}=1/2` `a_{2,1}={a_{1,1}-a_{1,0}}/{u_6-u_5}={-1-1/2}/{4-3}=-3/2` `a_{2,2}=a_{1,1}/{u_7-u_6}=1/{4-4}=1/0` `N_{4,2}^{(1)}=2[a_{1,0}N_{4,1}(5/2)+a_{1,1}N_{5,1}(5/2)]` и`N_{4,2}^{(2)}=2[a_{2,0}N_{4,0}(5/2)+a_{2,1}N_{5,0}(5/2)+a_{2,2}N_{6,0}(5/2)]` Теперь `a_{1,1}`, `a_{2,1}`, и `a_{2,2}` все используют различные узлы, которые не находится в массиве, но они умножаются соответственно на `N_{5,1}(5⁄2)`, `N_{5,0}(5⁄2)`, и `N_{6,0}(5⁄2)`, которые также не в массиве. Эти термины определены равными нулю, и мы остались с`N_{4,2}^{(1)}=2a_{1,0}N_{4,1}(5/2)=1/2` `N_{4,2}^{(2)}=2a_{2,0}N_{4,0}(5/2)=1` Чтобы проверить эти значения, напомним, из Раздела 2.2 следует, что `N_{4,2}(u)=1⁄2(u-2)^2` на `u∈[2,3)`. Вычисление `N_{3,2}^{(1)}(5⁄2)`, `N_{3,2}^{(2)}(5⁄2)`, `N_{2,2}^{(1)}(5⁄2)`, и `N_{2,2}^{(2)}(5⁄2)` аналогично. |
Основываясь на этих наблюдениях (и Пример2.4), это не трудно разработать Алгоритм A2.3, который вычисляет ненулевые базисные функции и их производных, вплоть до `n`ой производной (`n≤p`). Выход в двумерном массиве, ders. ders[k][j] является kой производной функции `N_{i-p+j,p}`, где `0≤k≤n` и `0≤j≤p`. Используются два локальных массива:
DersBasisFuns(i,u,p,n,U,ders)
{ /* Вычислить ненулевые базисные функции и их */
/* производные. Первый раздел это изменённый A2.2 */
/* для хранения функций и узлов различия. */
/* Вход: i,u,p,n,U */
/* Выход: ders */
ndu[0][0] = 1.0;
for (j=1; j<=p; j++)
{
left[j] = u – U[i+1-j];
right[j] = U[i+j] – u;
saved = 0.0;
for (r=0; r<j; r++)
{ /* Нижний треугольник */
ndu[j][r] = right[r+1] + left[j-r];
temp = ndu[r][j-1]/ndu[j][r];
/* Верхний треугольник */
ndu[j][r] = saved + right[r+1]*temp;
saved = left[j-r]*temp;
}
ndu[j][j] = save
}
for (j=0; j<=p; j++) /* Загрузим базисные функции */
ders[0][j] = ndu[j][p];
/* В этом разделе вычисляем производные (уравнение [2.9]) */
for (r=0; r<p; r++) /* Цикл по индексу функции */
{
s1 = 0; s2 = 1; /* Альтернативные строки в массиве */
a[0][0] = 1.0;
/* Цикл вычисления k-ой производной */
for (k=1; k<=n; k++)
{
D = 0.0;
Rk = r – k; pk = p – k;
if (r >= k)
{
a[s2][0] = a[s1][0]/ndu[pk+1][rk];
d = a[s2][0]*ndu[rk][pk];
}
if (rk >= -1) j1 = 1;
else j1 = -rk;
if (r - 1 <= pk) j2 = k - 1;
else j2 = p - r;
for (j=j1; j<=j2; j++)
{
a[s2][j] = (a[s1][j] – a[s1][j-1])/ndu[pk+1][rk+j];
d += a[s2][k]*ndu[rk+j][pk];
}
if (r <= pk)
{
a[s2][k] = -a[s1][k-1]/ndu[pk+1][r];
d += a[s2][k]*ndu[r][pk];
}
ders[k][r] = d;
j = s1; s1 = s2; s2 = j; /* Поменять местами строки */
}
}
/* Умножьте через его в определенных факторов */
/* (Уравнение [2.9]) */
r = p;
for (k=1; k<=n; k++)
{
for (j=0; j<=p; j++) derk[k][j] *=r;
r *= (p - k);
}
}
Обратим внимание теперь на последние два алгоритма, а именно вычисления одной базисной функции, `N_(i,p)(u)`, или производные, `N_(i,p)^((k))(u)`, единичной базисной функции. Решения этих задач приводят к треугольным таблицам вида
| `N_{i,0}` | |||
| `N_{i,1}` | |||
| `N_{i+1,0}` | `N_{i,2}` | ||
| `ldots` | `vdots` | `N_{i,p}` | |
| `N_{i+p-1,0}` | `N_{i+p-2,2}` | ||
| `N_{i+p-1,1}` | |||
| `N_{i+p,0}` |
Примеры
| Пример2.5 | Пусть `p=2`, `U={0,0,0,1,2,3,4,4,5,5,5}`, и `u=5⁄2`. Расчет `N_{3,2}(5⁄2)` дает
|
OneBasisFuns(p,m,U,i,u,Nip)
{ /* Вычислить базисную функцию */
/* Вход: p,m,U,i,u */
/* Выход: Nip */
if ((i == 0 && u == U[0]) || /* Частный */
(i == m-p-1 && u == U[m]) /* случай */
{
Nip = 1.0; return;
}
If (u < U[i] || u >= U[i+p+1]) /* Локальное свойство */
{
Nip = -.0; return;
}
for (j=0; j<=p; j++) /* Инициализация функций нулевой степени */
if (u >= U[i+j] && u < U[i+j+1]) N[j] = 1.0;
else N[j] = 0.0;
for (k=1; k<=p; k++) /* Вычислить треугольную таблицу */
{
if (N[0] == 0.0) saved = 0.0;
else saved = ((u-U[i])*N[0])/(U[i+k]-U[i]);
for (j=0; j<=p-k+1; j++)
{
Uleft = U[i+j+1];
Uright = U[i+j+k+1]
If (N[j+i] == 0.0)
{
N[j] = saved; saved = 0.0;
}
else
{
temp = N[j+1]/(Uright-Uleft);
N[j] = saved+(Uright-u)*temp;
saved = (u- Uleft)*temp;
}
}
}
Nip = N[0];
}
Теперь для фиксированной `i`, вычисление производных, `N_{i,p}^{(k)}(u)`, для `k=0,…,n`, `n≤p`, использует уравнение (2.9). Например, если `p=3` и `n=3`, то
`N_{i,3}^{(1)}=3(N_{i,2}/(u_{i+3}-u_i)-N_{i+1,2}/(u_{i+4}-u_{i+1}))`
`N_{i,3}^{(2)}=3(N_{i,2}^{(1)}/(u_{i+3}-u_i)-N_{i+1,2}^{(1)}/(u_{i+4}-u_{i+1}))`
`N_{i,3}^{(3)}=3(N_{i,2}^{(2)}/(u_{i+3}-u_i)-N_{i+1,2}^{(2)}/(u_{i+4}-u_{i+1}))`
Используя треугольные таблицы, мы должны вычислить| `k=0:` | |||||||||||||||||
|
|||||||||||||||||
| `k=1:` | |||||||||||||||||
|
|||||||||||||||||
| `k=2:` | |||||||||||||||||
|
|||||||||||||||||
| `k=3:` | |||||||||||||||||
|
На словах алгоритм:
Алгоритм A2.5 вычисляет `N_{i,p}^{(k)}(u)` для `k=0,…,n`,`n≤p`. `k`-ая производная возвращается в ders[k].
DersOneBasisFuns(p,m,U,i,u,n,ders)
{ /* Вычислить производных базисной функции Nip */
/* Вход: p,m,U,i,u,n */
/* Выход: ders */
if (u < U[i] || u >= U[i+p+1]) /* Локальное свойство */
{
for (k=0; k<=n; k++) ders[k] = 0.0;
return;
}
for (j=0; j<=p; j++) /* Инициализация функций нулевой степени */
if (u >= U[i+j] && u < U[i+j+1]) N[j][0] = 1.0;
else N[j][0] = 0.0;
for (k=0; k<=p; k++) /* Вычислить полную треугольную таблицу */
{
if (N[0][k-1] == 0.0) saved = 0.0;
else saved = ((u – U[i])*N[0][k-1])/(U[i+k] - U[i]);
for (j=0; j<p-k+1; j++)
{
Uleft = U[i+j+1];
Uright = U[i+j+k+1];
if (N[j+1][k-1] == 0.0)
{
N[j][k] = saved; saved = 0.0;
}
else
{
temp = N[j+1][k-1]/(Uright-Uleft);
N[j][k] = saved+(Uright-u)*temp;
Saved = (u-Uleft)*temp;
}
}
}
ders[0] = N[0][p]; /* Значение функций */
for (k=1; k<=n; k++) /* Вычисление производных */
{
for (j=0; j<k; j++) /* Загрузка соответствующей колонки */
ND[j] = N[j][p-k];
for (jj=1; jj<=k; jj++) /* Вычислить таблицу ширины k */
{
if (ND[0] == 0.0) saved = 0.0;
else saved = ND[0]/(U[i+p-k+jj]-U[i]);
for (j=0; j<k-jj+1; j++)
{
Uleft = U[i+j+1];
Uright = U[i+j+p+jj+1];
if (ND[j+1] == 0.0)
{
ND[j] = (p-k+jj)*saved; saved = 0.0;
}
else
{
temp = ND[j+1]/(Uright-Uleft);
ND[j] = (p-k+jj)*( saved-temp);
saved =temp;
}
}
ders[k] = ND[0]; /* k-ая производная */
}
}
Наконец, отметим, что Алгоритмы A2.3 и A2.5 вычисляют производные с права, если `u` есть узел. Тем не менее, уравнения (2.5), (2.9), (2.10), и другие в этой главе можно было бы определить с помощью интервалов вида и `u∈(u_i,u_{i+1}]`. Это не изменит Алгоритмы A2.2 через A2.5. Другими словами, производные слева можно найти, просто имея алгоритм использующий интервалы в форме `(u_i,u_{i+1}]`, вместо `[u_i,u_{i+1})`. В предыдущем примере, с `p=2` и `U={0,0,0,1,2,3,4,4,5,5,5}`, если `u=2`, то охватывают `i=3` следует производные от левой, и `i=4` урожайность производные с правой стороны.
Упражнения