Форум » Программистские штучки » Средства Матлаба в сплайн-картировании » Ответить
Средства Матлаба в сплайн-картировании
administrator: Начну новую тему,т.к. в теме "ку" разбираются несколько иные вопросы. Сплайны - они и в африке сплайны, метод конечных элементов также и в африке метод конечных элементов, и я верю, что весь набор программ Андрея Михайловича можно реализовать средствами Матлаба без написания дополнительных программ, надо только очень хорошо ориентироваться в его возможностях и одновременно понимать, что за математические объекты и методы вычислений мы на самом деле используем. Сегодня речь пойдет о построении одномерной матрицы A. По существу, это вычисление значений B-сплайна третьей степени. Бэ, а не бета. От слова базисный. Андрей Михайлович ошибался. B-сплайн может быть любой степени, т.к. он абстрактно выводится например как преобразование Фурье функции (sin x/x)^p, где p - нужная степень B-сплайна+1. (см. Марчук, Агошков. Введение в проекционно-сеточные методы. 1981). В матлабе, в spline toolbox, есть функция spcol, решающая задачу, аналогичную программам a0, a0_grd и их векторизованым вариантам (a0rc, a0_grdrc). Только spcol более универсальна и позволяет работать с нерегулярными сетками и с B-сплайнами любой степени, не только третьей. По матлабовской справке spcol генерирует "матрицу коллокации". О связи сплайн-коллокации и нашей возни я уже писал в "Когда просто хочется что-то сказать". Отличия spcol и a0 сводятся к следующему: -a0 работает с "областью построения", а spcol - с координатами узлов, возможна и нерегулярная сетка. Поэтому нужно задавать координату каждого узла, т.е. подовать на вход программы вектор координат узлов. -spcol может построить сплайн любой степени (средний параметр вызова). Для того, чтобы строился наш любимый кубический B-сплайн, нужно указать 4 (на единицу больше). -вектор, генерируемый spcol (или матрица, если задано несколько точек), не содержит расширенной области. Если я укажу n узлов, программа построит вектор длиной n-m, где m - средний параметр (степень сплайна +1). a0 же, наоборот, построит вектор длиной n+2. Суммарная разница составит шесть шагов сетки для кубического сплайна. -при увеличении степени b-сплайна в spcol он ползет (удлиняется) влево, т.е. например если для сплайна 1-й степени ненулевые значения были в 7-м и 6-м узле, то для сплайна 2-й степени - в 5-м, 6, 7-м, а для 3-й степени - в 4,5,6,7. Т.е. целая часть координаты точки отражает положение крайнего _правого_ полинома. В a0 наоборот целая часть координаты точки отражает положение крайнего _левого_ полинома. Поэтому для получения одинаковых результатов нужно к координатам узлов, прииемлемых для a0, прибавить m-1 единиц шага сетки. Окончательно имеем: вызову [pre2]a0rc([a b],x,h)[/pre2] где x может быть как одно число, так и вектор чисел, соответствует: [pre2]spcol([a:h:b+6*h],4,x+3*h])[/pre2] здесь 4-ка определяет степень сплайна+1, а цифры 3 и 6 возникают из-за различия в нумерации узлов в a0 и spcol и того, что spcol не делает расширения области. В более общем виде, для сплайна степени m-1 можно записать: [pre2]spcol([a:h:b+(m+2)*h],m,x+(m-1)*h)[/pre2] но это выражение верно только тогда, когда расширенная область для других степеней задается так же, как для третьей (как на самом деле - не знаю). Вызову a0rc_grd, генерирующую массив значений сплайнов без нулей при выполнении визуализации, соответствует: [pre2]spcol([1:1:8],4,[0:(1/dd):1-(1/dd)]+4),[/pre2] где dd - коэффициент сгущения. Или для сплайна произвольной степени m-1 [pre2]spocl([1:1:2m],m,[0:(1/dd):1-(1/dd)]+m) [/pre2] Ту же самую программу spcol, естественно, можно применить и для B-сплайнов первой степени (функций-крышек). Средний параметр будет 2.
Ответов - 2
administrator: Это были значения сплайнов, но той же программой spcol можно вычислить и производные, применив прочитанное у К. де Бора (К. де Бор. Практическое руководство по сплайнам. Радио и связь, 1985) замечательное свойство B-сплайнов: n-я производная от B-сплайна m-й степени получается как линейная комбинация n+1 B-сплайнов (m-n)-й степени. Если B-сплайн (исходный) был вычислен от точки с коодинатой x, то сплайны линейной комбинации берутся от точек x, x+h,.. x+nh. Например: 2-я производная кубического B-сплайна - это линейная комбинация трех B-сплайнов первой степени от точки x, x+h, x+2h. Не прочитал о коэффициентах линейной комбинации, но подобрал их сам. Для 2-й производной коэффициенты линейной комбинации: 1 -2 1. 1-я производная - комбинация двух параболических сплайнов. Для нее коэффициенты линейной комбинации -1 1 от точек x и x+h. Коэффициенты я нашел довольно просто. Крайние полиномы в производных от B-сплайна третьей степени - не линейная комбинация, а крайний полином от B-сплайна меньшей степени. Сравнивая их с со сплайном меньшей степени от того же числа, определяем знак крайних полиномов. Минус два определил по вычислению a2 от целого числа. Одномерные сплайны могли быть только [1 0], а ответ получился [1 -2 1]. Тогда вырисовывается такой алгоритм замены программ a1 и a2: а) для первой производной, вызову [pre2]sp1=a1([a b],x,h)[/pre2] соответствует последовательность: [pre2] sp11=spcol([a:h:b+5*h],3,x+2*h) sp12=spcol([a:h:b+5*h],3,x+3*h) sp1=sp12-sp21[/pre2] б)для второй производной, вызову [pre2]sp2=a2([a b],x,h)[/pre2] соответствует последовательность: [pre2] sp21=spcol([a:h:b+4*h],2,x+h) sp22=spcol([a:h:b+4*h],2,x+2*h) sp23=spcol([a:h:b+4*h],2,x+3*h) sp2=sp21+sp23-2*sp22. [/pre2] Этот алгоритм сразу векторизованный. Однако он имеет недостаток - перерасход памяти по сравнению с традиционным расчетом, для первой производной - в два раза,для второй - в три. Так же, как и реализации составления матриц Q через произведение Кронекера. Кроме того, лишнее сложение малых чисел (в отличие от прямого вычисления) не увеличивает точность, а совсем наоборот. Отмечу, что spcol и без таких ухищрений может вычислять производные от B-сплайна, но эта программа предназначена для составления матрицы коллокации, и производные она вычисляет для удобства коллокации, что для нас неудобно, а именно: для одной координаты в выходе spcol будут идти подряд несколько строк: значение B-сплайна и его производные до указанной. То есть если использовать spcol в таком варианте, нужно будет фильтровать ее выход, брать каждую вторую или третью строку матрицы. Экономии памяти по сравнению с моим вариантом нет никакой, возможно есть увеличение точности и уменьшение времени счета. Чтобы получить производные этим способом, нужно вектор координат (третий параметр spcol) прогнать через программу brk2knt, у нее 2 параметра, второй - до какой производной нужно проводить расчеты+1 (если до второй, то параметр =3). Т.е. например [pre2]spcol([a:h:b+(m+2)*h],m,brk2knt(x+(m-1)*h,3))[/pre2]
administrator: Кроме того оказалось, что в Матлабе таки реализована поддержка расширенной области. Для этого существует функция augknt. Через нее нужно прогонять вектор узлов сетки. С ее учетом вызову [pre2]a0([a b],h,x)[/pre2] соответствует [pre2]spcol(augknt(a:h:b,m),m,x)[/pre2] где m - степень B-сплайна +1. Но соответствует не совсем, т.к. для точек на краях диапазона (от a до a+(2*h) и от b-(2*h) до b) результат получается другим. Не знаю, где правильно - у нас или через augknt. Возможно, это и не то немного. Но для точек в середине диапазона все сходится.
полная версия страницы