/*  -*-  Mode: Maxima -*- */

/*
**
** Copyright (C) 2015 Marek Rychlik <rychlik@email.arizona.edu>
**
** 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.
**
** You should have received a copy of the GNU General Public License
** along with this program; if not, write to the Free Software
** Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
**
*/

/* Variable holding diagnostic information about ODE; examine after calling LINODE */
define_variable(linode_method,[],any_check)$

/* When this switch is set, when complex characteristic values are
 present, trig functions are introduced to represent the oscillatory
 part of the solution; when the switch is off, homogenous equation is
 solved using complex exponential functions. */
define_variable(linode_prefer_real_basis,true,boolean)$


/* Solve linear ode with constant coefficients or Euler-Cauchy type.
*  The equation may be homogenous or non-homogenous.
*/
linode(de,y,x,[iconds]):=block([programmode:true,backsubst:true,m1pbranch:true,domain:complex,
  de1,de2,lambd,r,char_eqn0,char_eqn,basis,order,yp,yh,gsoln,ratio,char_nums,char_pairs,W,D,U],
  linode_method:[],
  /* if not symbolp(y) then error("Second argument",y,"is not a symbol."),*/
  if not symbolp(x) then error("Third argument",x,"is not a symbol."),
  /* Convert equation to expression */
  if inpart(de,0)="=" then de1:lhs(de)-rhs(de) else de1:de,
  /* Determine non-homogenous equation right-hand side */
  r:-ev(subst(0,y,de1),diff),
  /* Homogenous part of the equation*/
  de2:de1+r,
  /*================================================================
  * SOLVE HOMOGENOUS EQUATION 
  *================================================================*/
  /* ==== METHOD1: TRY CONSTANT COEFFICIENTS ==== */
  /* Build characteristic equation */
  char_eqn:subst(exp(lambd*x),y,de2)*exp(-lambd*x),
  char_eqn:ev(char_eqn,[diff,expand]),
  /* Check if this eliminated x */
  if freeof(x,char_eqn) then (
    push('constant_coefficients,linode_method),
    /* Find characteristic numbers */
    char_nums:solve(char_eqn,lambd),
    /* Make root/multiplicity pairs */
    char_pairs:makelist([ev(lambd,char_nums[j]),multiplicities[j]],
      j,1,length(char_nums)),
    /* Construct the basis for the homogenous equation */
    basis:[],
    while char_pairs#[] do block([p,re,im,tmp,char_num,mult],
      '[char_num,mult]::pop(char_pairs),
      re:realpart(char_num),
      im:imagpart(char_num),
      p:[conjugate(char_num),mult],
      tmp:sublist(char_pairs,lambda([q],is(equal(p-q,[0,0])))),
      if linode_prefer_real_basis and im # 0 and tmp#[] then (
        char_pairs:sublist(char_pairs,lambda([q],not member(q,tmp))),
        if is(im<0) then im:-im,
        basis:append(basis,
          apply(append,
            makelist([x^l*exp(re*x)*cos(im*x),
              x^l*exp(re*x)*sin(im*x)],l,0,mult-1)))
        )
      else (
        basis:append(basis,makelist(x^l*exp(char_num*x),l,0,mult-1))
        )
      ),
    /* Divide by the coefficient at the highest derivative */
    if r#0 then r:r/coeff(char_eqn,lambd,hipow(char_eqn,lambd)),
    go(nonhomogenous)
    ),
  /* ==== METHOD2: EULER-CAUCHY ==== */
  char_eqn0:subst(x^lambd,y,de2)*x^(-lambd),
  char_eqn0:ev(char_eqn0,[diff,expand]),
  /* Handles the case of a modified Euler-Cauchy, multiplied by a power of x */
  char_eqn:ev(char_eqn0,x=1),
  /* The ratio may be a function of x */
  ratio:factor(char_eqn0/char_eqn),
  if freeof(lambd,ratio) and freeof(x,char_eqn) then (
    push('cauchy_euler,linode_method),
    /* Find characteristic numbers */
    char_nums:solve(char_eqn,lambd),
    /* Make root/multiplicity pairs */
    char_pairs:makelist([ev(lambd,char_nums[j]),multiplicities[j]],
      j,1,length(char_nums)),
    /* Construct the basis for the homogenous equation */
    basis:[],
    while char_pairs#[] do block([p,re,im,tmp,char_num,mult],
      '[char_num,mult]::pop(char_pairs),
      re:realpart(char_num),
      im:imagpart(char_num),
      p:[conjugate(char_num),mult],
      tmp:sublist(char_pairs,lambda([q],is(equal(p-q,[0,0])))),
      if linode_prefer_real_basis and im # 0 and tmp#[] then (
        char_pairs:sublist(char_pairs,lambda([q],not member(q,tmp))),
        if is(im<0) then im:-im,
        basis:append(basis,
          apply(append,
            makelist([log(x)^l*x^re*cos(im*log(x)),
              log(x)^l*x^re*sin(im*log(x))],l,0,mult-1)))
        )
      else (
        basis:append(basis,makelist(log(x)^l*x^char_num,l,0,mult-1))
        )
      ),
    /* Right-hand side needs to be scaled before variation of parameters */
    if r#0 then r:expand(r/x^hipow(char_eqn,lambd)/ratio),
    go(nonhomogenous)
    ),
  push('unknown_type,linode_method),
  return(false),
  /*================================================================
  * SOLVE NONHOMOGENOUS EQUATION 
  *================================================================*/
  nonhomogenous,
  order:length(basis),
  /* General solution of the homogenous equation */
  yh:makelist(concat('%k,j),j,1,order).basis,
  if r=0 then (
    push('homogenous,linode_method),
    yp:0,
    go(ivp)
    ) ,
  push('nonhomogenous,linode_method),
  if order=1 then (
    push('first_order,linode_method),
    yp:integrate(r/basis[1],x)*basis[1],
    go(ivp)
    ),
  /*================================================================
  * VARIATION OF PARAMETERS 
  *================================================================*/
  push(variation_of_parameters,linode_method),
  /* Build Wronskian matrix */
  W:apply(matrix,makelist(diff(basis,x,j),j,0,order-1)),
  /* Solve W*U' = r * [0,0,...,1]^T by Cramer's Rule */
  /* Wronskian determinant. */
  D:ev(determinant(W),[expand]),
  /* For real solutions invovlving trigs, this may simplify Wronskian determinant significantly */
  if linode_prefer_real_basis then D:ev(D,[trigreduce]),
  /* Derivatives of varying parameters */
  U:expand(makelist((-1)^(order+j)*determinant(minor(W,order,j)),j,1,order)/D*r),
  /* Varying parameters */
  U:map(lambda([u],integrate(u,x)),U),
  /* Particular solution */
  yp:expand(U.basis),
  /* Attempt simplification of trigs, if real basis is used */
  if linode_prefer_real_basis then yp:ev(yp,[trigreduce]),
  /*================================================================
  * SOLVE INITIAL VALUE PROBLEM
  *================================================================*/
  ivp,
  /* Form general solution */
  gsoln:y=yp+yh,
  /* Solve initial conditions, if given */
  if iconds#[] then (
    apply(bc,append([gsoln],iconds))
    )
  else  gsoln
  )$


noteqn(x):=if atom(x) or not inpart(x,0)="="
	     then (disp(x), disp("not an equation"), error())$

/* Check if a symbolic variable is a parameter produced by LINODE, i.e. %k1, %k2, ... */
paramp(var):=block([x],
  if not symbolp(var) then return(false),
  x:string(var),
  slength(x) >= 2 and substring(x,1,3) = "%k")$

/* Solves an initial/boundary/multiple point value problem starting
from the general solution obtained with LINODE.  It tries to simplify
the answer, and produce a real solution, if possible.  */
bc(gsoln,[bdryconds]):= block([programmode:true,backsubst:true,singsolve:true,
  eqns:[],xa,yavals,eqnsa,sub,vars],
  if oddp(length(bdryconds)) then error("Odd number of initial/boundary condition arguments."),
  while bdryconds # [] do (
    xa:first(bdryconds),
    yavals:second(bdryconds),
    bdryconds:rest(rest(bdryconds)),
    noteqn(xa),
    every(noteqn, yavals), 
    eqnsa:subst(gsoln,yavals),
    eqnsa:ev(eqnsa,[diff,expand]),
    eqnsa:at(eqnsa,xa),
    /* NOTE: This is a workaround due to a bug in pquotient in file rat3a.lisp.
    It only seems to affect the form involving trig functions */
    eqns:append(eqns, eqnsa)
    ),
  /* Get unbound parameters %k1, %k2, ... */
  vars:sublist(listofvars(gsoln),paramp),
  /* Construct solution to the linear system */
  sub:linsolve(eqns,vars),
  /* Substitute solutions for constants into supplied solution */
  if sub=[] then [] else expand(ratsimp(demoivre(subst(sub,gsoln))))
  )$