ssrkc4:=proc(neqn,t0,u0,h,u,e) # # Procedure ssrkc4 computes an ODE solution by the classical fourth # order RK method for one step along the solution (by calls to derv # to define the ODE derivative vector). It also estimates the # truncation error of the solution, and applies this estimate as # a correction to the solution vector. # # Argument list # # neqn number of first order ODEs # # t0 initial value of independent variable # # u0 initial condition vector of length neqn # # h integration step # # u ODE solution vector of length neqn after # one modified Euler step # # e estimate of truncation error of the solu- # tion vector # # Type variables local ut0, ut, t, i, k1, k2, k3, k4, sum2, sum4: # # Declare arrays ut0:=array(1..neqn): ut:=array(1..neqn): k1:=array(1..neqn): k2:=array(1..neqn): k3:=array(1..neqn): k4:=array(1..neqn): sum2:=array(1..neqn): sum4:=array(1..neqn): # # Derivative vector at initial (base) point read "c:\\odelib\\maple\\pdelin\\derv.txt": derv(neqn,t0,u0,ut0): # # k1, advance of dependent variable vector and # independent variable for calculation of k2 for i from 1 to neqn do k1[i]:=h*ut0[i]: u[i]:=u0[i]+0.5*k1[i]: end do: t:=t0+0.5*h: # # Derivative vector at new u, t derv(neqn,t,u,ut): # # k2, advance of dependent variable vector and # independent variable for calculation of k3 for i from 1 to neqn do k2[i]:=h*ut[i]: u[i]:=u0[i]+0.5*k2[i]: end do: t:=t0+0.5*h: # # Derivative vector at new u, t derv(neqn,t,u,ut): # # k3, advance of dependent variable vector and # independent variable for calculation of k4 for i from 1 to neqn do k3[i]:=h*ut[i]: u[i]:=u0[i]+k3[i]: end do: t:=t0+h: # # Derivative vector at new u, t derv(neqn,t,u,ut): # # k4 for i from 1 to neqn do k4[i]:=h*ut[i]: end do: # # Second order step for i from 1 to neqn do sum2[i]:=u0[i]+k2[i]: end do: # # Fourth order step for i from 1 to neqn do sum4[i]:=u0[i]+(1.0/6.0)*(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i]): end do: t:=t0+h: # # Truncation error estimate for i from 1 to neqn do e[i]:=sum4[i]-sum2[i]: end do: # # Fourth order solution vector (from 2,4 RK pair); # two ways to the same result are listed for i from 1 to neqn do # u[i]:=sum2[i]+e[i]: u[i]:=sum4[i]: end do: # # End of ssrkc4 end: