2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




Начать новую тему Ответить на тему
 
 глюк IMSL???
Сообщение15.12.2010, 15:08 


28/12/08
74
Помогите пожалуйста.
MS Visual Studio 2005 + Intel Fortran 9.1 + IMSL FNL 5.0
Никак не могу разобраться с IMSL-евской функцией PPVAL. Ведёт себя непредсказуемо.

Пример:
Хочу аппроксимировать функцию f(x)=dsin(Pi/2.d0*x) на промежутке [0,2] с помощью B-spline 4-го порядка. Для этого использую:

BSNAK - Computes the “not-a-knot” spline knot sequence.
BSINT - Computes the spline interpolant, returning the B-spline coefficients.
BSVAL - Evaluates a spline, given its B-spline representation.
В принципе, этого достаточно для построения сплайна. Но так как мне надо будет вычислять сплайн во многих промежуточных точках и важна скорость вычисления, то надо преобразовать сплайн в полином.

BSCPP - Converts a spline in B-spline representation to piecewise polynomial representation.
PPVAL - Evaluates a piecewise polynomial.
Код:
program spline
    use BSINT_INT
    use BSNAK_INT
    use BSCPP_INT
    use PPVAL_INT
    use BSVAL_INT
    implicit none
    real(8) f
    real(8) :: xi1=0.0d0, xi2=2.0d0
    integer, parameter  :: korder=4, ndata=20, ncoef=20, nknot=korder+ncoef
    real(8) :: break(ncoef), ppcoef(korder,ncoef), bscoef(ncoef)
    real(8), dimension(nknot) :: xknot
    integer i,nppcf
    real(8) x, xdata(ndata), fdata(ndata)
    real(8) z1,z2,z3
    real(8), parameter :: Pi=3.141592653589793238d0
   
    f(x)=dsin(Pi/2.d0*x)
   
    do i=1, ndata
      x=xi1+(xi2-xi1)*dfloat(i-1)/dfloat(ndata-1)
      xdata(i) = x
      fdata(i) = f(x)
    enddo
    CALL D_BSNAK (NDATA, XDATA, KORDER, XKNOT)
    CALL D_BSINT (NCOEF, XDATA, FDATA, KORDER, XKNOT, BSCOEF)
    CALL D_BSCPP (KORDER, XKNOT, NCOEF, BSCOEF, NPPCF, BREAK, PPCOEF)
   
!   Первый блок
    I=2*ndata-3
    x = dFLOAT(I)/dFLOAT(NDATA-1)
    z1=f(x)
    z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
    z3=D_PPVAL(x,BREAK,PPCOEF)
    print*,'x=',x
    print*,z1,'  ',z2,'  ',z3
    print*   

!   Второй блок
    x = dFLOAT(2*ndata-3)/dFLOAT(NDATA-1)
    z1=f(x)
    z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
    z3=D_PPVAL(x,BREAK,PPCOEF)   
    print*,'x=',x
    print*,z1,'  ',z2,'  ',z3
    print*

!   Третий блок
    I=2*ndata-3
    x = dFLOAT(I)/dFLOAT(NDATA-1)
    z1=f(x)
    z2=D_BSVAL(x,KORDER,XKNOT,NCOEF,BSCOEF)
    z3=D_PPVAL(x,BREAK,PPCOEF)
    print*,'x=',x
    print*,z1,'  ',z2,'  ',z3
    print*

    pause 'FINISH: press ENTER.'
end program spline

Результат:
Код:
x=1.94736842105263
8.257934547233223E-002     8.258268483722089E-002    [b]0.000000000000000E+000[/b]

x=1.94736842105263
8.257934547233223E-002     8.258268483722089E-002    8.258268483722092E-002

x=1.94736842105263
8.257934547233223E-002     8.258268483722089E-002    8.258268483722092E-002

Вот те жирные нули в первой строке результата выглядят совершенно странно. Должно быть то же самое, что строкой ниже, то-есть второй столбик (рассчитанным по BSVAL) должен практически совпадать с третьим (рассчитанный по PPVAL). У мене IMSL 5.0. Интересно, что покажет IMSL 6.0.

Что интересно, при некоторых x функция PPVAL в каждом блоке считает правильно (например x=1.1234d0)! А при некоторых нет.

Обратите внимание, первый и третий блоки идентичны, а результат расчёта PPVAL в них разный!
При компиляции оптимизация отсутствовала. Было бы хорошо, если бы кто-то запустил этот код у себя и закинул сюда результат. Было бы хорошо знать, исправлен ли этот глюк в последующих версиях фортрана и ИМСЛя. И вообще следует разобраться, чей глюч: фортрана, ИМСЛя или мой.

PS.
Я уже поднял этот вопрос на ру-борде, но ребята там сильно отклонилиь от темы.

 Профиль  
                  
 
 Re: глюк IMSL???
Сообщение24.12.2010, 01:04 


28/12/08
74
Проблема решена. Решение можно найти здесь и здесь.

 Профиль  
                  
 
 Re: глюк IMSL???
Сообщение24.12.2010, 12:35 


28/12/08
74
Коротко.
Следует явно указать опциональные параметры в PPVAL:
Код:
D_PPVAL(x,BREAK,PPCOEF,korder=korder,nintv=nppcf)

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 3 ] 

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group