Artifact f9d6912199006f7d689baf3b29b1a313dbce742783c2a59ffef49aa25c6e2815:
- Executable file
r38/packages/odesolve/odenon1.red
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 33834) [annotate] [blame] [check-ins using] [more...]
module odenon1$ % Special form nonlinear ODEs of order 1 % F.J.Wright@maths.qmw.ac.uk, Time-stamp: <11 August 2001> % Original author: Malcolm MacCallum % Basic layout is to test first for well-known special forms, namely: % separable % quasi-separable (separable after linear transformation) % (algebraically) homogeneous % quasi-homogeneous (homogeneous after linear transformation) % Bernoulli % Riccati % solvable for x or y, including Clairaut and Lagrange % exact (FJW: general exact ode now handled elsewhere) % and at a later stage of development to add more general methods, % such as Prelle Singer % Note that all cases of first order ODEs can be considered equivalent % to searches for integrating factors. (In particular Lie methods do % not help here.) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% algebraic procedure odesolve(ode, y, x); %% %% Temporary definition for test purposes. %% begin scalar !*precise, !*odesolve!-solvable!-xy, solution; %% ode := num !*eqn2a ode; % returns ode as expression %% if (solution := ODESolve!-NonLinear1(ode, y, x)) then %% return solution %% else %% write "***** ODESolve cannot solve this ODE!" %% end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global '(ODESolve_Before_Non1Grad_Hook ODESolve_After_Non1Grad_Hook)$ algebraic procedure ODESolve!-NonLinear1(ode, y, x); %% Top-level solver for non-linear first-order ODEs. begin scalar odecoeffs, gradient, solution, p, ode_p; traceode1 "Entering top-level non-linear first-order solver ..."; %% traceode "This is a nonlinear ODE of order 1."; p := gensym(); ode_p := num sub(df(y,x) = p, ode); %% If ode contains exponentials then the above sub can produce a %% denominator when there is none in ode! odecoeffs := coeff(ode_p, p); if length odecoeffs = 2 and not smember(p, odecoeffs) then << % first DEGREE ODE gradient := -first odecoeffs / second odecoeffs; symbolic if (solution := or( ODESolve!-Run!-Hook( 'ODESolve_Before_Non1Grad_Hook, {gradient, y, x}), ODESolve!-Separable(gradient, y, x), ODESolve!-QuasiSeparable(gradient, y, x), ODESolve!-Homogeneous(gradient, y, x), ODESolve!-QuasiHomog(gradient, y, x), ODESolve!-Bernoulli(gradient, y, x), ODESolve!-Riccati(gradient, y, x), ODESolve!-Run!-Hook( 'ODESolve_After_Non1Grad_Hook, {gradient, y, x}))) then return solution >>; %% If ode degree neq 1 or above solvers fail then ... %% A first-order ode is "solvable-for-y" if it can be put into %% the form y = f(x,p) where p = y'. This form includes %% Clairaut and Lagrange equations as special cases. The %% Lagrange form is y = xF(y') + G(y'). If F(p) = p then this %% is a Clairaut equation. It is "solvable-for-x" if it can be %% put into the form x = f(y,p). if (solution := ODESolve!-Clairaut(ode, ode_p, p, y, x)) then return solution; %% Avoid infinite loops: symbolic if !*odesolve!-solvable!-xy then return; symbolic(!*odesolve!-solvable!-xy := t); %% "Solvable for y" includes Lagrange as a special case: symbolic return ODESolve!-Solvable!-y(ode_p, p, y, x) or ODESolve!-Solvable!-x(ode_p, p, y, x) end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Support routines algebraic procedure ODENon!-Linear1(ode, y, x); %% Solve the linear ODE: dy/dx = gradient(x,y) = P(x)*y + Q(x) %% Split into 2 procedures by FJW. begin scalar gradient; gradient := coeff(num ode, df(y,x)); % {low, high} gradient := -first gradient/second gradient; traceode!* "This is a first-order linear ODE solved by "; return if smember(y, gradient) then begin scalar P, Q; traceode "the integrating factor method."; P := lcof(num gradient,y)/den gradient; Q := gradient - P*y; return { y = ODENon!-Linear1PQ(P, Q, x) } end else << traceode "quadrature."; % FJW: Optionally turn off final integration: { y = ODESolve!-Int(gradient, x) + newarbconst() } >> end$ algebraic procedure ODENon!-Linear1PQ(P, Q, x); %% Solve the linear ODE: dy/dx = P(x)*y + Q(x) %% Called directly by ODESolve!-Bernoulli begin scalar intfactor, !*combinelogs; %% intfactor simplifies better if logs in the integral are %% combined: symbolic(!*combinelogs := t); intfactor := exp(int(-P, x)); %% Optionally turn off final integration: return (newarbconst() + ODESolve!-Int(intfactor*Q,x))/intfactor end$ %% algebraic procedure unfactorize factorlist; %% %% Multiply out a factor list of the form %% %% { {factor, multiplicity}, ... }: %% for each fac in factorlist product first fac^second fac$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Separable ODEs algebraic procedure ODESolve!-Separable(gradient, y, x); %% The ODE has the form dy/dx = gradient = F(x,y). If F(x,y) = %% f(x)g(y) then the ODE is separable, in which case return the %% solution; otherwise return nil. begin scalar f, g, !*redefmsg; traceode1 "Testing for a separable ODE ..."; %% Handle implicit dependence on x (ignoring that via y): symbolic (<< depend1(y, x, nil); %% Hack sub to handle implicit dependences: copyd('ODESolve!-old!-subsublis, 'subsublis); copyd('subsublis, 'ODESolve!-subsublis); g := errorset!*( {'ODESolve!-Separable1, mkquote gradient, mkquote x}, nil); copyd('subsublis, 'ODESolve!-old!-subsublis); if errorp g then RedErr {"(in ODESolve!-Separable1)", emsg!*}; g := car g; if depends(g, x) then g := nil; >> where depl!* = depl!*); if not g then return; %% Then F(alpha,y) = f(alpha)g(y), and F(x,y)/F(alpha,y) = %% f(x)/f(alpha) is free of y: if depends(f := gradient/g, y) then return; traceode "It is separable."; %% Any separation of F(x,y) will suffice! %% gradient := int(1/g, y) - int(f, x) + newarbconst(); %% Try to optimize structure of solution to avoid a likely %% logarithm of a function of y: gradient := int(1/g, y); gradient := if part(gradient,0) = log then part(gradient,1) - newarbconst()*exp(int(f, x)) else gradient - int(f, x) + newarbconst(); return { num gradient = 0 } end$ algebraic procedure ODESolve!-Separable1(gradient, x); %% Find a small constant alpha such that F(alpha,y) exists and %% F(alpha,y) neq 0, and return F(alpha,y), where F = gradient: begin scalar numer, denom, alpha, d, n; numer := num gradient; denom := den gradient; alpha := 0; while (d:=sub(x=alpha,denom)) = 0 or (n:=sub(x=alpha,numer)) = 0 do alpha := alpha + 1; return n/d end$ symbolic procedure ODESolve!-subsublis(u,v); % NOTE: This definition assumes that with the exception of *SQ and % domain elements, expressions do not contain dotted pairs. % This is the standard `subsublis' plus support for one level of % indirect dependence (cf. freeof), in which case sub converts the % dependent variable into an independent variable with a different % but related name. % u is an alist of substitutions, v is an expression to be % substituted: begin scalar x; return if x := assoc(v,u) then cdr x % allow for case of sub(sqrt 2=s2,atan sqrt 2). else if eqcar(v,'sqrt) and (x := assoc(list('expt,cadr v,'(quotient 1 2)),u)) then cdr x else if atom v then if (x := assoc(v,depl!*)) and % FJW %% Run through variables on which v depends: << while (x := cdr x) and not assoc(car x,u) do; x >> % FJW then mkid(v, '!!) else % FJW v else if not idp car v then for each j in v collect ODESolve!-subsublis(u,j) else if x := get(car v,'subfunc) then apply2(x,u,v) else if get(car v,'dname) then v else if car v eq '!*sq then ODESolve!-subsublis(u,prepsq cadr v) else for each j in v collect ODESolve!-subsublis(u,j) end$ algebraic procedure ODESolve!-QuasiSeparable(gradient, y, x); %% The ODE has the form dy/dx = gradient = F(x,y). If F(x,y) = %% f(y+kx) then the ODE is quasi-separable, in which case return %% the solution; otherwise return nil. begin scalar k; traceode1 "Testing for a quasi-separable ODE ..."; %% F(x,y) = f(y+kx) iff df(F,x)/df(F,y) = k = constant (using %% PARTIAL derivatives): k := (df(gradient,x)/df(gradient,y) where df(y,x) => 0); if depends(k, x) then return; % sufficient since y = y(x) traceode "It is separable after letting ", y+k*x => y; %% Setting u = y+kx gives du/dx = dy/dx+k = f(u)+k: gradient := sub(x=0, gradient) + k; %% => int(1/(f(u)+k),u) = x + arbconst. gradient := sub(y=y+k*x, int(1/gradient, y)) - x + newarbconst(); return { num gradient = 0 } end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Algebraically homogeneous ODEs algebraic procedure ODESolve!-Homogeneous(gradient, y, x); %% The ODE has the form dy/dx = gradient = F(x,y). If F(x,y) = %% f(y/x) then the ODE is algebraically homogeneous. %% Setting y = vx => v + x dv/dx = F(x,vx) = f(v) begin scalar v; v := gensym(); traceode1 "Testing for a homogeneous ODE ..."; gradient := sub(y=v*x, gradient); if depends(gradient, x) then return; traceode "It is of algebraically homogeneous type ", "solved by a change of variables of the form `y = vx'."; %% Integrate to give int(1/(f(v)-v),v) = int(1/x, x) + arbconst %% Hence exp int(1/(f(v)-v),v) = arbconst * x gradient := exp(int(1/(gradient-v), v)); gradient := (gradient where (sqrt ~a)*(sqrt ~b) => sqrt(a*b)); gradient := sub(v=y/x, x/gradient) + newarbconst(); return {num gradient = 0} end$ %% A quasi-homogeneous first-order nonlinear ODE has the form %% dy/dx = F(..., ((a1*x + b1*y + c1)/(a2*x + b2*y + c2))^p, ...) %% where F is an arbitrary composition of functions and p is an %% arbitrary power. F may have other arguments that do not depend on %% x or y or that depend in the same way (e.g. F may be a sum or %% product). The argument of F that depends on x or y will be a %% quotient of expanded polynomials if p is a positive integer. %% Otherwise, it will be a power of a quotient (expt (quotient ... )), %% in which case the `expt' can be treated as part of the composition %% F, or a quotient of powers (quotient (expt ... ) (expt ... )) which %% must be treated separately. algebraic procedure ODESolve!-QuasiHomog(gradient, y, x); %% The ODE has the form dy/dx = gradient = F(x,y). If F(x,y) = %% f((a1*x + b1*y + c1)/(a2*x + b2*y + c2)) where the function f %% may be arbitrary then the ODE is reducible to algebraically %% homogeneous. Find the first linear-rational sub-expression, and %% use it to try to make the ODE homogeneous: begin scalar tmp, n, d, soln; traceode1 "Testing for a quasi-homogeneous ODE ..."; %% First, find an "argument" that is a rational function, with %% numerator and denominator that both depend on x. if not(tmp := symbolic ODESolve!-QuasiHomog1(reval gradient, x)) then return; n := num tmp; d := den tmp; %% Now check that numerator and denominator have the same degree %% in x (and y): if (tmp := deg(n,x)) neq deg(d,x) then return; %% If that degree > 1 then extract the squarefree factor: if tmp = 1 then % => deg(d,x) = 1 %% ( if deg(n,y) neq 1 or deg(d,y) neq 1 then return ) else << %% Use partial squarefree factorization to find p'th root of %% numerator and denominator: %% If f = poly = (a x + b y + c)^p %% then f' = a p(a x + b y + c)^(p-1) %% => gcd(f, f') = (a x + b y + c)^(p-1) %% => f / gcd(f, f') = a x + b y + c. n := n / gcd(n, df(n, x)); % must be linear in x and y %% if deg(n,x) neq 1 or deg(n,y) neq 1 then return; d := d / gcd(d, df(d, x)); % must be linear in x and y %% if deg(d,x) neq 1 or deg(d,y) neq 1 then return >>; %% Check that numerator and denominator are really LINEAR %% POLYNOMIALS in x and y (where y depends on x): if length(tmp := coeff(n, y)) neq 2 or depends(tmp, y) then return; if depends(second tmp, x) or % coeff of y^1 length(tmp := coeff(first tmp, x)) neq 2 or % coeff of y^0 depends(tmp, x) then return; if length(tmp := coeff(d, y)) neq 2 or depends(tmp, y) then return; if depends(second tmp, x) or % coeff of y^1 length(tmp := coeff(first tmp, x)) neq 2 or % coeff of y^0 depends(tmp, x) then return; %% The degenerate case a1*b2=a2*b1 is now treated as quasi-separable. traceode "It is quasi-homogeneous if ", "the result of shifting the origin is homogeneous ..."; soln := first solve({n,d},{x,y}); n := rhs first soln; d := rhs second soln; gradient := sub(x=x+n, y=y+d, gradient); %% ODE was quasi-homogeneous iff the new ODE is homogeneous: if (soln := ODESolve!-Homogeneous(gradient,y,x)) then return sub(x=x-n, y=y-d, soln); traceode "... which it is not!" end$ %%% The calls to `depends' below are inefficient! symbolic procedure ODESolve!-QuasiHomog1(u, x); %% Assumes "algebraic" form! Get the first argument of any %% composition of functions that is a quotient of polynomials or %% symbolic powers (expt forms) that both depend on `x'. if atom u then nil % not QH, else operator ... else if car u eq 'quotient and % quotient, in which depends(cadr u, x) and % numerator depends on x, and depends(caddr u, x) then % denominator depends on x. if eqcar(cadr u, 'expt) then % symbolic powers ( if eqcar(caddr u, 'expt) and (caddr cadr u eq caddr caddr u) then {'quotient, cadr cadr u, cadr caddr u} ) else ( if eqcar(cadr u, 'plus) and eqcar(caddr u, 'plus) then u ) % polynomials (?) else % Process first x-dependent argument of operator u: begin a: if (u := cdr u) then if depends(car u, x) then return ODESolve!-QuasiHomog1(car u, x) else go to a end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Bernoulli ODEs symbolic operator ODESolve!-Bernoulli$ symbolic procedure ODESolve!-Bernoulli(rhs, y, x); %% The ODE has the form df(y,x) = rhs. If rhs has the Bernoulli %% form P(x)*y + Q(x)*y^n then extract P(x), Q(x), n and return the %% solution else return nil. ( begin scalar num_rhs, den_rhs, C1, C2, d, d1, d2, d3, P, Q, n; traceode1 "Testing for a Bernoulli ODE ..."; %% Degrees will be constructed in true prefix form. %% Need sum of two terms, both with main var (essentially) y. depend1(x, y, nil) where !*msg = nil; << updkorder y; rhs := simp rhs >> where kord!* = kord!*; num_rhs := numr rhs; den_rhs := denr rhs; if domainp num_rhs or not(mvar num_rhs eq y) then return; %% Num may have the form y^d * (y^d1 C1(x) + y^d2 C2(x)) ... if null red num_rhs then << d := ldeg num_rhs; num_rhs := lc num_rhs; if domainp num_rhs then return >>; %% Now num must have the form y^d1 C1(x) + y^d2 C2(x), %% where d1 > d2 and d2 = 0 is allowed (if d <> 0 or d3 <> 0). if (C1 := get!!y!^n!*C(num_rhs, y)) then << d1 := car C1; C1 := cdr C1 >> else return; num_rhs := red num_rhs; %% Allow d2 = 0 => num_rhs freeof y if not smember(y, num_rhs) then << d2 := 0; C2 := num_rhs >> else if red num_rhs then return else if (C2 := get!!y!^n!*C(num_rhs, y)) then << d2 := car C2; C2 := cdr C2 >> else return; %% Den must have the form C3(x) or y^d3 C3(x). %% In the latter case, combine the powers of y. if smember(y, den_rhs) then if null red den_rhs and (den_rhs := get!!y!^n!*C(den_rhs, y)) then << d3 := car den_rhs; den_rhs := cdr den_rhs; d1 := {'difference, d1, d3}; d2 := {'difference, d2, d3} >> else return; %% Simplify the degrees of y and find which is 1: if d then << d1 := {'plus, d1, d}; d2 := {'plus, d2, d} >>; d1 := aeval d1; d2 := aeval d2; if d1 = 1 then << P := C1; Q := C2; n := d2 >> else if d2 = 1 then << P := C2; Q := C1; n := d1 >> else return; %% A final check that P, Q, n are valid: if Bernoulli!-depend!-check(P, y) or Bernoulli!-depend!-check(Q, y) or Bernoulli!-depend!-check(den_rhs, y) or Bernoulli!-depend!-check(n, x) then return; %% ( Last test implies Bernoulli!-depend!-check(n, y). ) %% For testing: %% return {'list, mk!*sq(P ./ den_rhs), mk!*sq(Q ./ den_rhs), n}; P := mk!*sq(P ./ den_rhs); Q := mk!*sq(Q ./ den_rhs); return ODESolve!-Bernoulli1(P, Q, y, x, n) end ) where depl!* = depl!*$ symbolic procedure Bernoulli!-depend!-check(f, xy); %% f is a standard form, xy is an identifier (kernel). if numr difff(f, xy) then << % neq 0 (nil) traceode "It is not of Bernoulli type because ..."; MsgPri(nil, !*f2a f, "depends (possibly implicitly) on", get(xy, 'odesolve!-depvar) or xy, nil); %% (y might be gensym -- 'odesolve!-depvar set in odeintfc) t >>$ symbolic procedure get!!y!^n!*C(u, y); %% Assume that u is a standard form representation of y^n * C(x) %% with y the leading kernel. Return (n . C) or nil begin scalar n, C; if mvar u eq y then << % ?? y^n * C(x), n nonnegint n := ldeg u; C := lc u; return if not domainp C and smember(y, mvar C) then ( if (C := get!!y!^n!*C1(C, y)) then % y^n * y^n1 * C(x), n nonnegint {'plus, n, car C} . cdr C ) else n . C >> else % (y^n1)^n * C(x), n nonnegint return get!!y!^n!*C1(u, y) end$ symbolic procedure get!!y!^n!*C1(u, y); % u = (y^n1)^n * C(x), n nonnegint begin scalar n, C; n := mvar u; if not(eqcar(n, 'expt) and cadr n eq y) then return; n := {'times, caddr n, ldeg u}; C := lc u; return n . C end$ algebraic procedure ODESolve!-Bernoulli1(P, Q, y, x, n); begin scalar !*odesolve_noint; % Force integration? traceode "It is of Bernoulli type."; n := 1 - n; return if symbolic !*odesolve_explicit then { y = ODENon!-Linear1PQ(n*P, n*Q, x)^(1/n)* newroot_of_unity(n) } % explicit form else { y^n = ODENon!-Linear1PQ(n*P, n*Q, x) } % implicit form end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Riccati ODEs algebraic procedure ODESolve!-Riccati(rhs, y, x); %% The ODE has the form df(y,x) = rhs. If rhs has the Riccati form %% a(x)y^2 + b(x)y + c(x) then transform to a reduced linear %% second-order ODE and attempt to solve it. symbolic if not !*odesolve_fast then % heuristic solution algebraic begin scalar a, b, c, soln, !*ratarg; traceode1 "Testing for a Riccati ODE ..."; %% rhs may have a denominator that depends on y, so ... symbolic(!*ratarg := t); c := coeff(rhs, y); if length c neq 3 or depends(c, y) then return; a := third c; b := second c; c := first c; c := a*c; b := -(df(a,x)/a + b); traceode "It is of Riccati type ", "and transforms into the linear second-order ODE: ", num(df(y,x,2) + b*df(y,x) + c*y) = 0; soln := {c, b, 1}; % low .. high soln := ODESolve!-linear!-basis(soln, 0, y, x, 2, if c then 0 else if b then 1 else 2); % min_order if not soln then << traceode "But ODESolve cannot solve it!"; return >>; %% The solution of the linear second-order ode must have the %% form y = arbconst(1)*y1 + arbconst(2)*y2, in which only the %% ratio of the arbconst's is relevant here, so ... soln := first soln; soln := newarbconst()*first soln + second soln; return {y = sub(y = soln, -df(y,x)/(a*y))} %% BEWARE: above minus sign missing in Zwillinger, first edn. end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ODEs that are "solvable for y or x", including Clairaut and Lagrange % as special cases. algebraic procedure ODESolve!-Clairaut(ode, ode_p, p, y, x); %% Assuming that ode is first order, determine whether it is of %% Clairaut type f(xy'-y) = g(y'), and if so return its general %% solution together with a singular solution if one exists. begin scalar sing_soln; traceode1 "Testing for a Clairaut ODE ..."; sing_soln := sub(y = x*p, ode_p); if depends(sing_soln, x) or depends(sing_soln, y) then return; % Not Clairaut traceode "It is of Clairaut type."; %% Look for a singular solution: sing_soln := solve(df(ode, x), df(y,x)); sing_soln := for each soln in sing_soln join %% NB: Cannot use algebraic code to check soln because it may %% contain derivatives that evaluate to 0. if symbolic eqcar(caddr soln, 'root_of) then {} else {num sub(soln, ode) = 0}; return (sub(p = newarbconst(), ode_p) = 0) . sing_soln end$ symbolic operator ODESolve!-Solvable!-y$ symbolic procedure ODESolve!-Solvable!-y(ode_p, p, y, x); %% ode_p has the form f(x,y,p) where p = y'. ( algebraic begin scalar c, lagrange, ode_y, ode_x; traceode1 "Testing for a ""solvable for y"" or Lagrange ODE ..."; %% Is the ODE "solvable for y"? c := coeff(ode_p, y); if length c neq 2 or depends(c, y) then return; ode_y := -first c / second c; %% ode_y has the form f(x,p) where p = y' and ode is y = ode_y. %% If f(x,p) = xF(p) + G(p) then ode is a Lagrange equation. if not depends(den ode_y, x) and length(c:=coeff(num ode_y, x)) = 2 and not depends(c, x) then lagrange := 1; % Lagrange form symbolic depend1(p, x, t); ode_x := num(p - df(ode_y, x)); if lagrange then << %% ODE is a Lagrange equation, hence ode_x is a LINEAR ODE %% for x(p) that can be solved EXPLICITLY (using an %% integrating factor) for a single value of x. symbolic depend1(x, p, t); ode_x := num(ode_x where df(p,x) => 1/df(x,p)); symbolic depend1(p, x, nil); traceode "It is of Lagrange type and reduces to this ", "subsidiary ODE for x(y'): ", ode_x = 0; ode_x := ODENon!-Linear1(ode_x, x, p) >> else if symbolic !*odesolve_fast then return traceode "Sub-solver terminated: fast mode, no heuristics!" else << %% ode_x is an arbitrary first-order ODE for p(x), so ... traceode "It is ""solvable for y"" and reduces ", "to this subsidiary ODE for y'(x):"; ode_x := ODESolve!-FirstOrder(ode_x, p, x); if not ode_x then << traceode "But ODESolve cannot solve it!"; return >> >>; traceode "The subsidiary solution is ", ode_x, " and the main ODE can be solved parametrically ", "in terms of the derivative."; return if symbolic(!*odesolve_implicit or !*odesolve_explicit) then %% Try to eliminate p between ode_y and ode_x else fail. %% Assume that the interface code will try to actually solve %% this for y: Odesolve!-Elim!-Param(ode_y, y, ode_x, p, y) else %% Return a parametric solution {y(p), x(p), p}: if lagrange then % soln explicit for x for each soln in ode_x collect ODESolve!-Simp!-ArbParam sub(p=newarbparam(), {y = sub(soln, ode_y), soln, p}) else for each soln in ode_x join << %% Make solution as explicit as possible: soln := solve(soln, x); for each s in soln collect ODESolve!-Simp!-ArbParam sub(p=newarbparam(), if symbolic eqcar(caddr s, 'root_of) then {y=ode_y, sub(part(rhs s, 2)=x, part(rhs s, 1)), p} else {y=sub(s, ode_y), s, p}) >> end ) where depl!* = depl!*$ symbolic operator ODESolve!-Solvable!-x$ symbolic procedure ODESolve!-Solvable!-x(ode_p, p, y, x); %% ode_p has the form f(x,y,p) where p = y'. not !*odesolve_fast and % heuristic solution ( algebraic begin scalar c, ode_x, ode_y; traceode1 "Testing for a ""solvable for x"" ODE ..."; %% Is the ODE "solvable for x"? c := coeff(ode_p, x); if length c neq 2 or depends(c, x) then return; ode_x := -first c / second c; %% ode_x has the form f(y,p) where p = y' and ode is x = ode_x. symbolic depend1(p, y, t); ode_y := num(1/p - df(ode_x, y)); %% ode_y is an arbitrary first-order ODE for p(y), so ... traceode "It is ""solvable for x"" and reduces ", "to this subsidiary ODE for y'(y):"; ode_y := ODESolve!-FirstOrder(ode_y, p, y); if not ode_y then << traceode "But ODESolve cannot solve it! "; return >>; traceode "The subsidiary solution is ", ode_y, " and the main ODE can be solved parametrically ", "in terms of the derivative."; return if symbolic(!*odesolve_implicit or !*odesolve_explicit) then %% Try to eliminate p between ode_x and ode_y else fail. %% Assume that the interface code will try to actually solve %% this for y: Odesolve!-Elim!-Param(ode_x, x, ode_y, p, y) else for each soln in ode_y join << %% Return a parametric solution {y(p), x(p), p}: %% Make solution as explicit as possible: soln := solve(soln, y); for each s in soln collect ODESolve!-Simp!-ArbParam sub(p=newarbparam(), if symbolic eqcar(caddr s, 'root_of) then {sub(part(rhs s, 2)=y, part(rhs s, 1)), x=ode_x, p} else {s, x=sub(s, ode_x), p}) >> end ) where depl!* = depl!*$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The arbitrary parameters in solutions: algebraic operator arbparam$ algebraic (!!arbparam := 0)$ algebraic procedure newarbparam(); arbparam(!!arbparam:=!!arbparam+1)$ algebraic procedure Odesolve!-Elim!-Param(ode_y, y, ode_x, p, depvar); %% ode_y is an expression for y and ode_x is an rlist of odesolve %% solutions for x, both in terms of a parameter p. Return a list %% of equations corresponding to the equations in ode_x but with p %% eliminated with ode_y, or nil if this is not possible. %% depvar is the true dependent variable. begin scalar soln, result; if polynomialp(ode_y := num(y - ode_y), p) then << result := {}; while ode_x neq {} and polynomialp(soln := num !*eqn2a first ode_x, p) do << result := resultant(ode_y, soln, p) . result; ode_x := rest ode_x >>; if ode_x = {} then return Odesolve!-Tidy!-Implicit(result, y) >>; ode_y := solve(ode_y, p); %% solve here may return a one_of construct (zimmer (4) & (19)). %% I don't see why, but ... %% (expand_cases not defined until solve loaded.) ode_y := expand_cases ode_y; if not smember(root_of, ode_y) then << result := for each soln in ode_y join for each s in ode_x join if rhs s = 0 then % s is f(x,y) = 0 if (s:=sub(soln, num lhs s)) neq 0 then {num s} else {} else {num(sub(soln, rhs s) - x)}; % s is x = f(x,y) return Odesolve!-Tidy!-Implicit(result, depvar) >>; traceode "But cannot eliminate parameter ", "to make solution explicit." end$ algebraic procedure Odesolve!-Tidy!-Implicit(solns, depvar); %% Remove repeated and irrelevant factors from implicit solutions. for each soln in solns join for each fac in factorize soln join if smember(depvar, fac:=first fac) then {fac = 0} else {}$ switch odesolve_simp_arbparam$ % DEFAULT OFF. TEMPORARY? symbolic operator ODESolve!-Simp!-ArbParam$ symbolic procedure ODESolve!-Simp!-ArbParam u; %% Simplify arbparam expressions within parametric solution u %% (cf. ODESolve!-Simp!-ArbConsts) begin scalar !*precise, x, y, ss_x, ss_y, arbexprns_x, arbexprns_y, arbexprns, param; if not(rlistp u and length u = 4) then TypErr(u, "parametric ODE solution"); if not !*odesolve_simp_arbparam then return u; % TEMPORARY? x := lhs cadr u; y := lhs caddr u; if not(ss_x := ODESolve!-Structr(caddr cadr u, x, y, 'arbparam)) then return u; if not(ss_y := ODESolve!-Structr(caddr caddr u, x, y, 'arbparam)) then return u; ss_x := cdr ss_x; ss_y := cdr ss_y; arbexprns_x := for each s in cdr ss_x collect caddr s; arbexprns_y := for each s in cdr ss_y collect caddr s; if null(arbexprns := intersection(arbexprns_x, arbexprns_y)) then return u; arbexprns_x := cdr ss_x; ss_x := car ss_x; arbexprns_y := cdr ss_y; ss_y := car ss_y; arbexprns_x := for each s in arbexprns_x join if member(caddr s, arbexprns) then {s} else << ss_x := subeval{s, ss_x}; nil >>; arbexprns_y := for each s in arbexprns_y join if member(caddr s, arbexprns) then {s} else << ss_y := subeval{s, ss_y}; nil >>; traceode "Simplifying the arbparam expressions in ", u, " by the rewrites ..."; param := cadddr u; for each s in arbexprns_x do << %% s has the form ansj = "expression in arbparam(n)" traceode rhs s => param; %% Remove other occurrences of arbparam(n): ss_x := algebraic sub(solve(s, param), ss_x); %% Finally rename ansj as arbparam(n): ss_x := subeval{{'equal, cadr s, param}, ss_x} >>; for each s in arbexprns_y do << ss_y := algebraic sub(solve(s, param), ss_y); ss_y := subeval{{'equal, cadr s, param}, ss_y} >>; %% Try a further heuristic simplification: if smember(param, den ss_x) and smember(param, den ss_y) then begin scalar ss_x1, ss_y1; ss_x1 := algebraic sub(param = 1/param, ss_x); if smember(param, den ss_x1) then return; ss_y1 := algebraic sub(param = 1/param, ss_y); if smember(param, den ss_y1) then return; traceode "Simplifying further by the rewrite ", param => 1/param; ss_x := ss_x1; ss_y := ss_y1 end; return makelist {{'equal, x, ss_x}, {'equal, y, ss_y}, param} end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% symbolic operator Polynomialp$ %% symbolic procedure polynomialp(pol, var); %% %% Returns true if numerator of pol is polynomial in var. %% %% Assumes on exp, mcd. %% begin scalar kord!*; kord!* := {var}; %% pol := numr simp!* pol; %% while not domainp pol and mvar pol eq var and %% (domainp lc pol or not smember(var, mvar lc pol)) do %% pol := red pol; %% return not smember(var, pol) %% end$ symbolic procedure Polynomialp(pol, var); %% Returns true if numerator of pol is polynomial in var. %% Assumes on exp, mcd. Polynomial!-Form!-p(numr simp!* pol, !*a2k var)$ symbolic procedure Polynomial!-Form!-p(sf, y); %% A standard form `sf' is polynomial if each of its terms is %% polynomial: domainp sf or (Polynomial!-Term!-p(lt sf, y) and Polynomial!-Form!-p(red sf, y))$ symbolic procedure Polynomial!-Term!-p(st, y); %% A standard term `st' is polynomial if either (a) its %% leading power is polynomial and its coefficient is free of y, or (b) %% its leading power is free of y and its coefficient is polynomial: if tvar st eq y then not smember(y, tc st) else if not smember(y, tvar st) then Polynomial!-Form!-p(tc st, y)$ endmodule; end;