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