rkc4b:=proc(neqn,t0,tf,u0,nsteps,abserr,relerr,u) # # Procedure rkc4b computes an ODE solution by a variable step # classical fourth order RK method for a series of points # along the solution by repeatedly calling procedure ssrkc4 # for a single classical fourth order RK 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 of rkc4 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 rkc4 step read "c:\\odelib\\maple\\pdenon\\ssrkc4.txt": ssrkc4(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)/16.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 rkc4b end: