module odenonn$ % Special form nonlinear ODEs of order > 1
% F.J.Wright@maths.qmw.ac.uk, Time-stamp: <14 August 2001>
% Trivial order reduction.
% Special cases of Lie symmetry, namely
% autonomous, equidimensional and scale invariant equations.
% Simplification of arbitrary constants.
% TO DO:
% avoid computing orders in both reduce and shift
algebraic procedure ODESolve!-nonlinearn(ode, y, x);
%% `symbolic' mode here is ESSENTIAL, otherwise the code generated
%% is as if this were a macro, except then it does not get eval'ed!
symbolic ODENon!-Reduce!-Order(ode, y, x)$
%% The following defines are used to allow easy changes to the calling
%% sequence.
define ODENon!-Reduce!-Order!-Next = ODESolve!-Shift$
%% Shifting currently NEEDED for Zimmer (8) (only)!
define ODESolve!-Shift!-Next = ODESolve!-nonlinearn!*1$
switch odesolve_equidim_y$ % TEMPORARY?
symbolic(!*odesolve_equidim_y := t)$ % TEMPORARY?
algebraic procedure ODESolve!-nonlinearn!*1(ode, y, x);
%% The order here seems to be important in practice:
symbolic or(
ODESolve!-Autonomous(ode, y, x),
ODESolve!-ScaleInv(ode, y, x), % includes equidim in x
!*odesolve_equidim_y and
ODESolve!-Equidim!-y(ode, y, x) )$
algebraic procedure ODENon!-Reduce!-Order(ode, y, x);
%% If ode does not explicitly involve y and low order derivatives
%% then simplify by reducing the effective order (unless there is
%% only one) and then try to solve the reduced ode directly to give
%% a first integral. Applies only to odes of order > 1.
begin scalar deriv_orders, min_order, max_order;
traceode1 "Trying trivial order reduction ...";
deriv_orders := get_deriv_orders(ode, y);
%% Check for purely algebraic factor from some simplification,
%% such as autonomous reduction:
if deriv_orders = {} or deriv_orders = {0} then
return {ode = 0}; % purely algebraic!
%% Avoid reduction to a purely algebraic equation:
if (min_order := min deriv_orders) = 0 or
length deriv_orders = 1 then return
ODENon!-Reduce!-Order!-Next(ode, y, x);
max_order := max deriv_orders;
ode := sub(df = odesolve!-df, ode);
for ord := min_order : max_order do
ode := if ord = 1 then (ode where odesolve!-df(y,x) => y)
else (ode where odesolve!-df(y,x,ord) =>
odesolve!-df(y,x,ord-min_order));
ode := sub(odesolve!-df = df, ode);
traceode "Performing trivial order reduction to give ",
"the order ", max_order - min_order, " nonlinear ODE: ",
ode = 0;
ode := symbolic(
(if max_order - min_order = 1 then % first order
ODESolve!-nonlinear1(ode, y, x)
else
ODENon!-Reduce!-Order!-Next(ode, y, x))
where !*odesolve_explicit = t);
if not ode then <<
traceode "Cannot solve order-reduced ODE!";
return % abandon solution
>>;
%% ode := sub(y = df(y,x,min_order), ode);
traceode "Solution of order-reduced ODE is ", ode;
traceode "Restoring order, ", y => df(y,x,min_order),
", to give: ", sub(y = df(y,x,min_order), ode),
" and re-solving ...";
ode := for each soln in ode join
%% Each `soln' here is an EQUATION for y that may be
%% implicit.
if lhs soln = y then % explicit
{ y = ODESolve!-multi!-int!*(rhs soln, x, min_order) }
else <<
soln := solve(soln, y);
for each s in soln collect
if lhs s = y then % explicit
y = ODESolve!-multi!-int!*(rhs s, x, min_order)
else % implicit
%% leave unsolved for now
sub(y = df(y,x,min_order), s)
>>;
return ODESolve!-Simp!-ArbConsts(ode, y, x)
end$
algebraic procedure ODESolve!-multi!-int!*(y, x, m);
%% Integate y wrt x m times and add arbitrary constants:
ODESolve!-multi!-int(y, x, m) +
%% << >> below is NECESSARY to force immediate evaluation!
for i := 0 : m-1 sum <<newarbconst()>>*x^i$
% Internal wrapper function for ODESolve!-Shift:
algebraic operator odesolve!-sub!*$
algebraic procedure ODESolve!-Shift(ode, y, x);
%% A first attempt at canonicalizing an ODE by shifting the
%% independent variable.
symbolic if not !*odesolve_fast then % heuristic solution
algebraic begin scalar deriv_orders, a, c, d;
traceode1 "Looking for an independent variable shift ...";
deriv_orders := get_deriv_orders(ode, y);
deriv_orders := sort(deriv_orders, >);
%% Try to find a non-trivial "coefficient" polynomial
%% constituent c that is linear in x.
while deriv_orders neq {} and
(c := lcof(ode, df(y,x,first deriv_orders))) freeof x do
deriv_orders := rest deriv_orders;
if deriv_orders = {} then % not shiftable
return ODESolve!-Shift!-Next(ode, y, x);
if (d := deg(c, x)) neq 1 then <<
c := decompose c;
while (c := rest c) neq {} and deg(rhs first c, x) neq 1 do;
%% << null loop body >>
if c neq {} then c := rhs first c;
if deg(c, x) neq 1 then % not shiftable
return ODESolve!-Shift!-Next(ode, y, x)
>>;
%% c = ax + b is a linear component polynomial of the ode
%% coefficients.
if not(c freeof y) or not((c := coeff(c,x)) freeof x) or
first c = 0 then % not shiftable
return ODESolve!-Shift!-Next(ode, y, x);
c := first c / (a := second c);
%% ode is a function of ax + b (= a(x + c)), so ...
ode := sub(x = x - c, ode) / a^d;
%% This will leave sub(..., df(...)) symbolic, so ...
ode := num sub(sub = odesolve!-sub!*, ode);
ode := (ode where odesolve!-sub!*(~a! ,~b! ) => b! );
traceode "This ODE can be simplified by the ",
"independent variable shift ", x => x - c,
" to give: ", ode = 0;
ode := ODESolve!-Shift!-Next(ode, y, x);
if ode then return
for each soln in ode collect
if symbolic rlistp soln then % parametric solution
for each s in soln collect
if symbolic eqcar(s, 'equal) and lhs s = x
then x = rhs s - c else s
else sub(x = x + c, soln)
end$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Autonomous, equidimensional and scale-invariant ODEs
% ====================================================
algebraic procedure ODESolve!-Autonomous(ode, y, x);
%% If ODE is autonomous, i.e. x does not appear explicitly, then
%% reduce the order by using y as independent variable and then try
%% to solve the reduced ODE directly. Applies only to ODEs of
%% order > 1. Do not apply to a linear ODE, because it will become
%% nonlinear!
begin scalar ode1, u, soln;
traceode1 "Testing whether ODE is autonomous ...";
ode1 := (ode where df(y,x) => 1, df(y,x,~n) => 1);
if smember(x, ode1) then return; % not autonomous
u := gensym();
symbolic depend1(u,x,t); symbolic depend1(u,y,t);
ode := (ode where df(y,x) => u, df(y,x,~n) => df(u,x,n-1));
ode := (ode where df(u,x,~n) => u*df(df(u,x,n-1),y) when n > 1,
%% above condition n > 1 is NECESSARY!
df(u,x) => u*df(u,y));
symbolic depend1(u,x,nil);
traceode
"This ODE is autonomous -- transforming dependent variable ",
"to derivative to give this ODE of order 1 lower: ", ode = 0;
ode := symbolic(ODESolve!*1(ode, u, y)
where !*odesolve_explicit = t);
if not ode then <<
symbolic depend1(u,y,nil);
traceode "Cannot solve transformed autonomous ODE!";
return
>>;
ode := sub(u = df(y,x), ode);
symbolic depend1(u,y,nil);
traceode "Restoring order to give these first-order ODEs ...";
soln := {};
a: if ode neq {} then
if (u := ODESolve!-FirstOrder(first ode, y, x)) then <<
soln := append(soln, u);
ode := rest ode;
go to a
>> else <<
traceode "Cannot solve one of the first-order ODEs ",
"arising from solution of transformed autonomous ODE!";
return
>>;
return ODESolve!-Simp!-ArbConsts(soln, y, x)
end$
algebraic procedure ODESolve!-ScaleInv(ode, y, x);
%% If ODE is scale invariant, i.e. invariant under x -> a x, y ->
%% a^p y, then transform it to an equidimensional-in-x ODE and try
%% to solve it. If p = 0 then it is already equidimensional-in-x
%% as a special case. Returns a solution or nil if this method
%% does not lead to a solution. PROBABLY NOT USEFUL FOR LINEAR
%% ODES.
begin scalar u, p, ode1, pow, !*allfac;
traceode1 "Testing whether ODE is scale invariant or ",
"equidimensional in the independent variable ", x, " ...";
u := gensym(); p := gensym();
ode1 := (ode where df(y,x,~n) => mkid(u,n)*x^(p-n),
df(y,x) => mkid(u,1)*x^(p-1));
%% mkid's to avoid spurious cancellations.
ode1 := num sub(y = u*x^p, ode1);
%% Try to choose p to make ode1 proportional to some single
%% power of x. Assume ode1 is a sum of terms.
begin scalar part1, n_parts;
part1 := part(ode1, 1);
n_parts := arglength ode1; % must be at least 2 terms
for i := 2 : n_parts do <<
parti := part(ode1, i)/part1;
pow := df(parti, x)*x/parti;
if pow then << pow := solve(pow, p); n_parts := 0 >>
>>;
if n_parts then
%% Scale invariant for ANY p =>
%% equidimensional in both x and y
return pow := {p=0};
ode1 := (ode1 - part1)/part1;
while pow neq {} and % check scale invariance
(symbolic eqcar(caddr cadr pow, 'root_of) or
not(sub(first pow, ode1) freeof x)) do
pow := rest pow
end;
if pow = {} then return; % not scale invariant
if not(p := rhs first pow) then
%% Scale invariant for p=0 =>
%% equidimensional in x ...
return ODESolve!-ScaleInv!-Equidim!-x(ode, y, x);
%% ode is scale invariant (with p neq 0)
symbolic depend1(u, x, t);
ode := sub(y = x^p*u, ode);
traceode "This ODE is scale invariant -- applying ", y => x^p*u,
" to transform to the simpler ODE: ", ode = 0;
ode := ODESolve!-ScaleInv!-Equidim!-x(ode, u, x);
symbolic depend1(u, x, nil);
if ode then return sub(u = y/x^p, ode);
traceode "Cannot solve transformed scale invariant ODE!"
end$
algebraic procedure ODESolve!-ScaleInv!-Equidim!-x(ode, y, x);
%% ODE is equidimensional in x, i.e. invariant under x -> ax, so
%% transform it to an autonomous ODE and try to solve it. (This
%% includes "reduced" Euler equations as a special case. Could
%% ignore terms independent of y in testing equidimensionality; if
%% there are such terms then the simplified ODE will not be
%% autonomous. This generalization includes ALL Euler equations.)
%% Returns a solution or nil if this method does not lead to a
%% solution.
begin scalar tt, exp_tt;
tt := gensym();
%% ode is equidimensional in x; x -> exp(tt):
exp_tt := exp(tt);
symbolic depend1(y,tt,t);
ode := (ode where df(y,x) => df(y,tt)/exp_tt,
df(y,x,~n) => df(df(y,x,n-1),tt)/exp_tt when
numberp n and n > 0); % n > 0 condition is necessary!
ode := num sub(x = exp_tt, ode);
traceode
"This ODE is equidimensional in the independent variable ",
x, " -- applying ", x => exp_tt,
" to transform to the simpler ODE: ",
ode = 0;
symbolic depend1(y,x,nil); % Necessary to avoid dependence loops
%% ode should be autonomous PROVIDED no term independent of y
ode := symbolic ODESolve!-Autonomous(ode, y, tt);
symbolic depend1(y,x,t); %%% ???
symbolic depend1(y,tt,nil);
if ode then return sub(tt = log x, ode);
traceode "Cannot solve transformed equidimensional ODE!"
end$
algebraic procedure ODESolve!-Equidim!-y(ode, y, x);
%% If ODE is equidimensional in y, i.e. invariant under y -> ay,
%% then simplify the ODE and try to solve the result. Returns a
%% solution or nil if this method does not lead to a solution. Do
%% not apply to a linear ODE, which is trivially equidimensional in
%% y, because it will become nonlinear!
begin scalar ode1, u, exp_u;
traceode1 "Testing whether ODE is equidimensional in ",
"the dependent variable ", y, " ...";
u := gensym(); % to avoid spurious cancellations
ode1 := (ode where df(y,x,~n) => y*mkid(u,n),
df(y,x) => y*u);
%% ode1 must be proportional to some single positive integer
%% power of y:
if reduct(ode1, y) or depends(lcof(ode1, y), y) then return;
%% ode is equidimensional in y; y -> exp(u):
exp_u := exp(u);
symbolic depend1(u,x,t);
ode := lcof(num sub(y = exp_u, ode), exp_u);
%% (Lcof above to remove irrelevant factor of a power of y.)
traceode
"This ODE is equidimensional in the dependent variable ",
y, " -- applying ", y => exp_u,
" to transform to the simpler ODE: ",
ode = 0;
%% ode here could be ANY kind of ode. It should be less
%% nonlinear, but I don't think there is any guarantee that it
%% will be linear -- is there? Hence we must call the full
%% general ode solver again:
ode := ODESolve!*1(ode, u, x);
symbolic depend1(u,x,nil);
if not ode then <<
traceode "Cannot solve transformed equidimensional ODE!";
return
>>;
return for each soln in ode collect
if lhs soln = u then y = exp rhs soln % retain explicit soln
else sub(u = log y, ode)
end$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Simplification of Arbitrary Constants
% =====================================
algebraic procedure ODESolve!-Simp!-ArbConsts(solns, y, x);
%% To be applied to non-parametric solutions of ODEs containing
%% arbconsts introduced earlier.
for each soln in solns collect
if symbolic rlistp soln then soln else
(ODESolve!-Simp!-ArbConsts1(lhs soln, y, x) =
ODESolve!-Simp!-ArbConsts1(rhs soln, y, x))$
algebraic procedure ODESolve!-Simp!-ArbConsts1(soln, y, x);
%% Simplify arbconst expressions within soln. Messy arbconst
%% expressions can be introduced by the integrator from simple
%% arbconsts, and there would appear to be no way to avoid this
%% other than to remove them after the event.
begin scalar !*precise, ss, acexprns;
if not(ss := ODESolve!-Structr(soln, x, y, 'arbconst))
then return soln;
acexprns := rest ss; ss := first ss;
traceode "Simplifying the arbconst expressions in ", soln,
" by the rewrites ...";
for each s in acexprns do <<
%% s has the form ansj = "expression in arbconst(n)"
%% MAY NEED TO CHECK ONLY 1 ARBCONST?
%% n!* must be a global algebraic variable!
rhs s where arbconst(~n) => (n!* := n);
traceode rhs s, " => ", arbconst(n!*); % to evaluate `rhs'
%% Remove other occurrences of arbconst(n!*):
ss := sub(solve(s, arbconst(n!*)), ss);
%% Finally rename ansj as arbconst(n):
ss := sub(lhs s = arbconst(n!*), ss)
>>;
return ss
end$
symbolic operator ODESolve!-Structr$
%% symbolic procedure ODESolve!-Structr(u, x, y, arbop);
%% %% Return an rlist consisting of an expression involving variables
%% %% ansj representing sub-structures followed by equations of the
%% %% form ansj = sub-structure, where the sub-structures depend
%% %% non-trivially on the arbitrary opertor arbop, essentially in the
%% %% format returned by structr, or nil if this decomposition is not
%% %% possible.
%% begin scalar !*savestructr, !*precise, ss, arbexprns;
%% !*savestructr := t;
%% ss := cdr structr u; % rlistat; ss = (exprn eqns)
%% if null cdr ss then return;
%% %% Ignore "structures" of the form ansj = arbop(i)
%% ss := car ss . for each s in cdr ss join
%% if eqcar(caddr(s:=reval s), arbop)
%% then << arbexprns := s . arbexprns; nil >>
%% else {s};
%% %% by substituting them back into the structure list:
%% if arbexprns then
%% ss := cdr subeval nconc(arbexprns, {makelist ss});
%%
%% %% Get simplifiable arbop expressions:
%% arbexprns := nil;
%% ss := car ss . for each s in cdr ss join
%% begin scalar rhs_s; rhs_s := caddr s;
%% return if smember(arbop, rhs_s) and
%% not(depends(rhs_s, x) or depends(rhs_s, y))
%% then << arbexprns := s . arbexprns; nil >>
%% else {s}
%% end;
%% if null arbexprns then return;
%% %% Rebuild the rest of the stucture as ss:
%% ss := if cdr ss then
%% subeval nconc(cdr ss, {car ss})
%% else car ss;
%% return makelist(ss . arbexprns)
%% end$
symbolic procedure ODESolve!-Structr(u, x, y, arbop);
%% Return an rlist representing U that consists of an expression
%% involving variables `ansj' representing sub-structures followed
%% by equations of the form `ansj = sub-structure', where the
%% sub-structures depend non-trivially on the arbitrary opertor
%% ARBOP and do not depend on X or Y, essentially in the format
%% returned by structr, or nil if this decomposition is not
%% possible.
begin scalar !*savestructr, !*precise, ss, arbexprns;
!*savestructr := t;
ss := cdr structr u; % rlistat; ss = (exprn eqns)
if null cdr ss then return;
%% Ignore trivial structure of the form ansj = arbop(i) by
%% substituting it back into the structure list:
ss := car ss . for each s in cdr ss join
if eqcar(caddr(s:=reval s), arbop)
then << arbexprns := s . arbexprns; nil >>
else {s};
if null cdr ss then return;
if arbexprns then
ss := cdr subeval nconc(arbexprns, {makelist ss});
%% Ignore structure that does not depend on arbop by
%% substituting it back into the structure list:
arbexprns := nil;
ss := car ss . for each s in cdr ss join
if not smember(arbop, s)
then << arbexprns := s . arbexprns; nil >>
else {s};
if null cdr ss then return;
if arbexprns then
ss := cdr subeval nconc(arbexprns, {makelist ss});
%% Ignore all structure that depends on x or y by repeatedly
%% substituting it back into the structure list:
arbexprns := t;
while arbexprns and cdr ss do <<
arbexprns := nil;
ss := car ss . for each s in cdr ss join
if depends(s, x) or depends(s, y)
then << arbexprns := s . arbexprns; nil >>
else {s};
if arbexprns then
ss := cdr subeval nconc(arbexprns, {makelist ss})
>>;
if null cdr ss then return;
return makelist ss
end$
endmodule$
end$