##               COPYRIGHT NOTICE

##  Copyright (C) 2014, 2014,  Edwin L. (Ted) Woollett
##  http://www.csulb.edu/~woollett
##  myode.R 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.R  July, 2014

##  euler1 for one first order o.d.e. --
##  one dependent variable. See example of output below code.
##  An example:   func corresponding to ode: dy/dt = -y, 
##   deriv = function(t, y) {-y}
##   grid of times at which y(t) is desired:  tL = seq(0,1,0.1)
##   soln for y(0) = 1:   yL = euler1(1,tL,deriv) 
## euler1(init,grid,func) calls func to calculate the Euler
##  solution vector. func(t,y) corresponds to independent 
##  variable = t, dependent variable = y  

euler1 = function(init, grid ,func)  {   
   n = length(grid)
   h = grid[2] - grid[1]
   w.num = vector(length = n)
   w.num[1] = init
   for (j in 1:(n-1)) {
      w.num [j+1] =  w.num [j] +
         h*func(grid[j], w.num[j]) }
   w.num}

##  Example using a coarse grid with dt = 0.1
##  > deriv = function(t, y) {-y}
##  > tL = seq(0,1,0.1)
##  > yL = euler1(1,tL,deriv)
##  > yL
##   [1] 1.0000000 0.9000000 0.8100000 0.7290000 0.6561000 0.5904900 0.5314410
##   [8] 0.4782969 0.4304672 0.3874205 0.3486784


##  myeuler(init,grid,func) for arbitrary number  of
##  first order o.d.e.'s. 

##  for one dimension: func returns a one dimen vector c(a) == a
##  tL = seq(0,3,0.1); deriv = function(t,y){ -t*y };
##  yL = myeuler(1,tL,deriv) solves dy/dt = -t*y
##  with y(0) = 1.

##  for two dimensions: func returns a two dimen vector c(a,b)
##  Example: simple harmonic oscillator: y[1] is x, y[2] is vx = dx/dt
##  solution vector for x(0) = 1, vx(0) = 0 case:
##  tL = seq(0,1,0.1)
##  sho = function(t,y) {
##     with(as.list(y), {
##        dx = y[2]
##        dvx = -4*pi^2*y[1]
##        c(dx,dvx)})}
##  out = myeuler(c(1,0), tL, sho)
##  xL = out[[1]]
##  yL = out[[2]]
##  > xL = out[[1]]
##  > xL
##   [1]  1.0000000  1.0000000  0.6052158 -0.1843525 -1.2128505 -2.1685690
##   [7] -2.6454734 -2.2662610 -0.8426575  1.4756299  4.1265851
##   myeuler returns a R list.
##    in  myeuler: each element of 'solnList' is a vector which contains the grid values
##    of one of the dependent variables.

myeuler = function(init, grid ,func)  {
   num.var = length(init)
   solnList = list()      
   n.grid = length(grid)
   for (k in 1:num.var) {
        solnList[[k]] = vector(length = n.grid)
        solnList[[k]][1] = init[k] }
   h = grid[2] - grid[1]  # step size
   yL = vector(length = num.var)  # solution at beginning of each step
   for (j in 1:(n.grid-1)) {
      for (k in 1:num.var) yL[k] = solnList[[k]][j]
      dyL = func(grid[j], yL)  #  vector of derivatives
      for (k in 1:num.var) solnList[[k]][j+1] = solnList[[k]][j] +
                                               h*dyL[k]}
   if (num.var==1) solnList[[1]] else solnList}
   
##  myrk4( init, grid, func)
##  for a single first order ode (one dependent variable), myrk4 returns a R vector containing
##    the values of that dependent variable (and func should return a single derivative value).
##  For two or more first order odes (two or more dependent variables), myrk4 returns
##    a R list of vectors, one for each dependent variable (and func should return a
##    R vector of derivative values).

