kill(all)$

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

numerical solution :second nonlinear ODE  'the van der Pol equations'
================================

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

How to understand this code ?
=======================

       equ:'diff(x,t,2)+4*(x^2-1)*'diff(x,t)+x = 0;/* ODE */
       equ2:['diff(x1,t) = x2,
       ev(solve(equ,'diff(x,t,2))[1],'diff(x,t,2) = 'diff(x2,t),
       'diff(x,t) = x2,x = x1)]; /*two first order ODEs */
        equ2:map(rhs,equ2);      
       

analyse
------------

1e) second order ODE : x1,x2,highest derivative 
=====================================

p1:'diff(x,t,2)+4*(x^2-1)*'diff(x,t)+x = 0;

p1aa:listofvars(p1);

p1a:solve(p1,'diff(x,t,2));

p1b:p1a[1];

rem :third order ODE : x1,x2,x3,...............'diff(x,t,3) = highest derivatives

2e)  second order ODE: 'numerical solution '
==================================

need : two first order equation (to use rkf45)

first step : relation between x1,x2
----------------------------------------------

→ first ODE 

p2a:'diff(x1,t) = x2;

p2a1:'diff(lhs(p2a),t)='diff(rhs(p2a),t);

→ second ODE (use p1b)

a)  first replace derivatives
-------------------------------------

p2b1:ev(p1b,'diff(x,t)=x2);

p2b2:ev(p2b1,'diff(x,t,2)='diff(x2,t));

b) second :variable  : x ,see listofvars 'above'
---------------------------------------------------------------

p2b3:ev(p2b2,x=x1);

c) list of 2 ODEs
-----------------------

p3:[p2a,p2b3];

→ only rhs needed: 'rkf45'

p3a:map(rhs,p3);

rem:independent varable don't change  , t 

solve : rkf45  ( need split ODE = first order ODEs)
=========

a) simulation domain (time =t)
-------------------------------------------

tstart:0$
tend  :20$

sol:rkf45(p3a,[x1,x2],[0.75,0],[t,0,20],report = true)$

b)  initial values ,see sol above
-------------------------------------------

p4:sol[1];

c)  solution sol
---------------------
datastructure → sol :[[tstart,x10,x20], ...........[.time., x1.,x2.]............]

p4[1]:/* time value */
p4[2];/* x1 value   */
p4[3];/* x2 value   */

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

 wxplot2d([[discrete,map(lambda([u],part(u,[1,2])),sol)],
               [discrete,map(lambda([u],part(u,[1,3])),sol)]],
              [style,[lines,2]],[xlabel,"t"],[ylabel,"x, x'"],
   [legend,"x(t)","x'(t)"],[title,"the van der Pol equations"]);

rem :[discrete,map(lambda([u],part(u,[1,j])),sol)]
-----------------------------------------------------------------
j=2 → [1,2]  reason :[time,x1]  =  > x(t), blue
j=3 → [1,3]  reason  [time,x1]  =  > x'(t),red

rem :[style,[lines,2]], 2 ,tickness of line. = 'option plot2d'