Runge–Kutta–Fehlberg RKF System of ODEs Working CAS Function
11-04-2022, 04:35 AM (This post was last modified: 11-04-2022 05:12 AM by Ioncubekh.)
Post: #1
 Ioncubekh Junior Member Posts: 42 Joined: Oct 2022
Runge–Kutta–Fehlberg RKF System of ODEs Working CAS Function
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

Attached File(s) Thumbnail(s)

HP Prime G1
Python via Android Termux...
 « Next Oldest | Next Newest »

User(s) browsing this thread: 1 Guest(s)