Working Runge–Kutta–Fehlberg RKF for one equation

[Ioncubekh]

Two things that I would like users help on
1] How to remove crtlst() function; any other alternative so that list of values generated by python (x3 lists) can be directly assigned to global C1, C2 & C3 of 2var-Stats APP
2] Suggestions to use it for 'system of equation' provided no numpy on Prime

I am attaching a PDF. Example one is solved by the provided code & results are matching. So one function f(x,y) at least can be solved by this RKF code

Instructions:

Code:

rkf(f,a,b,y0,tol,hmax,hmin,inHigher,inLower) -> [list(),list()]

  Runge-Kutta-Fehlberg which combines RK 4th and RK 5h
  order for solving an IVP in the form:
        / y'(x) = f(x,y)
        \ y(a) = y0

  RK5 for truncation error:
        y5_i+1 = y5_i + d1*k1 + d3*k3 + d4*k4 + d5*k5 + d6*k6

  RK4 for local error estimation:
        y4_i+1 = y4_i + c1*k1 + c3*k3 + c4*k4 + c5*k5

  Evaluations:
        k1 = h * f( x, y )
        k2 = h * f( x + a2 * h , y + b21 * k1)
        k3 = h * f( x + a3 * h , y + b31 * k1 + b32 * k2)
        k4 = h * f( x + a4 * h , y + b41 * k1 + b42 * k2 + b43 * k3 )
        k5 = h * f( x + a5 * h , y + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4 )
        k6 = h * f( x + a6 * h , y + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5 )


  Params
  ----------
      f        -> (function) to be solved
      a        -> (float) initial interval value
      b        -> (float) end interval value
      y0       -> (float) value of the function at point a
      tol      -> (float) tolerance for truncation error
      hmax     -> (float) max step size
      hmin     -> (float) min step size
      inHigher -> (float) higher increment for step variation
      inLower  -> (float) lower increment for step variation
Return
  -------
      X     -> (list) independent variable
      Y     -> (list) solution value of the function
      P     -> (list) of steps used according to X

RKS algorithm, solving example 1 from attached document

Code:

EXPORT crtlst(dt1,dt2,dt3)
BEGIN
PRINT;                        
C1:={}; C1:=dt1;
C2:={}; C2:=dt2;
C3:={}; C3:=dt3;
END;

EXPORT pyRKF()
BEGIN
PYTHON(rkf);
END;



#PYTHON rkf()
#start

from math import *
import hpprime as hp

# some defaults for arguments provided or user can define latter
def rkf(f,a,b,y0,tol=0.001,hmax=0.2,hmin=0.1,inHigher=0.1,inLower=0.1):

  # Coefficients related to the independent variable of the evaluations
  a2  =   2.500000000000000e-01  #  1/4
  a3  =   3.750000000000000e-01  #  3/8
  a4  =   9.230769230769231e-01  #  12/13
  a5  =   1.000000000000000e+00  #  1
  a6  =   5.000000000000000e-01  #  1/2

  # Coefficients related to the dependent variable of the evaluations
  b21 =   2.500000000000000e-01  #  1/4
  b31 =   9.375000000000000e-02  #  3/32
  b32 =   2.812500000000000e-01  #  9/32
  b41 =   8.793809740555303e-01  #  1932/2197
  b42 =  -3.277196176604461e+00  # -7200/2197
  b43 =   3.320892125625853e+00  #  7296/2197
  b51 =   2.032407407407407e+00  #  439/216
  b52 =  -8.000000000000000e+00  # -8
  b53 =   7.173489278752436e+00  #  3680/513
  b54 =  -2.058966861598441e-01  # -845/4104
  b61 =  -2.962962962962963e-01  # -8/27
  b62 =   2.000000000000000e+00  #  2
  b63 =  -1.381676413255361e+00  # -3544/2565
  b64 =   4.529727095516569e-01  #  1859/4104
  b65 =  -2.750000000000000e-01  # -11/40

  # Coefficients related to the truncation error
  # Obtained through the difference of the 5th and 4th order RK methods:
  #     R = (1/h)|y5_i+1 - y4_i+1|
  r1  =   2.777777777777778e-03  #  1/360
  r3  =  -2.994152046783626e-02  # -128/4275
  r4  =  -2.919989367357789e-02  # -2197/75240
  r5  =   2.000000000000000e-02  #  1/50
  r6  =   3.636363636363636e-02  #  2/55

  # Coefficients related to RK 4th order method
  c1  =   1.157407407407407e-01  #  25/216
  c3  =   5.489278752436647e-01  #  1408/2565
  c4  =   5.353313840155945e-01  #  2197/4104
  c5  =  -2.000000000000000e-01  # -1/5

  # Init x and y with initial values a and y0
  # Init step h with hmax, taking the biggest step possible
  x = a
  y = y0
  h = hmax

  # Init vectors to be returned
  X = [x]
  Y = [y]
  P = [h]

  while x < b - h:

      # Store evaluation values
      k1 = h * f( x, y )
      k2 = h * f( x + a2 * h, y + b21 * k1 )
      k3 = h * f( x + a3 * h, y + b31 * k1 + b32 * k2 )
      k4 = h * f( x + a4 * h, y + b41 * k1 + b42 * k2 + b43 * k3 )
      k5 = h * f( x + a5 * h, y + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4 )
      k6 = h * f( x + a6 * h, y + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5 )

      # Calulate local truncation error
      r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
      # If it is less than the tolerance, the step is accepted and RK4 value is stored
      if r <= tol:
          x = x + h
          y = y + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
          X.append( x )
          Y.append( y )
          P.append( h )

      # Prevent zero division
      if r == 0: r = P[-1]
      
      # Calculate next step size
      h = h * min( max( 0.84 * ( tol / r )**0.25, inLower ), inHigher )

      # Upper limit with hmax and lower with hmin
      h = hmax if h > hmax else hmin if h < hmin else h

  return X,Y,P

def tolst(a):
  return  str(a).replace('[','{').replace(']','}')


#define function
f0= lambda t, u: u*(1-u)*(2-u)
# define arguments
t,f0_,step=rkf(f0, 0, 10, 0.5, 0.01, 0.3, 0.2, 0.1, 0.1 )

#string manipulations
hp.eval('crtlst('+tolst(t)+','+tolst(f0_)+','+tolst(step)+')')
print('Done')

#end

[Ioncubekh]

Ok I think I have been able to come up with a no-Numpy RKF script & it was really an headache taking care of dimensions & data flow in the absence of arrays. But what is good is that the numpy version (as opposed to what i thought) took greater time for execution

Code:

In [176]: %run rkfClassnoNump.py
Execution time: 0.037133216857910156 seconds
Number of data points: 914

In [177]: %run rkfClass.py
Execution time: 0.07887840270996094 seconds
Number of data points: 914

Lorentz system of x3 equations were used as a test drive. Both .py scripts (rename .txt) produce same image. See attachments, tomorrow will test no-Numpy version on Prime.

time() module has to be removed before checking on Prime