Artifact f682981ae7fe2843c5d7273553c3be204186790bb59f4259d419544380577cde:
- Executable file
r38/packages/odesolve/odelin.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: 26903) [annotate] [blame] [check-ins using] [more...]
module odelin$ % Simple linear ODE solver % F.J.Wright@Maths.QMW.ac.uk, Time-stamp: <14 September 2000> % Based in part on code by Malcolm MacCallum. % TO DO: % Polynomial solutions for polynomial coeffs. % Check the distinction between finding a solution technique for an % ODE that then fails internally to solve it (e.g. failing to solve % an auxiliary equation) and failing to find a solution technique. % (Should the former be handled by `odefailure'?) %% Techniques implemented %% ====================== %% First order (integrating factor) %% Constant coefficients %% Euler and shifted Euler %% Exact %% Trivial order reduction (dep var and low-order derivs missing) %% Second-order special function ODEs (module odespcfn) %% Notes: Overall factors are handled in most cases by making the ODE %% "monic". %% Internal representation %% ======================= %% A linear ode is represented by its list of coefficient functions %% (odecoeffs), driver term (driver), and dependent (y) and %% independent (x) variables. The maximum (ode_order) and minimum %% (min_order) derivative orders are included in the representation %% for efficiency/convenience. Its solution is represented as a basis %% for the solution space of the reduced ODE together with a %% particular integral of the full ODE. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% algebraic procedure odesolve(ode, y, x); %% %% Temporary definition for test purposes. %% begin scalar !*precise, solution; %% ode := num !*eqn2a ode; % returns ode as expression %% if (solution := ODESolve!-linear(ode, y, x)) then %% return solution %% else %% write "***** ODESolve cannot solve this ODE!" %% end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global '(ODESolve_Before_Lin_Hook ODESolve_After_Lin_Hook)$ algebraic procedure ODESolve!-linear(ode, y, x); %% MAIN LINEAR SOLVER %% Assumes ODE is an algebraically irreducible "polynomial" expression. begin scalar reduced_ode, auxvar, auxeqn, odecoeffs, first_arb, solution, driver; %% The following decomposition is needed FOR ALL linear ODEs. %% The DRIVER is the part of the ODE independent of y, such that %% the ODE can be expressed as REDUCED_ODE = DRIVER. %% driver := if part(ode, 0) = plus then %% -select(~u freeof y, ode) else 0; driver := -sub(y=0, ode); reduced_ode := ode + driver; auxvar := symbolic gensym(); %% df(y, x, n) => m^n, where m = auxvar auxeqn := sub(y=e^(auxvar*x), reduced_ode)/e^(auxvar*x); odecoeffs := coeff(auxeqn, auxvar); % low .. high traceode "This is a linear ODE of order ", high_pow, "."; first_arb := !!arbconst + 1; symbolic if not(solution := ODESolve!-linear!-basis (odecoeffs, driver, y, x, high_pow, low_pow) or (not !*odesolve_fast and %% Add a switch to control access to this thread? %% It is currently necessary for Zimmer (8). << traceode "But ODESolve cannot solve it using linear techniques, so ..."; %% NB: This will probably produce a NONLINEAR ODE! %% But, in desperation, try it anyway ... (ODESolve!-Interchange(ode, y, x) where !*odesolve_basis = nil) >>)) then return; %% Return solution as BASIS or LINEAR COMBINATION, assuming a %% SINGLE solution since the ODE is linear: return if symbolic !*odesolve_basis then % return basis if (part(solution, 1, 0) = equal) then if lhs first solution = y and % solution is explicit (auxeqn := ODESolve!-LinComb2Basis (rhs first solution, first_arb, !!arbconst)) then auxeqn else << write "***** Cannot convert nonlinear combination solution to basis!"; solution >> else solution else % return linear combination if part(solution, 1, 0) = list then {y = ODESolve!-Basis2LinComb solution} else solution end$ algebraic procedure ODESolve!-linear!-basis (odecoeffs, driver, y, x, ode_order, min_order); %% Always returns the solution in basis format. %% Called by ODESolve!-Riccati in odenon1. symbolic if ode_order = 1 then ODESolve!-linear1(odecoeffs, driver, x) else or( ODESolve!-Run!-Hook( 'ODESolve_Before_Lin_Hook, {odecoeffs,driver,y,x,ode_order,min_order}), ODESolve!-linearn (odecoeffs, driver, y, x, ode_order, min_order, nil), ODESolve!-Run!-Hook( 'ODESolve_After_Lin_Hook, {odecoeffs,driver,y,x,ode_order,min_order}))$ algebraic procedure ODESolve!-linear!-basis!-recursive (odecoeffs, driver, y, x, ode_order, min_order); %% Always returns the solution in basis format. %% Internal linear solver called recursively. symbolic if ode_order = 1 then ODESolve!-linear1(odecoeffs, driver, x) else ODESolve!-linearn (odecoeffs, driver, y, x, ode_order, min_order, t)$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solve a linear first-order ODE by using an integrating factor. % Based on procedure linear1 from module ode1ord by Malcolm MacCallum. algebraic procedure ODESolve!-linear1(odecoeffs, driver, x); %% Solve the linear ODE reduced_ode = driver, where %% reduced_ode = A(x)*(dy/dx + P(x)*y), driver = A(x)*Q(x). %% Uses Odesolve!-Int to optionally turn off final integration. begin scalar A, P, Q; A := second odecoeffs; P := first odecoeffs/A; Q := driver/A; return if P then % dy/dx + P(x)*y = Q(x) begin scalar intfactor, !*combinelogs; traceode "It is solved by the integrating factor method."; %% intfactor simplifies better if logs are combined: symbolic(!*combinelogs := t); P := (P where tan(~x) => sin(x)/cos(x)); intfactor := exp(int(P, x)); return if Q then { {1/intfactor}, Odesolve!-Int(intfactor*Q,x)/intfactor } else { {1/intfactor} } end else << % dy/dx = Q(x) traceode "It is solved by quadrature."; if Q then {{1}, Odesolve!-Int(Q, x)} else {{1}} >> end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Try to solve a linear ODE of order > 1. % Based on procedure linearn from module linearn by Malcolm MacCallum. % If the first integral of an exact ODE has constant coefficients, is % of Euler type or has trivally reducible order then so does the % original ODE. Also, trivial order reduction preserves constant % coefficients or Euler type, and the reduced ODE is not further % reducible. Hence, when the linear ODE solver is called recursively % to solve a first integral of an exact ODE or a trivially order- % reduced ODE, the argument `recursive' is used to avoid checking % again whether it has constant coefficients, is of Euler type or has % trivially reducible order. algebraic procedure ODESolve!-linearn (odecoeffs, driver, y, x, ode_order, min_order, recursive); %% Solve the linear ODE: reduced_ode = driver. begin scalar lcoeff, odecoeffs1, driver1, solution; %% Make the ODE "monic" as assumed by some solvers: %% (Note that this makes algebraic factorization largely %% irrelevant!) if (lcoeff := part(odecoeffs, ode_order+1)) = 1 then << odecoeffs1 := odecoeffs; driver1 := driver >> else << odecoeffs1 := for each c in odecoeffs collect c/lcoeff; % low .. high %% Could discard last element of odecoeffs1 because it must be 1! driver1 := driver/lcoeff >>; if recursive then goto a; %% Test for constant coefficients: if odecoeffs1 freeof x then return ODESolve!-LCC(odecoeffs1, driver1, x, ode_order); traceode "It has non-constant coefficients."; %% Test for Euler form: if (solution := ODESolve!-Euler(odecoeffs1, driver1, x, ode_order)) then return solution; %% Test for trivial order reduction. The result cannot have %% constant coeffs or Euler form, but it could be first order or %% exact or ... if min_order neq 0 and ode_order neq min_order and %% else would reduce to purely algebraic equation (solution := ODELin!-Reduce!-Order (odecoeffs, driver, y, x, ode_order, min_order)) then return solution; a: %% Non-trivial solution techniques for recursive calls ... %% Test for exact form - try monic then original form: if (solution := ODELin!-Exact(odecoeffs1, driver1, y, x, ode_order)) then return solution; if lcoeff neq 1 and (solution := ODELin!-Exact(odecoeffs, driver, y, x, ode_order)) then return solution; %% Add other methods here ... %% Null return implies failure. %% FINALLY, test for a second-order special-function equation: if ode_order = 2 and (solution := ODESolve!-Specfn(odecoeffs1, driver1, x)) then return solution end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Convert between basis and linear combination formats % Solution basis output format: { {B1, B2, ...}, PI } % where {Bi} is a basis for the reduced ODE and PI is a particular % intergral for the full ODE, which may be absent. % This corresponds to the linear-combination output format % y = ( for each B in {B1, B2, ...} sum newarbconst()*B ) + PI. % This agrees with Maple, e.g. Maple V Release 5 gives % > dsolve(diff(y(x),x) + y(x) = x, output=basis); % [[exp(-x)], -1 + x] % > dsolve(diff(y(x),x) + y(x) = x); % y(x) = x - 1 + exp(-x) _C1 algebraic procedure ODESolve!-Basis2LinComb solution; %% Convert basis { {B1, B2, ...}, PI } to linear combination: begin scalar lincomb; lincomb := for each B in first solution sum <<newarbconst()>>*B; %% << >> above is NECESSARY to force immediate evaluation! if length solution > 1 then lincomb := lincomb + second solution; return lincomb end$ algebraic procedure ODESolve!-LinComb2Basis (lincomb, first_arb, last_arb); %% Convert linear combination to basis { {B1, B2, ...}, PI }: ODESolve!-LinComb2Basis1({}, lincomb, first_arb, last_arb)$ algebraic procedure ODESolve!-LinComb2Basis1 (basis, lincomb, first_arb, last_arb); %% `basis' is a LIST of independent reduced_ode solutions. %% Algorithm is to recursively move components from lincomb to %% basis. begin scalar coeffs, C; C := arbconst last_arb; coeffs := coeff(lincomb, C); if high_pow > 1 or smember(C, coeffs) then return % cannot convert else if high_pow = 1 then << basis := second coeffs . basis; lincomb := first coeffs >>; %% else independent of this arbconst return if first_arb >= last_arb then { basis, lincomb } else ODESolve!-LinComb2Basis1 (basis, lincomb, first_arb, last_arb-1) end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solve a linear, constant-coefficient ODE. % Based on code by Malcolm MacCallum. % There are (at least) 4 ways to get the particular integral (P.I.) % for a given driving term on the right of the equation: % 1. The method of undetermined coefficients: this is similar to the % integrator in that for a given driving term one has to find the % functional form of the P.I. and then solve for the numerical % coefficients in it. Making it really general is as big a task as % rewriting the integrator. % 2. The method of variation of parameters: this expands the P.I. as a % sum of functions of `x' times the linearly independent solutions in % the complementary function (C.F.). % 3. Factorise the linear operator (done anyway for the C.F.) and then % apply for each root `m' the operation % ans := exp(m*x) * int(ans * exp(-m*x)) % This is a form of the "D-operator method". N.B. Some `m' are % complex. % 4. Use Laplace transforms (and some kind of table lookup for the % inverse transforms). % The current implementation first tries to use the "D-operator % method", but as soon as any integral fails to evaluate it switches % to "variation of parameters". algebraic procedure ODESolve!-LCC(odecoeffs, driver, x, ode_order); % Returns a solution basis or nil (if it fails). begin scalar auxvar, auxeqn, i, auxroots, solutions, PI; traceode "It has constant coefficients."; %% TEMPORARY HACK -- REBUILD AUXEQN: auxvar := symbolic gensym(); i := -1; auxeqn := for each c in odecoeffs sum c*auxvar^(i:=i+1); % First we solve for the auxiliary roots: auxroots := solve(auxeqn, auxvar); % and check the solution carefully: if ode_order neq (for each multi in multiplicities!* sum multi) then return traceode "But insufficient roots of auxiliary equation!"; solutions := auxroots; a: if lhs first solutions neq auxvar then return traceode "But auxiliary equation could not be solved!"; if (solutions := rest solutions) neq {} then goto a; % Now we find the complementary solution: solutions := ODESolve!-LCC!-CompSoln(auxroots, x); % Next the particular integral: if driver = 0 then return { solutions }; if not (PI := ODESolve!-LCC!-PI(auxroots, driver, x)) then %% (Cannot use `or' as an algebraic operator!) PI := ODESolve!-PI(solutions, driver, x); return { solutions, PI } end$ algebraic procedure ODESolve!-LCC!-CompSoln(auxroots, x); %% Construct the complimentary solution (functions) from the roots %% of the auxiliary equation for a linear ODE with constant %% coefficients. Pairs of complex conjugate roots are converted to %% real trigonometric form up to the minimum of their %% multiplicities (regardless of complex switch and parameters). %% `auxroots' is a list of equations with a temporary variable on %% the left and an auxilliary root on the right. The root %% multiplicities are stored as a list in the global variable %% `multiplicities!*'. `x' is the independent variable. begin scalar multilist, crootlist, ans, multi, imroot, exppart; %% crootlist will be a list of lists of the form %% {unpaired_complex_root, multiplicity}. multilist := multiplicities!*; crootlist := {}; ans := {}; for each root in auxroots do << root := rhs root; multi := first multilist; multilist := rest multilist; % Test for complex roots: imroot := impart!* root; if imroot = 0 then << exppart := exp(root*x); for j := 1 : multi do ans := (x**(j-1)*exppart) . ans >> else begin scalar conjroot, conjmulti; %% Cannot assume anything about the order of the roots in %% auxroots, so build a list of the complex roots found to %% avoid using complex conjugate pairs twice. conjroot := conj!* root; conjmulti := 0; %% Essentially do assoc followed by delete if found: crootlist := for each root in crootlist join if first root = conjroot then << conjmulti := second root; {} >> else {root}; if conjmulti then % conjugate pair found: begin scalar minmulti; exppart := exp(repart!* root*x); minmulti := min(multi, conjmulti); imroot := abs imroot; % to avoid spurious minus sign imroot := (imroot where abs ~x => x); % to avoid spurious abs! for j := 1 : minmulti do ans := (x**(j-1)*cos(imroot*x)*exppart) . (x**(j-1)*sin(imroot*x)*exppart) . ans; if multi neq conjmulti then << %% Skip this unlikely case if possible minmulti := minmulti + 1; exppart := exp(root*x); for j := minmulti : multi do ans := (x**(j-1)*exppart) . ans; exppart := exp(conjroot*x); for j := minmulti : conjmulti do ans := (x**(j-1)*exppart) . ans >> end else crootlist := {root, multi} . crootlist end >>; %% Finally include unpaired complex roots: for each root in crootlist do << exppart := exp(first root*x); multi := second root; for j := 1 : multi do ans := (x**(j-1)*exppart) . ans >>; return ans end$ % The following procedures process complex-valued expressions with % regard only to their EXPLICIT complexity, i.e. assuming that all % symbolic quantities are pure real. They need to work with the % complex switch both on and off, which is slightly tricky! algebraic(vars!-are!-real := {repart ~x => x, impart ~x => 0})$ algebraic procedure repart!* u; << u := repart u; u where vars!-are!-real >>$ algebraic procedure impart!* u; << u := impart u; u where vars!-are!-real >>$ algebraic procedure conj!* u; << u := conj u; u where vars!-are!-real >>$ algebraic procedure ODESolve!-LCC!-PI(auxroots, driver, x); % Try to construct a particular integral using the `D-operator % method'. Factorise the linear operator (done anyway for the % C.F.) and then apply for each root m the operation % ans := exp(m*x) * int(ans * exp(-m*x)); % N.B. Some m may be complex. % See e.g. Stephenson, section 21.8, p 410. % Returns nil if any integral cannot be evaluated. begin scalar exp_mx, multiplicities, multi; traceode "Constructing particular integral using `D-operator method'."; multiplicities := multiplicities!*; while driver and auxroots neq {} do << exp_mx := exp((rhs first auxroots)*x); driver := driver/exp_mx; multi := first multiplicities; while driver and multi >= 1 do << driver := int(driver, x); if driver freeof int then multi := multi - 1 else driver := 0 >>; driver := exp_mx*driver; auxroots := rest auxroots; multiplicities := rest multiplicities >>; if driver = 0 then traceode "But cannot evaluate the integrals, so ..." else return driver end$ algebraic procedure ODESolve!-PI(solutions, R, x); % Given a "monic" forced linear nth-order ODE % y^(n) + a_(n-1)(x)y^(n-1) + ... + a_1(x)y = R(x) % and a set of n linearly independent solutions {yi(x)} to the % unforced ODE, construct a particular solution of the forced ODE % in the form of a (single) integral representation by the method % of variation of parameters. begin scalar n; traceode "Constructing particular integral using `variation of parameters'."; return if (n := length solutions) = 2 then begin scalar y1, y2, W; y1 := first solutions; y2 := second solutions; %% The Wronskian, kept separate to facilitate tracing: W := trigsimp(y1*df(y2, x) - y2*df(y1, x)); traceode "The Wronskian is ", W; R := R/W; return -ode!-int(y2*R, x)*y1 + ode!-int(y1*R, x)*y2 end else begin scalar Wmat, ys, W, i; %% Construct the (square) Wronskian matrix of the solutions: Wmat := {ys := solutions}; for i := 2 : n do Wmat := (ys := for each y in ys collect df(y,x)) . Wmat; load_package matrix; % to define mat Wmat := list2mat reverse Wmat; %% The Wronskian (determinant), kept separate for tracing: W := trigsimp det Wmat; traceode "The Wronskian is ", W; R := R/W; i := 0; return for each y in solutions sum ode!-int(cofactor(Wmat, n, i:=i+1)*R, x) * y end end$ % This facility should be in the standard matrix package! symbolic operator list2mat$ symbolic procedure list2mat M; % Input: (list (list A B ...) (list C D ...) ...) % Output: (mat (A B ...) (C D ...) ...) 'mat . for each row in cdr M collect cdr row$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Special cases of non-constant coefficients: algebraic procedure ODESolve!-Euler(odecoeffs, driver, x, ode_order); %% Solve a (MONIC) ODE having (essentially) the form %% reduced_ode = x^n df(y,x,n) + ... + a_{n-1} x df(y,x) + a_n y = driver %% odecoeffs = {a_n, a_{n-1} x, ..., a_0 x^n} / (a_0 x^n) %% = {a_n/(a_0 x^n), a_{n-1}/(a_0 x^{n-1}), ..., a_1/(a_0 x), 1} begin scalar tmp, shift, i, c, solution; odecoeffs := reverse odecoeffs; % high .. low %% Check for possible Euler or "shifted Euler" form. %% Find second non-zero ode coeff: tmp := rest odecoeffs; i := 1; while first tmp = 0 do << tmp := rest tmp; i := i+1 >>; tmp := first tmp; % second non-zero ode coeff tmp := den tmp; % ax^i or a(x+b)^i tmp := reverse coeff(tmp, x); % high .. low if high_pow neq i then return; % not Euler if second tmp then << % "shifted Euler" shift := second tmp/(i*first tmp); % b driver := sub(x=x-shift, driver) % x -> x-b >>; tmp := {first odecoeffs}; i := 0; odecoeffs := rest odecoeffs; a: if odecoeffs neq {} then << c := first odecoeffs * (x+shift)^(i:=i+1); if not(c freeof x) then return; % not Euler tmp := c . tmp; odecoeffs := rest odecoeffs; go to a >>; odecoeffs := tmp; traceode "It is of the homogeneous (Euler) type ", if shift then "(with shifted coefficients) " else "", "and is reducible to a simpler ODE ..."; i := -2; tmp := for each c in odecoeffs sum << i := i + 1; c * for j := 0 : i product (x-j) >>; odecoeffs := coeff(tmp, x); % TEMPORARY HACK! driver := sub(x=e^x, driver*x^ode_order); solution := ODESolve!-LCC(odecoeffs, driver, x, ode_order); solution := sub(x=log x, solution); if shift then solution := sub(x=x+shift, solution); return solution end$ algebraic procedure ODELin!-Exact(P_list, driver, y, x, n); %% Computes a (linear) first integral if ODE is an exact linear %% n'th order ODE P_n(x) df(y,x,n) + ... + P_0(x) y = R(x). begin scalar P_0, C, Q_list, Q, const, soln, PI; P_0 := first P_list; P_list := reverse rest P_list; % P_n, ..., P_1 %% ODE is exact if C = df(P_n,x,n) - df(P_{n-1},x,{n-1}) + ... %% + (-1)^{n-1} df(P_1,x) + (-1)^n P_0 = 0. for each P in P_list do C := P - df(C,x); C := P_0 - df(C,x); % C = 0 if exact if C then return; Q_list := {}; for each P in P_list do Q_list := (Q := P - df(Q,x)) . Q_list; % Q_0, ..., Q_{n-1} driver := int(driver, x) + (const := symbolic gensym()); %% The first integral is the LINEAR (n-1)'th order ODE %% Q_{n-1}(x) df(y,x,n) + ... + Q_0(x) y = int(R(x),x). traceode "It is exact, and the following linear ODE of order ", n-1, " is a first integral:"; if symbolic !*trode then << C := y; soln := first Q_list*y + ( for each Q in rest Q_list sum Q*(C := df(C,x)) ); write soln = driver >>; %% Recurse on the order: C := Q_list; %% First-integral ODE must have min order 0, since input ODE was %% already order-reduced. soln := ODESolve!-linear!-basis!-recursive (Q_list, driver, y, x, n-1, 0); PI := second soln; % MUST exist since driver neq 0 PI := coeff(PI, const); % { real PI, extra basis fn } return if high_pow = 1 then if first PI then { second PI . first soln, first PI } else { second PI . first soln } else << %% This error should now be redundant! write "*** Internal error in ODELin!-Exact:", " cannot separate basis functions! "; write "(Probably caused by `noint' option.)"; soln >> end$ algebraic procedure ODELin!-Reduce!-Order (odecoeffs, driver, y, x, ode_order, min_order); %% If ODE does not explicitly involve y (and perhaps low order %% derivatives) then simplify by reducing the effective order %% (unless there is only one) and try to solve the reduced ODE to %% give a first integral. Applies only to ODEs of order > 1. begin scalar solution, PI; ode_order := ode_order - min_order; for ord := 1 : min_order do odecoeffs := rest odecoeffs; traceode "Performing trivial order reduction to give the order ", ode_order, " linear ODE with coefficients (low -- high): ", odecoeffs; solution := ODESolve!-linear!-basis!-recursive (odecoeffs, driver, y, x, ode_order, 0); if not solution then << traceode "But ODESolve cannot solve the reduced ODE! "; return >>; traceode "Solution of order-reduced ODE is ", solution; traceode "Restoring order, ", y => df(y,x,min_order), ", to give: ", df(y,x,min_order) = solution, " and re-solving ..."; %% = lin comb of fns of x, so just integrate min_order times: if length solution > 1 then % PI = particular integral PI := second solution; solution := append( for each c in first solution collect ODESolve!-multi!-int(c, x, min_order), %% and add min_order extra basis functions: for i := min_order-1 step -1 until 0 collect x^i ); return if PI then { solution, ODESolve!-multi!-int(PI, x, min_order) } else { solution } end$ endmodule$ end$