kill(all)$

multiple odes rk4 ,version 1.0b (x[0],x[1],x[2],....)
peter.vlasschaert@gmail.com,16/06/2020

ratprint :false$

load packages 'lratsubst'
----------------------------------

load(lrats)$

definition : function fi(x0,xi), i=1..m
-----------------------------------------------
q : number of the ode (= function of the : rhs)
pe: list of odes
pl  : list variables from odes
pv : list values for evaluation
r   : result for function

f(q,pe,pl,pv,r):=block([],
    s1:pe[q],
    s2:map("=",pl,pv),
    r:ev(lratsubst(s2,s1),nouns)
)$

→ optimization of the code

grind(%);

differential equations
================

rem : ode : (o)rdinary (d)ifferential (e)quations.

differential equations
================

q1:'diff(x[1],t)=x[1]-2*x[2];
q2:'diff(x[2],t)=x[2];

→1e) step : list rhs odes
----------------------------------

p1:[rhs(q1),rhs(q2)];

→2e) step : list rhs odes
----------------------------------

p1a:listofvars(p1);

p1b:append([x[0]],p1a);

p1b[2]=1;

x1;

algo : RK4
========

→1e) step : initialisation
----------------------------------

d1:makelist(0,i,1,length(p1));
d2:makelist(0,i,1,length(p1));
d3:makelist(0,i,1,length(p1));
d4:makelist(0,i,1,length(p1));
d5:makelist(0,i,1,length(p1));

→2e) step : time
----------------------------------

dt:0.1;

→2e) initial values: p1b
----------------------------------

stv:[0,1,0];

algo 

xi:makelist(0,i,1,length(p1));

→3e) matrix :output
---------------------------------

;

x[0]:0;

x[1]:1;
x[2]:0;

xx:matrix([x[0],x[1],x[2]]);

for k:1 thru 6 do
(
for j:1 thru length(p1) do(
    ma1:[x[0]],
    ma2:makelist(x[j],j,1,length(p1)),
    ma3:append(ma1,ma2),
d1[j]:f(j,p1,p1b,ma3,r)*dt  
),
 for j:1 thru length(p1) do(
  a1:  [x[0]+dt/2],
  a2:  makelist(x[j]+d1[j]/2,j,1,length(p1)),
  a3: append(a1,a2),
d2[j]:f(j,p1,p1b,a3,r)*dt
    ),

for j:1 thru length(p1) do
(
  b1:  [x[0]+dt/2],
  b2:  makelist(x[j]+d2[j]/2,j,1,length(p1)),
  b3: append(b1,b2),
d3[j]:f(j,p1,p1b,b3,r)*dt
      ),
for j:1 thru length(p1) do
(
  c1:  [x[0]+dt],
  c2:  makelist(x[j]+d2[j]/2,j,1,length(p1)),
  c3: append(c1,c2),
d4[j]:f(j,p1,p1b,c3,r)*dt
    ),    
for jj:1 thru length(p1) do
   (
    d5[jj]:ev(p1a[jj]+(d1[jj]+2*d2[jj]+2*d3[jj]+d4[jj])/6,nouns) , 
       makelist(x[jj]:d5[jj],jj,1,length(p1))
  ),    
    d6:[x[0]+k*dt,x[1],x[2]],   
   xx:addrow(xx,d6)

)$

xx;