Форум программистов [Powered by Invision Power Board]

Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Fortran > [Pro] Как задать производную в общем виде

Автор: f999t1 20.12.2008, 01:00

как задать в фортране производную в общем виде
или можно какую-нибудь библиотеку подключить

на самом деле, решаю следующее
есть выражение вида x(j+1)=1-a(x(j))^2

на каждом шаге надо вычислять производную dx(j+1)/dx(j),
но скажым шагом полином увеличивается

подскажите, примерный алгоритм

Автор: FatalError 22.12.2008, 15:02

Вы бы хоть функцию нормально записали, а то возведение в степень там осталось не у дел перед скобкой, да и не ясно, с чего бы это

Цитата

скажым шагом полином увеличивается

Тогда, может, кто чего и подскажет.

Автор: Lido4ka 22.12.2008, 19:50

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

Автор: f999t1 23.12.2008, 22:20

так с производными боле мене понятно
вот текст

Код


program r

    integer i,j,k

    real f ,fb, ri, ri2, xi,xi2,x,r0,xb,rb,e,delta,deltar,fi,fib,fai
      real faib
  
    x=2.0/3.0
    r0=1
      rb=2
      
    e=0.0001
    
      xi=x       
    write(*,*) xi,x   
      
    
    
    !fi2=1-ri*(fi)**2
    ri=r0
      delta=1
    !*****************************************
    do while (delta>e)   ! i=1,70
      
    fi=fa(ri,xi)
    fib=fa(rb,xi)
      write(*,*) fi,fib
    ri2=ri-( (fi-xi)/(fi-fib)  )*(ri-rb)
      
    delta=abs(ri2-ri)
    write(*,*)'delta=',delta,'ri2=',ri2
    ri=ri2
      end do

    
      !**********************************
       

      


    
    
    
    contains
    real function fa(rt,xt)
    real xt,rt

    fa=1-rt*(xt**2)
    end function fa
        

    end

как задать рекурсию?
x(j+1)=1-a( x(j) )^2

т.е в первом ходе вычислил а из уравнения
fa=1-a*(x^2)
во втором ходе надо это сделать из уравнения
fa2=1-a*(fa^2)=1-a+2(a^3)(x^2)-(a^2)(x^4)
и т.д.

Автор: popovda 27.12.2008, 14:37

Производные на фортране сложно записать аналитически, не парсер же прикручивать smile

Воспользуйтесь конечными разностями, как по определению производной:

Код


(f(x + h) - f(x))/h ! Правая конечная разность
(f(x) -  f(x - h))/h ! Левая ----------"-------"-------
0.5*(f(x+h) - f(x-h))/h ! Центральная

Вторая производная 

(f(x+h) -2*f(x) + f(x-h))/h**2


Шаг выбирается в зависимости от типа данных, на Intel x86 это ~10^-3 --- 10^4 для real(4)
h~10-6 --- 10^-8 для real(8)


В общем виде --- да напишите функцию

derivs(f,x) - писать не буду:) Но интерфейс такой

interface
      real(8) function derives(f,x)
                 interface 
                         real(8) function f(x) 
                                    real(8), intent(in)::x;
                         end function f
                 end interface
                 
                 real(8), intent(in)::x
      end function derives
end interface

А рекурсию (хотя лучше её избегать) делают с помощью объявлений функции типа
recurcive function f(x) result A........