myrk4 = function(init, grid ,func)  {
   num.var = length(init)
   solnList = list()      
   n.grid = length(grid)
   for (k in 1:num.var) {
        solnList[[k]] = vector(length = n.grid)
        solnList[[k]][1] = init[k] }
   h = grid[2] - grid[1]   # step size 
   yL = vector(length = num.var)  # solution at beginning of each step
   for (j in 1:(n.grid-1)) {
      for (k in 1:num.var) yL[k] = solnList[[k]][j]
      k1 = func(grid[j], yL)  #  vector of derivatives
      k2 = func(grid[j] + h/2, yL + h*k1/2)  #  vector of derivatives
      k3 = func(grid[j] + h/2, yL + h*k2/2)  #  vector of derivatives
      k4 = func(grid[j] + h, yL + h*k3)  #  vector of derivatives
      for (k in 1:num.var) solnList[[k]][j+1] = solnList[[k]][j] +
                              h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 }
   if (num.var==1) solnList[[1]] else solnList}
   
##  two dependent variables example of use of myrk4.
##  > source("myode.R")
##  > tL = seq(0,1,0.1)
##  > fll(tL)
##    0   1   11 
##  note that the func sho returns a R vector of derivative values
##  > sho = function(t,y) { c( y[2],  -4*pi^2*y[1] ) }
##  > out = myrk4(c(1,0), tL, sho)
##  > xL = out[[1]]
##  > xL
##   [1]  1.000000  0.809102  0.310104 -0.306633 -0.806047 -0.997964 -0.809517 -0.312810  0.302669
##  [10]  0.802336  0.995920
##  > vxL = out[[2]]
##  > vxL
##   [1]  0.0000000 -3.6880842 -5.9680715 -5.9724674 -3.7014458 -0.0220779  3.6627122  5.9490744
##   [9]  5.9670776  3.7117057  0.0440659
##  > tL
##   [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

##  one dependent variable example using myrk4:
##  dy/dt = -y, y(0)=1 solution vector using  function deriv and same tL grid of times. 
##  > deriv = function(t, y) {-y}
##  > out = myrk4(1,tL,deriv)
##  > is.vector(out)
##  [1] TRUE
##  > out
##   [1] 1.000000 0.904837 0.818731 0.740818 0.670320 0.606531 0.548812 0.496586 0.449329 0.406570
##  [11] 0.367880
##




##   rk4_step(init,t,dt,func) returns y(t+dt) if one dependent variable
##     or returns the vector c(x(t+dt), vx(t+dt)) if simple harmonic oscillator (two dependent
##     variables), etc. See examples of use below.

rk4_step = function(init, t, dt, func)  {
   num.var = length(init)   
   h = dt   # step size 
   yL = init   
   k1 = func(t, yL)  #  vector of derivatives
   k2 = func(t + h/2, yL + h*k1/2)  #  vector of derivatives
   k3 = func(t + h/2, yL + h*k2/2)  #  vector of derivatives
   k4 = func(t + h, yL + h*k3)  #  vector of derivatives
   for (k in 1:num.var) yL[k] = yL[k] + h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 
   if (num.var==1) yL[1] else yL}
   
##  example 1: one dependent variable, integrate dy/dt=-y, one step dt = 0.1
##  from t=0, y(0)=1, using same derivatives function deriv as used with euler1 above.
##  > rkstep = rk4_step(1, 0, 0.1, deriv)
##  > rkstep
##  [1] 0.9048375

## example2: two dependent variables, simple harmonic oscillator with unit period
##  same derivatives function sho as used above: x(0)=1, vx(0) = 0.
##  > rkstep = rk4_step(c(1,0), 0, 0.1, sho)
##  > rkstep
##  [1]  0.8091019 -3.6880842

   
##  fll is a home-made list utility: print out first, last and length of a vector
   
fll = function(xL) {
    xlen = length(xL)
    cat(" ",xL[1]," ",xL[xlen]," ",xlen,"\n") }