kill(all)$

maxima: numerical solution rkf45, 'example, fourth order ODE',version 1.0
peter.vlasschaert@gmail.com,03/10/2020

load(lrats)$ load(rkf45)$ ratprint:false$

ref:maxima : numerical solution,nonlinear ODEs with rkf45 ,version 1.0a
peter.vlasschaert@gmail.com,03/10/2020.

example : fourth order ODE( → reduction to system first order ODE)
--------------------------------------

given : ODE ,info
==============

p1:'diff(x,t,4)+7.2*'diff(x,t,3)-4*'diff(x,t,2)-5*'diff(x,t,1)-2*x-8=0;

p1a:listofvars(p1);
p1b:derivdegree(p1,p1a[1],p1a[2]);

1e step solve : highest degree of derivative.
------------------------------------------------------------

p2:solve(p1,'diff(p1a[1],p1a[2],p1b))[1];/* must be expr ≠ */

p2a:expand(p2);

2e x1,x2,x3,x4
---------------------

p3:makelist(concat(x,i),i,1,p1b);

p3a:makelist('diff(concat(p1a[1],i),p1a[2])=concat(p1a[1],i+1),i,1,p1b-1);

3e)  put in 'p2' values, x1,x2,x3,x4
-----------------------------------------------

→ rhs 'x1,x2,x3',x4

a)  derivatives,rhs

'diff(p1a[1],t,3) =concat(x,4);

p4:makelist('diff(p1a[1],p1a[2],i) = concat(x,i+1),i,1,p1b-1);

→ rhs(p2)

p4a:lratsubst(p4,rhs(p2));

p4b:lratsubst([p1a[1]=p3[1]],p4a);

→ lhs(p2)

p4c:'diff(concat(p1a[1],p1b),p1a[2]);

→lhs(p2) =rhs(p2)

p4d:p4c=p4b;

4e) 4 ODEs ,first order
--------------------------------

p5:append(p3a,[p4d]);

p5a:expand(p5);

5e) rhs(p5),for rkf45
-----------------------------

p6:map(rhs,p5);

6e)  system of ODEs
------------------------------

-> coefmatrix 'p5'

p7:(-1)*coefmatrix(p5,p3);

p71:(-1)*augcoefmatrix(p5,p3);

p7a:map(lhs,p5);

p7b:apply(matrix,[p7a]);

→ lhs (ODE system)

p7c:transpose(p7b);

→rhs (ODE system)

p7d:apply(matrix,[p3]);

p7e:transpose(p7d);

→ lhs =rhs  ('ODE system')

math : X' = A*X

p7f:p7c = p7.p7e;

rem:p7f : not * =>  dot

rem :
x''-3*x'+2*x =cos(t)
x1=x ,x2=x',x2'=x''

solve :rkf45 ,use(p6)
===============

sol:rkf45(p6,p3,[0.75,1.6,0.5,1.2],[t,0,10],report = true)$

find : intial conditions ODE system
===========================

sol[1];

output :rkf45 , 'wxplot2d'
===================

for i:1 thru p1b do(
q1[i]:[discrete,map(lambda([u],part(u,[1,3])),sol)])$

q2:[q1[1],q1[2],q1[3],q1[4]]$

wxplot2d(q2,
[style,[lines,2]],[xlabel,"t"],[ylabel,"x, x',x'',x''',x''''"],
[legend,"x(t)","x'(t)","x''(t)","x'''(t)"],
[title,"4e order ODE "],[x,0,10],[y,0,9000]);