/*               COPYRIGHT NOTICE

Copyright (C) 2014, 2015,  Edwin L. (Ted) Woollett
http://www.csulb.edu/~woollett
myode.mac is a utility file associated with Ch. 2 of
Computational Physics with Maxima or R, Initial Value
Problems. See  cp2.pdf for more info.

This program is free software; you can redistribute
it and/or modify it under the terms of the
GNU General Public License as published by
the Free Software Foundation; either version 2 
of the License, or (at your option) any later version. 

This program is distributed in the hope that it
will be useful, but WITHOUT ANY WARRANTY;
without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. See the 
GNU General Public License for more details at
http://www.gnu.org/copyleft/gpl.html
*/


/* myode.mac  april, 2014  */

/*  euler1(dxdt,x,xi,[t,ti,tf,dt]) for one dependent variable x(t),
       xi, ti, tf, dt should evaluate to numbers.
       dxdt is an expression which can contain, potentially,
       both t and x as symbols.  */

euler1(dvar,var,init,domain) :=
block([dt,t0,n,vs,dvar0,euler_soln,r,k1,numer],numer:true,
   init : float(init),
   domain : float(domain),
   local(dvdt),
   define(dvdt(domain[1],var),float(dvar)),
   dt : domain[4],  
   t0 : domain[2],
   n: fix((domain[3] - t0)/dt),  
   vs: init,
   dvar0 : dvdt(t0, vs),
   if (not(numberp(dvar0))) then
         error("Expecting a number when the initial state is
          replaced in dvdt, but instead found:",dvar0),
   euler_soln : [[t0,vs]],
  
  for i thru n do (
      r: errcatch (k1 : dvdt(t0,vs)),
      if length(r) = 0 then return()
      else vs : vs + k1*dt,
      t0: t0 + dt,    
      euler_soln : cons([t0,vs], euler_soln)),
      
  reverse(euler_soln))$
  
  
/* yfdiff (calls euler1) returns a list of (h,yerror) pairs generated by
   the Euler solution of the o.d.e. dy/dx = dydx(x,y) with
   the initial value y(0) = 1, integrating over the
   interval [x,0,xfinal]. This is used in cp2.pdf to
   integrate the o.d.e. dy/dx = 8.5 - 20*x + 12*x^2 - 2*x^3,
   with y(0) = 1, for which ytrue = -x^4/2+4*x^3-10*x^2+17*x/2+1,
   as well as the o.d.e. dy/dx = -x*y, with y(0) = 1,
   where ytrue = exp(-x^2/2).
   See section on use of euler1()  */
   
  
yfdiff(dydx,ytrue,xfinal,hL) := 
block([tval,yerrL:[],h,esoln,yerr,numer],numer:true,
    tval : float(subst(x = xfinal,ytrue)), /* true value */
    print(" tval = ",tval),
    for h in hL do (
         esoln : euler1(dydx,y,1,[x,0,xfinal,h]),
         yerr : second(last(esoln)) - tval,
         print(" ", h, yerr),
         yerrL : cons([h, yerr], yerrL)),
    reverse(yerrL))$
    
    
/* euler_errors(n,h) uses the Euler method to generate an approximate
   numerical solution (using step size dx = h) to the first order
   o.d.e.: dy/dx = 8.5 - 20*x + 12*x^2 - 2*x^3 with the initial
   condition y(0) = 1.  At each step, the value of x, ytrue, yeuler,
   gl-err, and l-err are recorded in a table (using printf).
   gl-err (g-err in the code) is the global percent relative -
   truncation mainly -error of the Euler solution (compared with the exact true
   solution) at each step. The local truncation (mainly) error is simply
   the difference between the successive (increasing) values of
   the global error */
    
    
euler_errors(n,h) :=
block([x,ye,xn,ytrue,g_err:0, gp_err:0,l_err,numer],numer:true,
   local(yt,dydx),
   /* since x is not bound yet, these func defs work */
   define(yt(x),-x^4/2 + 4*x^3 - 10*x^2 + 17*x/2 + 1),
   define(dydx(x), 8.5 - 20*x + 12*x^2 - 2*x^3),
   x:0,
   ye : yt(x),
   printf(true,"~& ~3tx ~15tytrue ~24tyeuler ~35tgl-err ~48tl-err  ~%"),
   
   printf(true,"~& ~5f ~10t ~9f ~9f ~%",x,ye,ye),
   
   for i thru n do (
       xn : x + h,
       ye : ye + dydx(x)*h,
       ytrue : yt(xn),
       g_err : (ytrue - ye)*100/ytrue,
       l_err : g_err - gp_err,
       printf(true,"~& ~5f ~10t ~9f ~9f ~34t ~6f ~47t ~6f ~%",xn,ytrue,ye,g_err,l_err),  
       gp_err : g_err,
       x : xn))$
       
/*
(%i36) euler_errors(3,0.5)$
   x           ytrue    yeuler     gl-err       l-err  
   0.0           1.0       1.0 
   0.5       3.21875      5.25     -63.11       -63.11 
   1.0           3.0     5.875     -95.83       -32.73 
   1.5       2.21875     5.125     -131.0       -35.15 
*/
 





/*    myeuler for arbitrary number of first order o.e.d.'s

      myeuler(dxdt,x,xinit,[t,tinit,tfinal,dt])  or
      myeuler([dxdt,dydt],[x,y],[xinit,yinit],[t,tinit,tfinal,dt])  etc.
    
        syntax: same as rk, rk4
        
      The Maxima code myeuler is adapted from the code
      design of the Maxima function rk(), which can be found in
    ...share/dynamics/dynamics.mac in Maxima v. 5.28.0
    copyright 2007 Jaime E. Villate <villate@fe.up.pt>    
     */
     
     
     
myeuler(ode, var, init, domain) :=
block([uvw,duvw,esoln,n,k1,t0,dt,
            r,numer:true,display2d:false],
  init : float(init),
  domain : float(domain),         
  if (not(listp(ode))) then (
      ode : [ode],
      var : [var],
      init : [init]),  
  
  local(efunc),
  define(funmake(efunc,cons(domain[1],var)),float(ode)),
  translate(efunc),
  dt : domain[4],  
  t0 : domain[2],
  n: fix((domain[3] - t0)/dt),  
  uvw: init,
  
  duvw : apply(efunc,cons(t0,uvw)),
  if (not(numberp(last(duvw)))) then
    error("Expecting a number when the initial state
      is replaced in the equations, but instead 
      found:",last(duvw)),
      
  esoln: [cons(t0, init)],
  for i thru n do (
      r: errcatch (k1: apply(efunc,cons(t0,uvw))),      
      if length(r)=0 then return()
      else uvw: uvw + k1*dt,
      t0: t0 + dt,
      esoln : cons(cons(t0,uvw), esoln)),
  reverse(esoln))$
  






/*
    The Maxima code rk4() in this file is adapted from rk() in 
    ...share/dynamics/dynamics.mac in Maxima v. 5.28.0
    copyright 2007 Jaime E. Villate <villate@fe.up.pt>    
    $Id: dynamics.mac,v 1.6 2010-04-20 12:55:46 villate Exp $
    Villate's rk() code was converted to Lisp in v. 5.30
    
    
    
    rk4  syntax:
    if the dxdt expression is a function of (t,x), then:
       for one 1st order o.d.e: rk4(dxdt,x,xinit,[t,tinit,tfinal,dt])     
	
	If the dx1dt expression and the dx2dt expression are functions of (t,x1,x2),
    then for two 1st order o.d.e.'s:
	
    rk4([dx1dt,dx2dt],[x1,x2],[x1init,x2init],[t,tinit,tfinal,dt])     
     and so on.
	 
	 output is the list
	 [[[t0,x1(t0),x2(t0)], [t0+h,x1(t0+h),x2(t0+h)],...]
	 
     */


rk4(ode, var, init, domain) :=
block([uvw,rksoln,n,k1,k2,k3,k4,t0,t1,dt,
            r,numer:true,display2d:false],
  init : float(init),
  domain : float(domain),         
  if (not(listp(ode))) then (
      ode : [ode],
      var : [var],
      init : [init]),  
  
  local(rkfunc),
  define(funmake(rkfunc,cons(domain[1],var)),float(ode)),
  translate(rkfunc),
  dt : domain[4],  
  t0 : domain[2],
  n: fix((domain[3] - t0)/dt),  
  uvw: init,
  
  if (not(numberp(last(apply(rkfunc,cons(t0,uvw)))))) then
    error("Expecting a number when the initial state is 
    replaced in the equations, but instead found:"
       ,last(apply(rkfunc,cons(t0,uvw)))),
       
  rksoln: [cons(t0, init)],
  for i thru n do (
    r: errcatch (
      t1: domain[2]+i*dt,
      k1: apply(rkfunc,cons(t0,uvw)),
      k2: apply(rkfunc,cons((t0+t1)/2, uvw+k1*dt/2)),
      k3: apply(rkfunc,cons((t0+t1)/2,uvw+k2*dt/2)),
      k4: apply(rkfunc,cons(t1,uvw+k3*dt))),
    if length(r)=0 then return()
    else uvw: uvw + dt*(k1+2*k2+2*k3+k4)/6,
    t0: t1,
    rksoln : cons(cons(t0,uvw), rksoln)),
  reverse(rksoln))$
  
/*
(%i30) (expr: t-x^2,a:1,b:0,c:8,dt:0.1)$
(%i31) pts : rk4(expr,x,a,[t,b,c,dt])$
(%i32) fll(pts);
(%o32) [[0.0,1.0],[8.0,2.796249427928215],81]

*/
    

  
/* one dimensional examples

solves dx/dt = t - x^2 with initial condition x(0) = 1
over the time interval[0,8]

(%i1) load(myode);
(%o1) "c:/k2/myode.mac"
(%i2) pts : rk4(t-x^2,x,1,[t,0,8,0.1])$
(%i3) fll(pts);
(%o3) [[0,1],[8.0,2.796249427928215],81]
(%i5) plot2d([discrete,pts],[xlabel,"t"],[ylabel,"x"],
          [style,[lines,3]])$
          
compare soln with built-in lisp code for rk()

(%i6) pts2 : rk(t-x^2,x,1,[t,0,8,0.1])$
(%i7) diff : pts - pts2$
(%i8) fll(diff);
(%o8) [[0.0,0.0],[0.0,0.0],81]
*/
/*
test robustness: bind parameters and
  dxdt to symbols:

(%i1) load(myode);
(%o1) "c:/k2/myode.mac"
(%i2) (expr: t-x^2,a:1,b:0,c:8,dt:0.1)$
(%i3) pts : rk4(expr,x,a,[t,b,c,dt])$
(%i4) fll(pts);
(%o4) [[0,1.0],[8.0,2.796249427928215],81]

 
 use %pi as time limit:
 
(%i9) c : %pi;
(%o9) %pi
(%i10) pts : rk4(expr,x,a,[t,b,c,dt])$
(%i11) fll(pts);
(%o11) [[0,1.0],[3.1,1.666175296833841],32] 

  include %pi in definition of dt:

(%i19) dt : %pi/100;
(%o19) %pi/100
(%i20) pts : rk4(expr,x,a,[t,b,c,dt])$
(%i21) fll(pts);
(%o21) [[0.0,1.0],[3.110176727053895,1.669468096099344],100]


*/

yfdiff_rk4(dydx,ytrue,xfinal,hL) := 
block([tval,yerrL:[],h,rsoln,yerr,numer],numer:true,
    tval : float(subst(x = xfinal,ytrue)), /* true value */
    print(" tval = ",tval),
    for h in hL do (
         rsoln : rk4(dydx,y,1,[x,0,xfinal,h]),
         yerr : second(last(rsoln)) - tval,
         print(" ", h, yerr),
         yerrL : cons([h, yerr], yerrL)),
    reverse(yerrL))$




head(mylist) :=
block([numL,nleft:6],
    numL : length(mylist),
    rest (mylist, -(numL - nleft)))$
    
tail(mylist) :=
block([numL,nleft:6],
    numL : length(mylist),
    rest (mylist, numL - nleft))$



take(%aL,%nn) := (map(lambda([x],part(x,%nn)), %aL))$

range(aaL) := print(" min = ",lmin(aaL),"  max = ", lmax(aaL))$

fll(x) := [first(x),last(x),length(x)]$
  
ratprint:false$
fpprintprec:8$
display2d:false$