# # Program that, using Runge-Kutta method, # solves mxm system of IVP of ODE # x[i]' = f[i](t,x[1],x[2],...,x[m]), t in [a,b] # x[i](0) = x0[i], # i=1,2,...,m # # input: # f --- a Vector of operators # for example: # > f:=Vector(2); # > f[1]:=(t,x,y)->x-0.1*x*y; # f[2]:=(t,x,y)->-.5*y+.02*x*y; # a, b --- end points of the t interval # x0 --- Vector of initial values # n --- number of points dividing [a,b] # # output: t --- the Vector of nodes dividing [a,b] # x --- the array(1..m), # each entry x[i] is a Vector # x[i][j] is the approximate value # of x[i](t[j]) # RungeKutta := proc( f :: Vector, # vector of RHS functions a :: numeric, # left end of the interval b :: numeric, # right end of the interval x0 :: Vector, # vector of initial values n :: integer # number of points for t ) local i, j, k, l, m, delta, delta2, t, x; m := LinearAlgebra[Dimension](f); delta:=evalf( (b-a)/(n-1) ); delta2:=0.5*delta; t:=Vector([ seq( evalf( a+(i-1)*delta ), i=1..n )]); x:=array([ seq( Vector(n), i=1..m) ] ); k:=array(1..4, [ seq( array(1..m), i=1..4) ] ); for i to m do x[i][1]:=evalf(x0[i]) end do; for i from 1 to n-1 do for j to m do k[1][j]:= evalf( delta*f[j](t[i], seq( x[l][i], l=1..m ) ) ) end do; for j to m do k[2][j]:= evalf( delta*f[j](t[i]+delta2, seq( x[l][i]+0.5*k[1][l], l=1..m ) ) ) end do; for j to m do k[3][j]:= evalf( delta*f[j](t[i]+delta2, seq( x[l][i]+0.5*k[2][l], l=1..m ) ) ) end do; for j to m do k[4][j]:= evalf( delta*f[j](t[i]+delta, seq( x[l][i]+k[3][l], l=1..m ) ) ) end do; for j to m do x[j][i+1]:= evalf( x[j][i] + (k[1][j]+2.0*(k[2][j]+k[3][j])+k[4][j])/6.0 ) end do; end do; print(` End of Runge-Kutta method for mxm IVP-ODE `); return t, x; end;