Artifact 840499cab0e268880719091adc657be192e3c809576b4c342e9f08110d725e1f:
- Executable file
r38/packages/odesolve/odeintfc.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: 35042) [annotate] [blame] [check-ins using] [more...]
module odeintfc$ % Enhanced ODE solver interface % F.J.Wright@Maths.QMW.ac.uk, Time-stamp: <30 October 2000> % Use: odesolve(ode, y, x, conds, options) % or dsolve(ode, y, x, conds, options) (cf. Maple) % The first argument must evaluate to an ODE or a list of ODEs. % Each dependent variable y may be either an identifier or an % operator, such as y(x), which is replaced by a new identifier % internally. If y(x) is specified and x is not specified then x is % extracted from y(x) (cf. Maple). Equations containing operators of % the form y(x) with different arguments x are trapped, currently as % an error, until differential-delay equation solving is implemented. % If a dependent variable (y) does not depend on the independent % variable (x) then y is automatically declared to depend on x, and a % warning message to this effect is output. Derivatives are not % evaluated until this dependence has been enforced. BUT NOTE THAT % THIS DOES NOT WORK IF THE FIRST ARGUMENT IS AN ASSIGNMENT! This is % because the assignment is performed BEFORE the ode solver takes % control. This is something of an inconsistency in the current % REDUCE algebraic processing model. % All arguments after the first are optional but the order must be % preserved. If the first argument is a list of ODEs then y is % expected to be a list of dependent variables. If x is specified % then y must also be specified (first). An empty list can be used as % a place-holder argument. If x and/or y are missing then they are % parsed out of the ODE. % Thus, possible argument combinations, each of which may optionally % be followed by conds, are: ode | ode, y | ode, y, x % Currently, conditions can be specified only for a single ODE. % If specified, conds must take the form of an unordered list of % (unordered lists of) equations with either y, x, or a derivative of % y on the left. A single list of conditions need not be contained % within an outer list. Combinations of conditions are allowed. % Conditions within one (inner) list all relate to the same x value. % For example: % Boundary conditions: % {{y=y0, x=x0}, {y=y1, x=x1}, ...} % Initial conditions: % {x=x0, y=y0, df(y,x)=dy0, ...} % Combined conditions: % {{y=y0, x=x0}, {df(y,x)=dy1, x=x1}, {df(y,x)=dy2, y=y2, x=x2}, ...} % Boundary conditions on the values of y at various values of x may % also be specified by replacing the variables by equations with % single values or matching lists of values on the right, of the form: % y = y0, x = x0 | y = {y0, y1, ...}, x = {x0, x2, ...} % The final argument may be one of the identifiers % implicit, explicit, laplace, numeric, series % specifying an option. The options "implicit" and "explicit" set the % switches odesolve_implicit and odesolve_explicit locally. The other % options specify solution techniques -- they are not yet implemented. % TO DO: % Improved condition code to handle eigenvalue-type BVPs. % Solve systems of odes, calling crack where appropriate %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % User interface % ============== put('odesolve, 'psopfn, 'odesolve!-eval)$ put('dsolve, 'psopfn, 'odesolve!-eval)$ % alternative (cf. Maple) listargp odesolve, dsolve$ % May have single list arg. symbolic procedure odesolve!-eval args; %% Establish a suitable global environment: (if !*div or !*intstr or !*factor or not !*exp or not !*mcd then NoInt2Int reval result else NoInt2Int result) where result = begin scalar !*evallhseqp, !*multiplicities, !*div, !*intstr, !*exp, !*mcd, !*factor, !*ifactor, !*precise, !*nopowers, !*algint, !*echo; %% Turn echo off to stop Win32 PSL REDUCE (only) %% outputting its trigsimp lap code at the end of %% odesolve.tst. (Don't ask!) !*evallhseqp := !*exp := !*mcd := t; return odesolve!-eval1 args end$ symbolic procedure odesolve(ode, y, x); %% Direct symbolic-mode interface equivalent to MAHM's original. %% Calls odesolve!-eval to ensure correct environment. odesolve!-eval{ode, y, x}$ global '(ODESolve!-tracing!-synonyms)$ ODESolve!-tracing!-synonyms := '(trode trace tracing)$ symbolic procedure odesolve!-eval1 args; %% args = (ode &optional y x conds) %% Parse variables from ode if necessary (like solve), %% automatically declare y to depend on x if necessary and %% optionally impose conditions on the general solution. %% Support for systems of odes partly implemented so far. ( begin scalar ode, system, y, x, yconds, xconds, conds, soln; if null args then RedErr "ODESolve requires at least one argument -- the ODE"; begin scalar df_simpfn, !*uncached; !*uncached := t; % Turn off simplification of df (in case y does not yet % depend on x) before evaluating args (which may be lists): df_simpfn := get('df, 'simpfn); put('df, 'simpfn, 'simpiden); args := errorset!*({'revlis, mkquote args}, t); put('df, 'simpfn, df_simpfn); if errorp args then error1(); args := car args end; ode := car args; args := cdr args; system := rlistp ode; %% Find dependent and independent variables: %% (rlistp is a smacro defined in rlisp.red) if args then << y := car args; if rlistp y then if null cdr y then % empty list - ignore y := 'empty else if rlistp cadr y or eqnp cadr y then y := nil % condition rlist else if system then y := makelist % rlist of dependent variables for each yy in cdr y collect !*a2k yy else MsgPri("ODESolve: invalid second argument", y, nil, nil, t) else if system then TypErr(y, "dependent var list") else if eqnp y then if cadr y memq 'output . ODESolve!-tracing!-synonyms then % option y := nil else << yconds := caddr y; y := !*a2k cadr y >> else if not smember(y:=!*a2k y, ode) then y := nil % option >>; if y then args := cdr args; if args and y then << x := car args; if rlistp x then if null cdr x then % empty list - ignore x := 'empty else x := nil % condition list else if eqnp x then if cadr x memq 'output . ODESolve!-tracing!-synonyms then % option x := nil else << xconds := caddr x; x := !*a2k cadr x >> else if not smember(x:=!*a2k x, ode) then x := nil % option >>; if x then args := cdr args; if y eq 'empty then y := nil; if x eq 'empty then x := nil; %% If x not given and y an operator (list) then extract x: if null x and y then if rlistp y then begin scalar yy; yy := cdr y; while yy and atom car yy do yy := cdr yy; if yy and cdar yy then x := cadar yy; if not idp x then x := nil end else if pairp y and cdr y then x := cadr y; %% Finally, attempt to parse variables from ode if necessary: if null y or null x then %%%%% NOTE: ODE ALREADY AEVAL'ED ABOVE %%%%% so some of the following is now redundant !!!!! begin scalar df_simpfn, k_list, df_list; % Turn off simplification of df (in case y does not yet % depend on x) before evaluating ode, which may be a list: df_simpfn := get('df, 'simpfn); put('df, 'simpfn, 'simpiden); k_list := errorset!*({'get_k_list, mkquote ode}, t); put('df, 'simpfn, df_simpfn); if errorp k_list then error1() else k_list := car k_list; df_list := get_op_knl('df, car k_list); for each knl in cdr k_list do df_list := union(df_list, get_op_knl('df, knl)); %% df_list is set of derivatives in ode(s). if null df_list then RedErr "No derivatives found -- use solve instead."; %% df_list = ((df y x ...) ... (df z x ...) ... ) if null y then << y := cadar df_list . nil; for each el in cdr df_list do if not member(cadr el, y) then y := cadr el . y; %% y is a list at this point. if system then if length ode < length y then RedErr "ODESolve: under-determined system of ODEs." else y := makelist y % algebraic list of vars else if cdr y then MsgPri("ODESolve -- too many dependent variables:", makelist y, nil, nil, t) else y := car y; % single var MsgPri("Dependent var(s) assumed to be", y, nil, nil, nil) >>; if null x then << x := caddar df_list; MsgPri("Independent var assumed to be", x, nil, nil, nil) >>; end; %% Process the ode (re-simplifying derivatives): EnsureDependency(y, x); ode := aeval ode; %% !*eqn2a is defined in alg.red if system then if length ode > 2 then %% RedErr "Solving a system of ODEs is not yet supported." %% Skip conditions TEMPORARILY! return ODESolve!-Depend( makelist for each o in cdr ode collect !*eqn2a o, y, x, nil) else << ode := !*eqn2a cadr ode; y := cadr y >> else ode := !*eqn2a ode; %% Process conditions (re-simplifying derivatives): if args then if rlistp(conds := aeval car args) then << args := cdr args; conds := if not rlistp cadr conds then conds . nil else cdr conds >> else conds := nil; %% Now conds should be a lisp list of rlists (of equations). if yconds then yconds := if rlistp yconds then cdr yconds else yconds . nil; if xconds then xconds := if rlistp xconds then cdr xconds else xconds . nil; %% Concatenate separate x & y conds onto conds list: while yconds and xconds do << conds := {'list, {'equal, x, car xconds}, {'equal, y, car yconds}} . conds; yconds := cdr yconds; xconds := cdr xconds >>; if yconds or xconds then RedErr "Different condition list lengths"; if conds then %% Move this into odesolve!-with!-conds? conds := makelist odesolve!-sort!-conds(conds, y, x); %% Process remaining control option arguments: while args do begin scalar arg; arg := car args; args := cdr args; if eqnp arg then % equation argument if cadr arg eq 'output then args := caddr arg . args else if cadr arg memq ODESolve!-tracing!-synonyms then !*trode := caddr arg else MsgPri("Invalid ODESolve option", arg, "ignored.", nil, nil) % keyword argument else if arg memq '(implicit explicit expand noint verbose basis noswap norecurse fast check) then set(mkid('!*odesolve_, arg), t) else if arg eq 'algint then on1 'algint else if arg eq 'full or !*odesolve_full then !*odesolve_expand := !*odesolve_explicit := t else if arg memq ODESolve!-tracing!-synonyms then !*trode := t else if arg memq '(laplace numeric series) then RedErr{"ODESolve option", arg, "not yet implemented."} %% Pass remaining args to routine called else RedErr{"Invalid ODESolve option", arg} end; if !*odesolve_verbose then algebraic << write "ODE: ", num ode=0; write "Dependent variable: ", y, "; independent variable: ", x; write "Conditions: ", symbolic(conds or "none"); >>; %% Rationalize conflicting options: %% Conditions override basis if conds then !*odesolve_basis := nil; %% %% Basis overrides explicit %% if !*odesolve_basis then !*odesolve_explicit := nil; %% Finally, solve the ode! if not getd 'ODESolve!*0 then % for testing return {'ODESolve, ode, y, x, conds}; %% soln := if conds then %% odesolve!-with!-conds(ode, y, x, conds) %% else odesolve!-depend(ode, y, x); if null(soln := ODESolve!-Depend(ode, y, x, conds)) then return algebraic {num ode=0}; %% Done as follows because it may be easier to solve after %% imposing conditions than before, and it would be necessary to %% remove root_of's before imposing conditions anyway. if !*odesolve_explicit and not ODESolve!-basisp soln then soln := ODESolve!-Make!-Explicit(soln, y, conds); if !*odesolve_expand then soln := expand_roots_of_unity soln; if !*odesolve_check then ODE!-Soln!-Check(if !*odesolve_noint then NoInt2Int soln else soln, ode, y, x, conds) where !*noint = t; return soln end ) where !*odesolve_implicit = !*odesolve_implicit, !*odesolve_explicit = !*odesolve_explicit, !*odesolve_expand = !*odesolve_expand, !*trode = !*trode, !*odesolve_noint = !*odesolve_noint, !*odesolve_verbose = !*odesolve_verbose, !*odesolve_basis = !*odesolve_basis, !*odesolve_noswap = !*odesolve_noswap, !*odesolve_norecurse = !*odesolve_norecurse, !*odesolve_fast = !*odesolve_fast, !*odesolve_check = !*odesolve_check$ symbolic procedure Odesolve!-Make!-Explicit(solns, y, conds); << %% SHOULD PROBABLY CHECK THAT Y IS NOT INSIDE AN UNEVALUATED %% INTEGRAL BEFORE TRYING TO SOLVE FOR IT -- IT SEEMS TO UPSET %% SOLVE! solns := for each soln in cdr solns join if cadr soln eq y then {soln} else << %% soln is an implicit solution of ode for y %% for each s in cdr reval aeval {'solve, soln, y} join %% %% Make this test optional? %% if eqcar(caddr s, 'root_of) or eval('and . %% mapcar(cdr expand_roots_of_unity subeval{s, ode}, %% 'zerop)) %% then {s} %% else if !*trode then algebraic write "Solution ", s, %% " discarded -- does not satisfy ODE"; traceode "Solution before trying to solve for dependent variable is ", soln; cdr reval aeval {'solve, soln, y} >>; %% It is reasonable to return root_of's here. %% Solving can produce duplicates, so ... %% solns := union(solns, nil); % union still necessary? %% Check that each explicit solution still satisfies any %% conditions: if conds then for each cond in cdr conds do % each cond is an rlist begin scalar xcond; xcond := cadr cond; cond := makelist for each c in cddr cond collect !*eqn2a c; solns := for each s in solns join if eqcar(caddr s, 'root_of) or union(cdr %% trig_simplify subeval{xcond, subeval{s, cond}}, nil) = {0} then {s} else algebraic traceode "Solution ", s, " discarded -- does not satisfy conditions"; end; makelist solns >>$ % Should now be able to use the standard package `trigsimp' instead! algebraic procedure trig_simplify u; u where tan_half_angle_rules$ algebraic(tan_half_angle_rules := { sin(~u) => 2tan(u/2)/(1+tan(u/2)^2), cos(~u) => (1-tan(u/2)^2)/(1+tan(u/2)^2) })$ %% Cannot include tan rule -- recursive! symbolic procedure get_k_list ode; %% Return set of all top-level kernels in ode or rlist of odes. %% (Do not cache to ensure derivatives are [eventually] evaluated %% properly!) begin scalar k_list, !*uncached; !*uncached := t; %% Do not make an assignment twice: if eqcar(ode, 'setk) then ode := caddr ode; if rlistp(ode := reval ode) then << k_list := get_k_list1 cadr ode; for each el in cddr ode do k_list := union(k_list, get_k_list1 el) >> else k_list := get_k_list1 ode; return k_list end$ symbolic procedure get_k_list1 ode; union(kernels numr o, kernels denr o) where o = simp !*eqn2a ode$ symbolic procedure get_op_knl(op, knl); %% Return set of all operator kernels within knl with op as car. if pairp knl then if car knl eq op then knl . nil else ( if op_in_car then union(op_in_car, op_in_cdr) else op_in_cdr ) where op_in_car = get_op_knl(op, car knl), op_in_cdr = get_op_knl(op, cdr knl)$ symbolic procedure EnsureDependency(y, x); for each yy in (if rlistp y then cdr y else y . nil) do if not depends(yy, x) then << MsgPri("depend", yy, ",", x, nil); depend1(yy, x, t) >>$ symbolic procedure odesolve!-sort!-conds(conds, y, x); %% conds is a lisp list of rlists of condition equations. %% Return a canonical condition list. %% Collect conditions at the same value of x, check them for %% consistency and sort them by increasing order of derivative. begin scalar cond_alist; for each cond in conds do begin scalar x_cond, y_conds, x_alist; if not rlistp cond then TypErr(cond, "ode condition"); %% Extract the x condition: y_conds := for each c in cdr cond join if not CondEq(c, y, x) then TypErr(c, "ode condition equation") else if cadr c eq x then << x_cond := c; nil >> else c . nil; if null x_cond then MsgPri(nil, x, "omitted from ode condition", cond, t); if null y_conds then MsgPri(nil, y, "omitted from ode condition", cond, t); %% Build the new condition alist, with the x condition as key: if (x_alist := assoc(x_cond, cond_alist)) then nconc(x_alist, y_conds) else cond_alist := (x_cond . y_conds) . cond_alist end; %% Now cond_alist is a list of lists of equations, each %% beginning with a unique x condition. %% Sort the lists and return a list of rlists: return for each cond in cond_alist collect makelist if null cddr cond then cond else car cond . begin scalar sorted, next_sorted, this, next, result; sorted := sort(cdr cond, 'lessp!-deriv!-ord); %% sorted is a list of equations. while sorted and (next_sorted := cdr sorted) do << if cadr(this := car sorted) eq cadr(next := car next_sorted) then %% Two conds have same lhs, so ... ( if caddr this neq caddr next then MsgPri("Inconsistent conditions:", {'list, this, next}, "at", car cond, t) ) % otherwise ignore second copy else result := this . result; sorted := next_sorted >>; return reversip(next . result) end end$ symbolic procedure CondEq(c, y, x); %% Return true if c is a valid condition equation for y(x). %% cf. eqexpr in alg.red eqexpr c and ( (c := cadr c) eq x or c eq y or (eqcar(c, 'df) and cadr c eq y and caddr c eq x %% Is the following test overkill? and (null cdddr c or fixp cadddr c)) )$ symbolic procedure lessp!-deriv!-ord(a, b); %% (y=y0) < (df(y,x)=y1) and df(y,x,m)=ym < df(y,x,n)=yn iff m < n %% But y might be a kernel rather than an identifier! if atom(a := cadr a) then % a = (y=?) not atom cadr b % b = (df(y,x,...)=?) else if atom(b := cadr b) then % b = (y=?) not atom cadr a % a = (df(y,x,...)=?) else if not(car a eq 'df) then % a = (y(x)=?) car b eq 'df % b = (df(y(x),x,...)=?) else % a = (df(y,x,...)=?), any y car b eq 'df and % b = (df(y,x,...)=?) if null(a := cdddr a) then % a = (df(y,x)=?) cdddr b % b = (df(y,x,n)=?), 1 < n else % a = (df(y,x,m)=?), m > 1 (b := cdddr b) and car a < car b$ % b = (df(y,x,n)=?), m < n %%% THE FOLLOWING PROCEDURE SHOULD PROBABLY INCLUDE THE CODE TO MAKE %%% SOLUTIONS EXPLICIT BEFORE RESTORING OPERATOR FORMS FOR Y. symbolic procedure ODESolve!-Depend(ode, y, x, conds); %% Check variables and dependences before really calling odesolve. %% If y is an operator kernel then check whether ode is a %% differential-delay equation, and if not solve ode with y %% replaced by an identifier. ( begin scalar xeqt, ylist, sublist; y := if rlistp ode then cdr y else y . nil; %% Using `t' as a variable causes trouble when checking %% dependence of *SQ forms, which may contain `t' as their last %% element, so... if x eq t then << xeqt := t; x := gensym(); for each yy in y do if idp yy then depend1(yy, x, t); %% Cannot simply use `sub' on independent variables in %% derivatives, so... ode := subst(x, t, reval ode); % reval subst? if conds then conds := subst(x, t, reval conds); % reval subst? sublist := (t.x) . sublist >>; for each yy in y do if idp yy and not(yy eq t) then << %% Locally and quietly remove any spurious inverse %% implicit dependence of x on y: ylist := yy . ylist; depend1(x, yy, nil) where !*msg = nil; >> else % replace variable begin scalar yyy; yyy := gensym(); depend1(yyy, x, t); ylist := yyy . ylist; put(yyy, 'odesolve!-depvar, yy); % for later access sublist := (yy.yyy) . sublist; if xeqt then yy := subeval{{'equal,t,x}, yy}; odesolve!-delay!-check(ode, yy); ode := subeval{{'equal,yy,yyy}, ode}; if conds then conds := subeval{{'equal,yy,yyy}, conds} end; ylist := reverse ylist; ode := if rlistp ode then ODESolve!-System(cdr ode, ylist, x) else if conds then odesolve!-with!-conds(ode, car ylist, x, conds) else ODESolve!*0(ode, car ylist, x); if null ode then return; if sublist then begin scalar !*NoInt; %% Substitute into derivatives and integrals %% (and turn off integration for speed). ode := reval ode; % necessary? for each s in sublist do ode := subst(car s, cdr s, ode); %% ode := reval ode; % necessary? end; return ode end ) where depl!* = depl!*$ symbolic procedure ODESolve!-System(ode, y, x); %% TEMPORARY {'ODESolve!-System, makelist ode, makelist y, x}$ algebraic operator ODESolve!-System$ symbolic procedure odesolve!-delay!-check(ode, y); %% Check that ode is not a differential-delay equation in y, %% i.e. check that every occurrence of the operator y = y(x) has %% the same argument (without any shifts). This could be used as a %% hook to call an appropriate solver -- if there were one! begin scalar odelist; odelist := if rlistp ode then cdr ode else ode . nil; for each ode in odelist do ( for each knl in kernels numr simp ode do for each yy in get_op_knl(y_op, knl) do if not(yy eq y) then MsgPri("Arguments of", y_op, "differ --", "solving delay equations is not implemented.", t) ) where y_op = car y end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Impose initial/boundary conditions % ================================== % A first attempt to impose initial/boundary conditions on the % solution of a single ode returned by odesolve. % Solving with conditions provides access to the general solution as % the value of the global algebraic variable ode_solution. % If the solution is explicit then the following code could be % simplified and should be slightly more efficient, but is it worth % testing for an explicit solution and adding the special case code? algebraic procedure odesolve!-with!-conds(ode, y, x, conds); % conds must be a list of ordered lists of the form % {x=x0, y=y0, df(y,x)=y1, df(y,x,2)=y2, ...}. % All conditions applied at the same value of x must be collected % into the same list. % More generality is allowed only by odesolve!-eval. % This code could perhaps be more efficient by building a list of % all required derivatives of the ode solution once and for all? begin scalar first!!arbconst, arbconsts; first!!arbconst := !!arbconst + 1; %% Find the general solution of the ode and assign it to the %% global algebraic variable ode_solution: %1.03% ode_solution := odesolve!-depend(ode, y, x); ode_solution := symbolic ODESolve!*0(ode, y, x); if not ode_solution then return; traceode "General solution is ", ode_solution; traceode "Applying conditions ", conds; arbconsts := for i := first!!arbconst : !!arbconst collect arbconst i; return for each soln in ode_solution join odesolve!-with!-conds1(soln, y, x, conds, arbconsts) end$ algebraic procedure odesolve!-with!-conds1(soln, y, x, conds, arbconsts); begin scalar arbconsteqns; %% Impose the conditions (care is needed if the solution is %% implicit): arbconsteqns := for each cond in conds join begin scalar xcond, ycond, dfconds, arbconsteqns; xcond := first cond; cond := rest cond; ycond := first cond; if lhs ycond = y then cond := rest cond else ycond := 0; %% Now cond contains only conditions on derivatives. arbconsteqns := if ycond then % Impose the condition on y: {sub(xcond := {xcond, ycond}, soln)} else {}; dfconds := {}; %% Impose the conditions on df(y, x, n). If the solution %% is implicit, then in general all lower derivatives will %% be introduced, so ... while cond neq {} do begin scalar dfcond, result; %% dfcond : next highest derivative %% result : of substituting for all derivatives %% dfconds : all derivatives so far including this one dfcond := first cond; cond := rest cond; dfconds := dfcond . dfconds; %% All conditions on derivatives are handled before %% conditions on x and y to protect against %% substituting for x or y in df(y,x,...): result := sub(dfconds, map(y => lhs dfcond, soln)); if not(result freeof df) then % see comment below RedErr "Cannot apply conditions"; arbconsteqns := sub(xcond, result) . arbconsteqns end; return arbconsteqns end; %% Solve for the arbitrary constants: arbconsts := solve(arbconsteqns, arbconsts); %% ***** SHOULD CHECK THAT THE SOLUTION HAS SUCCEEDED! ***** %% and substitute each distinct arbconst solution set into the %% general ode solution: return for each cond in arbconsts collect if rhs soln = 0 then % implicit solution num sub(cond, lhs soln) = 0 else sub(cond, soln) end$ %% The above df error can happen only if the solution is implicit and %% a derivative is missing from the sequence, which is unlikely. %% Should try to recover by computing the value of the missing %% derivative from the conditions on lower order derivatives, and %% letting solve eliminate them. Try this later IF it ever proves %% necessary. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Check the solution set % ====================== % This code checks that the solution satisfies the ODE and that a % general solution is general. % A `solution' is either a basis solution: % {{b1(x), b2(x), ..., bn(x)}, PI(x)}; PI may be omitted if zero % or a list of component solutions. algebraic procedure ODE!-Soln!-Check(soln, ode, y, x, conds); %% SOLN is a LIST of solutions; ODE is a differential EXPRESSION; Y %% is the dependent variable; X is the independent variable; CONDS %% is true if conditions were specified. begin scalar n, !*allowdfint, !*expanddf; symbolic(!*allowdfint := !*expanddf := t); ode := num !*eqn2a ode; % returns ode as expression %% Should compute order `on demand' in ODE!-Comp!-Soln!-Fails. n := ODE!-Order(ode, y); if symbolic ODESolve!-basisp soln then << %% Basis solution (of linear ODE). %% Only arises as general solution. %% Remove contribution from PI if there is one: if arglength soln = 2 and second soln then ode := num sub(y = y + second soln, ode); if length(soln := first soln) neq n then write "ODESolve warning - ", "wrong number of functions in basis!"; %% Test each basis function in turn: for each s in soln do if (s:=sub(y = s, ode)) and trigsimp s then write "ODESolve warning - ", "basis function may not satisfy ODE: ", s >> else << %% List of component solutions. %% Check generality: if not conds and ODESolve!-arbconsts soln < n then write "ODESolve warning - ", "too few arbitrary constants in general solution!"; for each s in soln do if ODE!-Comp!-Soln!-Fails(s, ode, y, x, n) then write "ODESolve warning - ", "component solution may not satisfy ODE: ", s; >> end$ % Each component solution may be % explicit: y = f(x) % implicit: f(x,y) = g(x,y); rhs MAY be 0 % unsolved: ode = 0, but CAN THIS CASE ARISE? % parametric: {y = g(p), x = f(p), p} algebraic procedure ODE!-Comp!-Soln!-Fails(soln, ode, y, x, n); %% SOLN is a SINGLE component solution; ODE is a differential %% EXPRESSION; Y is the dependent variable; X is the independent %% variable; N is the order of ODE. if symbolic eqnp soln then % explicit, implicit or unsolved if lhs soln = y and rhs soln freeof y then % explicit: y = f(x) (if (ode := sub(soln, ode)) then trigsimp ode) else if rhs soln = 0 and lhs soln = ode then 1 % unsolved: ode = 0 else % implicit: f(x,y) = 0 begin scalar derivs, deriv; %% Construct in `derivs' a list of successive derivatives of %% the implicit solution f(x,y) up to the order of the ODE in %% decreasing order; each expression is linear in the highest %% derivative. derivs := {soln := num !*eqn2a soln}; for i := 1 : n do derivs := (soln:=num df(soln,x)) . derivs; %% Substitute for each derivative in ODE in turn in %% decreasing order until the result is zero; if not the %% solution fails. while n > 0 and << deriv := solve(first derivs, df(y,x,n)); % linear if deriv = {} then 0 else ode := num sub(first deriv, ode) >> do << n := n - 1; derivs := rest derivs >>; if deriv = {} then << write "ODESolve warning - cannot compute ", df(y,x,n); return 1 >>; derivs := first derivs; ode := (ode where derivs => 0); % for tracing return ode % 0 for good solution end else if symbolic(rlistp soln and eqnp cadr soln) then % parametric: {y = g(p), x = f(p), p} begin scalar xx, yy, p, dp!/dx, deriv, derivs; yy := rhs first soln; % Should not depend on ordering! xx := rhs second soln; % Should not depend on ordering! p := third soln; % parameter %% Construct in `derivs' a list of successive derivatives of the %% parametric solution (yy,xx) up to the order of the ODE in %% decreasing order. dp!/dx := 1/df(xx,p); derivs := {deriv:=yy}; for i := 1 : n do derivs := (deriv:=dp!/dx*df(deriv,p)) . derivs; %% Substitute for each derivative in ODE in turn in %% decreasing order until the result is zero; if not the %% solution fails. while n > 0 and (ode := num sub(df(y,x,n)=first derivs, ode)) do << n := n - 1; derivs := rest derivs >>; return sub(y=yy, x=xx, ode) end else write "ODESolve warning - invalid solution type: ", soln$ %% Code to find the actual number of arbitrary constants in a solution: fluid '(ODESolve!-arbconst!-args)$ symbolic operator ODESolve!-arbconsts$ symbolic procedure ODESolve!-arbconsts u; %% Return the number of distinct arbconsts in any sexpr u. begin scalar ODESolve!-arbconst!-args; ODESolve!-arbconsts1 u; return length ODESolve!-arbconst!-args end$ symbolic procedure ODESolve!-arbconsts1 u; %% Collect all the indices of arbconsts in u into a set in the %% fluid variable ODESolve!-arbconst!-args. if not atom u then if car u eq 'arbconst then (if not member(cadr u, ODESolve!-arbconst!-args) then ODESolve!-arbconst!-args := cadr u . ODESolve!-arbconst!-args) else << ODESolve!-arbconsts1 car u; ODESolve!-arbconsts1 cdr u >>$ endmodule$ end$