RK1, RK2, RK3, RK4 one can use this application
https://www.hpmuseum.org/forum/thread-19012.html
RKF for single equation one can use
https://www.hpmuseum.org/forum/thread-19...#pid165807
CAUTION must read CAS Equations Internal Representation For Python:
https://www.hpmuseum.org/forum/thread-19071.html & avoid symbols such as √
Solved Examples from the document:
https://github.com/defencedog/jupyterNot...es_ODE.pdf
Usage: While in CAS, APP Statistics 2Var selected
~
t is reserved for algorithm & is not be used in [Variables] but can be used in [Symbolic Functions]
~
Flush CAS variables by typing reset as well as or
purge(x,y,z) purges x, y, and z or use different variables than previously defined
~ APP Statistics 2Var variables C1 .. C8 will be flushed
~ C1 will have independant variable t, C2 will be the function y, C3 y', C4 y'' ... So one usually plots C1 versus C2
Code:
RKF( [Symbolic Functions], [Variables], [Initial Variable Values], [Start End] )
Compare results with this jupyter notebook:
https://github.com/defencedog/jupyterNot...berg.ipynb
Example 2,3 & 4 attached
Code:
EXPORT RKF(a,b,c,d)
BEGIN
PRINT;
MSGBOX("Select App Statistics 2Var before proceeding");
PRINT("Done in "+TEVAL(algo(a,b,c,d)));
C1:={}; C1:=AVars("C1");
C2:={}; C2:=AVars("C2");
C3:={}; C3:=AVars("C3");
C4:={}; C4:=AVars("C4");
C5:={}; C5:=AVars("C5");
C6:={}; C6:=AVars("C6");
C7:={}; C7:=AVars("C7");
C8:={}; C8:=AVars("C8");
END;
#PYTHON algo(a,b,c,d)
import sys
import hpprime as hp
from math import *
class rkf():
def __init__(self,f, a, b, x0, atol=1e-8, rtol=1e-6, hmax=1e-1, hmin=1e-40):
self.f=f
self.a=a
self.b=b
self.x0=x0
self.atol=atol
self.rtol=rtol
self.hmax=hmax
self.hmin=hmin
def solve(self):
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
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
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
c1 = 1.157407407407407e-01 # 25/216
c3 = 5.489278752436647e-01 # 1408/2565
c4 = 5.353313840155945e-01 # 2197/4104
c5 = -2.000000000000000e-01 # -1/5
t = self.a
x = self.x0
h = self.hmax
T = [t]
X = [x]
while t < self.b:
if t + h > self.b:
h = self.b - t
# k1 = h * self.f(t, x)
k1 = list(map(lambda i:h*i, self.f(t, x)))
# k2 = h * self.f(t + a2 * h, x + b21 * k1 )
k2 = list(map(lambda i:h*i, self.f(t + a2 * h, [sum(i) for i in zip(x, [b21*i for i in k1] )] )))
# k3 = h * self.f(t + a3 * h, x + b31 * k1 + b32 * k2)
k3 = list(map(lambda i:h*i, self.f(t + a3 * h, [sum(i) for i in zip(x, [b31*i for i in k1], [b32*i for i in k2] )] )))
# k4 = h * self.f(t + a4 * h, x + b41 * k1 + b42 * k2 + b43 * k3)
k4 = list(map(lambda i:h*i, self.f(t + a4 * h, [sum(i) for i in zip(x, [b41*i for i in k1], [b42*i for i in k2], [b43*i for i in k3] )] )))
# k5 = h * self.f(t + a5 * h, x + b51 * k1 + b52 * k2 + b53 * k3 + b54 * k4)
k5 = list(map(lambda i:h*i, self.f(t + a5 * h, [sum(i) for i in zip(x, [b51*i for i in k1], [b52*i for i in k2], [b53*i for i in k3], [b54*i for i in k4] )] )))
# k6 = h * self.f(t + a6 * h, x + b61 * k1 + b62 * k2 + b63 * k3 + b64 * k4 + b65 * k5)
k6 = list(map(lambda i:h*i, self.f(t + a6 * h, [sum(i) for i in zip(x, [b61*i for i in k1], [b62*i for i in k2], [b63*i for i in k3], [b64*i for i in k4], [b65*i for i in k5] )] )))
# r = abs( r1 * k1 + r3 * k3 + r4 * k4 + r5 * k5 + r6 * k6 ) / h
r = list(map(lambda i:abs(i)/h, [sum(i) for i in zip([r1*i for i in k1], [r3*i for i in k3], [r4*i for i in k4], [r5*i for i in k5], [r6*i for i in k6] )] ))
# r = r / (self.atol+self.rtol*(abs(x)+abs(k1)))
r = [ (r[i] / ( self.atol+self.rtol*( abs(x[i]) + abs(k1[i]) ) ) ) for i in range(len(x)) ]
if len( r ) > 0:
r = max( r )
if r <= 1:
t = t + h
# x = x + c1 * k1 + c3 * k3 + c4 * k4 + c5 * k5
x = [sum(i) for i in zip(x, [c1*i for i in k1], [c3*i for i in k3], [c4*i for i in k4], [c5*i for i in k5] )]
T.append( t )
X.append( x )
h = h * min( max( 0.94 * ( 1 / r )**0.25, 0.1 ), 4.0 )
if h > self.hmax:
h = self.hmax
elif h < self.hmin or t==t-h:
raise RuntimeError("Error: Could not converge to the required tolerance.")
break
return (T,X)
def strlst(a):
a=a.lower()
if '{' in a:
_a=a.replace('{','').replace('}','').replace('^','**').replace('log','log10').replace('ln','log').replace('-','-').replace(' ','')
else:
_a=a.replace('[','').replace(']','').replace('^','**').replace('log','log10').replace('ln','log').replace('-','-').replace(' ','')
return (_a.split(','))
fns = strlst(sys.argv[0])
var_= strlst(sys.argv[1])
x0 = strlst(sys.argv[2])
sten= strlst(sys.argv[3])
x0 = [float(i) for i in x0]
sten = [float(i) for i in sten]
print('CAS string: ', sys.argv[0])
print('Functions: ', fns)
print('Variables: ', var_)
print('Initial Values: ', x0)
print('Start at: ',sten[0], 'End at: ', sten[1], '\n')
def functs(t,u):
for i in range(len(u)):
_i=var_[i]+'='+str(u[i])
exec(_i)
exec('t='+str(t))
return [eval(i) for i in fns]
t, sol = rkf( f=functs, a=sten[0], b=sten[1], x0=x0, atol=1e-8, rtol=1e-6 , hmax=1e-1, hmin=1e-40).solve()
def tolst(a):
return str(a).replace('[','{').replace(']','}')
for i in range(9):
var = 'C'+str(i)
hp.eval('AVars('+'"'+var+'"'+'):={}')
for i in range(1,len(x0)+2):
var = 'C'+str(i)
print('Column is populated:', var)
if (i==1):
hp.eval('AVars('+'"'+var+'"'+'):='+tolst(t))
else:
val = [val[i-2] for val in sol]
hp.eval('AVars('+'"'+var+'"'+'):='+tolst(val))
#end
Search
Member List
Calendar
Help
/
Runge–Kutta–Fehlberg RKF System of ODEs Working CAS Function
