## integration routines: rk4,... document.SUBSCRIPTION_OPTIONS = { "thing": "thread", "subscribed": false, "url": "subscribe", "icon": { "css": "fa fa-envelope-o" } };

Help
2011-05-31
2014-07-21
• Bernhard Seiwald - 2011-05-31

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. - 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

• Sylwester Arabas - 2014-07-21

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