This routine is called RK45 or Runge-Kutta-Fehlberg. It involves a
fifth order RK method where the difference between the fourth- and
fifth order evaluations is an estimate of the truncation error for the
classical rk4 procedure. The error estimate is used to change the
increment size. I've programmed it in Matlab that involves some
obscure matrix expressions, so I've converted it to pseudo pov.
This routine seems to perform just as well as Matlab's own ode45,
which I've used in many fluid and thermodynamics problems whitout
problems. Problems can arise with stiff, coupled systems but I don't
think cloths can be characterized as such.

It takes the function f(t,x) for the differential equation and finds
the value of x at t=t+h Input: t (current time), h (step size), x
(current variable value). At the end, t and x are updated to the next
time step.

#macro rk45(t,x,h)
c20=.25
c21=.25
c30=.375
c31=.09375
c32=.28125
c40=12/13
c41=1932/2197
c42=-7200/2197
c43=7296/2197
c51=439/216
c52=-8
c53=439/216
c54=-845/4104
c60=.5
c61=-8/27
c62=2
c63=-3544/2565
c64=1859/4104
c65=-.275
a1=25/216
a3=1408/2565
a4=2197/4104
a5=-.2
b1=16/135
b3=6656/12825
b4=28561/56430
b5=-.18
b6=2/55

k1=h*f(t,x)
k2=h*f(t+c20*h,x+c21*K1)
k3=h*f(t+c30*h,x+c31*k1+c32*k2)
k4=h*f(t+c40*h,x+c41*k1+c42*k2+c43*k3)
k5=h*f(t+h,x+c51*k1+c52*k2+c53*k3+c54*k4)
k6=h*f(t+c60*h,x+c61*k1+c62*k2+c63*k3+c64*k4+c65*k5)

x4=x+a1*k1+a3*k3+a4*k4+a5*k5
x=x+b1*k1+b3*k3+b4*k4+b5*k5+b6*k6
t=t+h
error=abs(x-x4)
#end rk45

The next routine uses rk45 and it's error output to increase or
decrease the time increment if the error goes beyond error_max or
below error_min. If its between the increment is kept unchanged.

sign(h) returns the sign of h. The following macro would do the trick:
#macro sign(a) a^2/(abs(a)*a) #end
I'm not sure how to code the exit_while_loop, but I'm sure it's
possible.

h is the initial time step size, tb is stop time, hmax and hmin are
upper and lower bounds on the  step size and errormin and errormax are
the error bounds. iflag has value 0 or 1: iflag=0 means successful
integration from ta-tb, iflag=1 means maximum number of iterations
reached.

#macro rk45Adaptive(t,x,h)
tb=2
itmax=1e4
errormax=1e-3
errormin=1e-4
hmin=1e-3
hmax=5
delta=1e-5
iflag=1
k=0
while k <= itmax
 k=k+1
 if (abs(h)<hmin) h=sign(h)*hmin end
 if (abs(h)>hmax) h=sign(h)*hmax end
 d=abs(tb-t)
 if (d<=abs(h))
  iflag=0  
  if (d<=delta*max(abs(tb),abs(t)) exit_while_loop
  h=sign(h)*d
 end
 xsave=x
 tsave=t
 rk45(t,x,h)
 if (iflag=0) then exit_while_loop
 if (error<error_min) h=2*h end
 if (error>error_max)
  h=h/2
  x=xsave
  t=tsave
  k=k-1
 end
end
#end rk45Adaptive
-- 
ICQ 74734588

Post a reply to this message