rkf45b:=proc(neqn,t0,tf,u0,nsteps,abserr,relerr,u) # # Procedure rkf45b computes an ODE solution by the variable # step RK Fehlberg 45 method for a series of points along # the solution by repeatedly calling procedure ssrkf45 for # a single RK Fehlberg 45 step. The truncation error is # estimated along the solution to adjust the integration # step according to a specified error tolerance. # # Argument list # # neqn number of first order ODEs # # t0 initial value of independent variable # # tf final value of independent variable # # u0 initial condition vector of length neqn # # nsteps maximum number number of rkf45 steps # # abserr absolute error tolerance # # relerr relative error tolerance # # u ODE solution vector of length neqn after # nsteps steps # # Type variables local e, h, i, j, tm, t, hmin, nfin1: # # Size arrays e:=array(1..neqn): # # Initial integration step h:=(tf-t0)/2.0: # # Minimum allowable step hmin:=(tf-t0)/nsteps: # # Start integration tm:=t0: # # While independent variable is less than the final # value, continue the integration while(tm <= tf*0.999) do # # If the next step along the solution will go past # the final value of the independent variable, set # the step to the remaining distance to the final # value if(tm+h > tf) then h:=tf-tm: end if: # # Single rkf45 step read "c:\\odelib\\maple\\ode1x1\\ssrkf45.txt": ssrkf45(neqn,tm,u0,h,u,e): t:=tm+h: # # Flag for the end of the integration nfin1:=1: # # Check if any of the ODEs have violated the error # criteria for i from 1 to neqn do if(abs(e[i]) > abs(u[i])*relerr+abserr) then # # Error violation, so integration is not complete. # Reduce integration step because of error violation # and repeat integration from the base point h:=h/2.0: nfin1:=0: break: end if: end do: # # If the current step is less than the minimum # allowable step, set the step to the minimum # allowable value and continue integration from # new base point if(h < hmin) then h:=hmin: nfin1:=1: end if: # # If there is no error violation, continue the integration # from new base point if(nfin1 = 1) then for i from 1 to neqn do u0[i]:=u[i]: end do: tm:=t: # # Test if integration step can be increased for i from 1 to neqn do if(abs(e[i]) > (abs(u[i])*relerr+abserr)/32.0) then # # Integration step cannot be increased nfin1:=0: break: end if: end do: # # Increase integration step if(nfin1 = 1) then h:=h*2.0: end if: # # Continue for no error violation (nfin1=1) end if: # # Continue while end do: # # End of rkf45b end: