Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project! See Demo

Close

integration routines: rk4,...

Help
2011-05-31
2014-07-21
  • Hi Sylwester,
    recently I had the impression that the integration method Result = RK4( Y, Dydx, X, H, Derivs  )
    has some troubles… 
    It seems that Result = LSODE( Y, X, H, Derivs    ) is missing.

    Let's do a simple check for RK4: calling the routine rkdemo (codes see below) results in:
    GDL> .compile rkdemo
    % Compiled module: RKDEMO.
    GDL> rkdemo
    % Compiled module: RHSL_FLINT.
    init
           0      -10.000000      -10.000000      -1000.0000       300.00000

    % RHSL_FLINT: Variable is undefined: X
    % Execution halted at: RHSL_FLINT          22 rhsl_flint.pro
    %                      RK4                    <INTERNAL_LIBRARY>
    %                      RKDEMO              23 rkdemo.pro
    %                      $MAIN$         
    GDL>

    Source:

    FUNCTION rhsl_flint, X, Y
    return,
    end

    pro rkdemo

      n=21
      xmin=-10.0
      xmax=10.0
      dx=(xmax-xmin)/(n-1)

      x=DINDGEN(n)*dx+xmin
      y=DINDGEN(n)*0.
      i=0
      y_=-1000.
      xx=x
      yy=y
      dydx = rhsl_flint(Xx,Yy)
      print, 'init'
      print, i, xx, x, y, dydx
      print

      h=dx
      for i=1,n-1 do begin

         dydx = rhsl_flint(Xx,Yy)
         yy=RK4(yy, dydx, Xx, H, 'rhsl_flint', /double)
         y=yy
         xx=x
         print, i, x, y
       
      endfor

      plot, x,y

    end

    Best, Bernhard
    _

     
  • Alain C.
    Alain C.
    2011-06-01

    You are right, we have a problem with RK4 in your specific case.
    (the cases we have in test suite seems to be OK)

    Just to help fixing the bug, I found the problem is located in line 275 (code "math_fun_ng.cpp")  
    BaseGDL* Steptwo = e->Interpreter()->call_fun(static_cast<DSubUD*>(newEnv->GetPro())->GetTree());

    As a temporary solution, you can use the RK4JMG()

    Maybe feeling a bug report would be better ?

    Alain

     
  • Hello,

    I believe the problem here is that rk4 is being supplied with a vector of X and not a scalar. I've reported it at the bug tracker:
    https://sourceforge.net/p/gnudatalanguage/bugs/609/

    Thanks for reporting it!
    Sylwester