Artifact 67fd65a57f30ecd1da0c51680c4dd95000c260ba4a7d1dae268fa941e0703796:
- Executable file
r38/packages/odesolve/odepatch.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: 14246) [annotate] [blame] [check-ins using] [more...]
module odepatch$ % Patches to standard REDUCE facilities % F.J.Wright@Maths.QMW.ac.uk, Time-stamp: <18 September 2000> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Integrator patches % ================== % Avoid trying to integrate an integral to speed up some odesolve % calls. % NB: subst copies, even if no substitution is made. It is therefore % very likely to destroy uniqueness of kernels! %% load_package int$ %% apply1('load!-package, 'int)$ % not at compile time! packages_to_load int$ % not at compile time! global '(ODESolveOldSimpInt)$ (if not(s eq 'NoIntInt_SimpInt) then ODESolveOldSimpInt := s) where s = get('int, 'simpfn)$ % to allow reloading put('int, 'simpfn, 'NoIntInt_SimpInt)$ fluid '(!*NoInt !*Odesolve_NoInt)$ symbolic procedure NoIntInt_SimpInt u; %% Patch to avoid trying to re-integrate symbolic integrals in the %% integrand, because it can take forever and achieve nothing! if !*NoInt then begin scalar v, varstack!*; % Based on int/driver/simpint1. % Varstack* rebound, since FORMLNR use can create recursive % evaluations. (E.g., with int(cos(x)/x**2,x)). u := 'int . u; return if (v := formlnr u) neq u then << v := simp subst('int!*, 'int, v); remakesf numr v ./ remakesf denr v >> else !*kk2q u end else if atom u or null cdr u or cddr u and (null cdddr u or cddddr u) then rerror(int,1,"Improper number of arguments to INT") else if cddr u then simpdint u % header from simpint else begin scalar car_u, result; %% put('int, 'simpfn, 'SimpNoInt); car_u := mk!*sq simp!* car u; %% car_u := subeval{{'equal, 'Int, 'NoInt}, car_u}; car_u := subst('NoInt, 'Int, car_u); u := car_u . !*a2k cadr u . nil; %% Prevent SimpInt from resetting itself! put('int, 'simpfn, ODESolveOldSimpInt); % assumed & RESET by simpint result := errorset!*({ODESolveOldSimpInt, mkquote u}, t); put('int, 'simpfn, 'NoIntInt_SimpInt); % reset INT interface if errorp result then error1(); return NoInt2Int car result; %% Does this cause non-unique kernels? end$ algebraic operator NoInt$ % Inert integration operator %% symbolic procedure SimpNoInt u; %% !*kk2q('NoInt . u)$ % remain symbolic symbolic operator Odesolve!-Int$ symbolic procedure Odesolve!-Int(y, x); %% Used in SolveLinear1 on ode1 to control integration. if !*Odesolve_NoInt then formlnr {'NoInt, y, x} else mk!*sq NoIntInt_SimpInt{y, x}$ % aeval{'int, y, x}$ %% put('Odesolve!-Int, 'simpfn, 'Simp!-Odesolve!-Int)$ %% symbolic procedure Simp!-Odesolve!-Int u; %% %% Used in SolveLinear1 on ode1 to control integration. %% if !*Odesolve_NoInt then !*kk2q('NoInt . u) % must eval u!!! %% else NoIntInt_SimpInt u$ % aeval{'int, y, x}$ symbolic procedure NoInt2Int u; %% Convert all NoInt's back to Int's, without algebraic evaluation. subst('Int, 'NoInt, u)$ switch NoIntInt$ !*NoIntInt := t$ put('NoIntInt, 'simpfg, '((nil (put 'int 'simpfn 'SimpInt) (rmsubs)) (t (put 'int 'simpfn 'NoIntInt_SimpInt))))$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Differentiator patches % ====================== % Differentiate integrals correctly! % NB: `ON' is flagged ignore and so not compiled, so... on1 'allowdfint$ % To replace the versions in `reduce/packages/poly/diff.red' once % tested. % deflist('((dfint ((t (progn (on1 'allowdfint) (rmsubs)))))), 'simpfg); symbolic procedure diffp(u,v); % U is a standard power, V a kernel. % Value is the standard quotient derivative of U wrt V. begin scalar n,w,x,y,z; integer m; n := cdr u; % integer power. u := car u; % main variable. if u eq v and (w := 1 ./ 1) then go to e else if atom u then go to f %else if (x := assoc(u,dsubl!*)) and (x := atsoc(v,cdr x)) % and (w := cdr x) then go to e % deriv known. % DSUBL!* not used for now. else if (not atom car u and (w:= difff(u,v))) or (car u eq '!*sq and (w:= diffsq(cadr u,v))) then go to c % extended kernel found. else if x := get(car u,'dfform) then return apply3(x,u,v,n) else if x:= get(car u,dfn_prop u) then nil else if car u eq 'plus and (w := diffsq(simp u,v)) then go to c else go to h; % unknown derivative. y := x; z := cdr u; a: w := diffsq(simp car z,v) . w; if caar w and null car y then go to h; % unknown deriv. y := cdr y; z := cdr z; if z and y then go to a else if z or y then go to h; % arguments do not match. y := reverse w; z := cdr u; w := nil ./ 1; b: % computation of kernel derivative. if caar y then w := addsq(multsq(car y,simp subla(pair(caar x,z), cdar x)), w); x := cdr x; y := cdr y; if y then go to b; c: % save calculated deriv in case it is used again. % if x := atsoc(u,dsubl!*) then go to d % else x := u . nil; % dsubl!* := x . dsubl!*; % d: rplacd(x,xadd(v . w,cdr x,t)); e: % allowance for power. % first check to see if kernel has weight. if (x := atsoc(u,wtl!*)) then w := multpq('k!* .** (-cdr x),w); m := n-1; % Evaluation is far more efficient if results are rationalized. return rationalizesq if n=1 then w else if flagp(dmode!*,'convert) and null(n := int!-equiv!-chk apply1(get(dmode!*,'i2d),n)) then nil ./ 1 else multsq(!*t2q((u .** m) .* n),w); f: % Check for possible unused substitution rule. if not depends(u,v) and (not (x:= atsoc(u,powlis!*)) or not depends(cadddr x,v)) and null !*depend then return nil ./ 1; % Derivative of a dependent identifier via the chain rule. % Suppose u(v) = u(a(v),b(v),...), i.e. given depend {u}, a, % b, {a, b}, v; then (essentially) depl!* = ((b v) (a v) (u b % a)) if !*expanddf and not(v memq (x:=cdr atsoc(u, depl!*))) then << w := nil ./ 1; for each a in x do w := addsq(w, multsq(simp{'df,u,a},simp{'df,a,v})); go to e >>; w := list('df,u,v); w := if x := opmtch w then simp x else mksq(w,1); go to e; h: % Final check for possible kernel deriv. if car u eq 'df then << % multiple derivative if cadr u eq v then % (df (df v x y z ...) v) ==> 0 if commutedf if !*commutedf and null !*depend then return nil ./ 1 else if !*simpnoncomdf and (w:=atsoc(v, depl!*)) and null cddr w % and (cadr w eq (x:=caddr u)) then % (df (df v x) v) ==> (df v x 2)/(df v x) etc. % if single independent variable << x := caddr u; % w := simp {'quotient, {'df,u,x}, {'df,v,x}}; w := quotsq(simp{'df,u,x},simp{'df,v,x}); go to e >> else if eqcar(cadr u, 'int) then % (df (df (int F x) A) v) ==> (df (df (int F x) v) A) ? % Commute the derivatives to differentiate the integral? if caddr cadr u eq v then % Evaluating (df u v) where u = (df (int F v) A) % Just return (df F A) - derivative absorbed << w := 'df . cadr cadr u . cddr u; go to j >> else if !*allowdfint and % Evaluating (df u v) where u = (df (int F x) A) % (If dfint is also on then this will not arise!) % Commute only if the result simplifies: not_df_p(w := diffsq(simp!* cadr cadr u, v)) then << % Generally must re-evaluate the integral (carefully!) w := 'df . reval{'int, mk!*sq w, caddr cadr u} . cddr u; go to j >>; % derivative absorbed if (x := find_sub_df(w:= cadr u . merge!-ind!-vars(u,v), get('df,'kvalue))) then <<w := simp car x; for each el in cdr x do for i := 1:cdr el do w := diffsq(w,car el); go to e>> else w := 'df . w >> else if !*expanddf and not atom cadr u then << % Derivative of an algebraic operator u(a(v),...) via the % chain rule: df(u(v),v) = u_1(a(v),b(v),...)*df(a,v) + ... x := intern compress nconc(explode car u, '(!! !! !_)); y := cdr u; w := nil ./ 1; m := 0; for each a in y do begin scalar b; m:=m#+1; if numr(b:=simp{'df,a,v}) then << z := mkid(x, m); put(z, 'simpfn, 'simpiden); w := addsq(w, multsq(simp(z . y), b)) >> end; go to e >> else w := {'df,u,v}; j: if (x := opmtch w) then w := simp x else if not depends(u,v) and null !*depend then return nil ./ 1 else w := mksq(w,1); go to e end$ symbolic procedure dfform_int(u, v, n); % Simplify a SINGLE derivative of an integral. % u = '(int y x) [as main variable of SQ form] % v = kernel % n = integer power % Return SQ form of df(u**n, v) = n*u**(n-1)*df(u, v) % This routine is called by diffp via the hook % "if x := get(car u,'dfform) then return apply3(x,u,v,n)". % It does not necessarily need to use this hook, but it needs to be % called as an alternative to diffp so that the linearity of % differentiation has already been applied. begin scalar result, x, y, dx!/dv; y := simp!* cadr u; % SQ form integrand x := caddr u; % kernel result := if v eq x then y % Special case -- just differentiate the integral: % df(int(y,x), x) -> y replacing the let rule in INT.RED else if not !*intflag!* and % not in the integrator % If used in the integrator it can cause infinite loops, % e.g. in df(int(int(f,x),y),x) and df(int(int(f,x),y),y) !*allowdfint and % must be on for dfint to work << % Compute PARTIAL df(y, v), where y must depend on x, so % if x depends on v, temporarily replace x: result := if numr(dx!/dv:=diffp(x.**1,v)) then %% (Subst OK because all kernels.) subst(x, xx, diffsq(subst(xx, x, y), v)) where xx = gensym() else diffsq(y, v); !*dfint or not_df_p result >> then % Differentiate under the integral sign: % df(int(y,x), v) -> df(x,v)*y + int(df(y,v), x) addsq( multsq(dx!/dv, y), simp{'int, mk!*sq result, x}) % MUST re-simplify it!!! % (Perhaps I should use prepsq - % kernels are normally true prefix?) else !*kk2q{'df, u, v}; % remain unchanged if n neq 1 then result := multsq( (((u .** (n-1)) .* n) .+ nil) ./ 1, result); return result end$ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Solve patches % ============= % Needed for definition of `mkrootsoftag' function. packages_to_load solve$ % not at compile time! endmodule$ end$ % ***** NOT INSTALLED AT PRESENT ***** %% An algebraic solver appropriate to ODESolve, that never returns %% implicit solutions and returns nil when it fails. Also, it solves %% quadratics in terms of plus_or_minus. %% *** NB: This messes up `root_multiplicities'. *** %% (Later also use the root_of_unity operator where appropriate. %% Could do this by a wrapper around `solvehipow' in solve/solv1.) % This switch controls `solve' globally once this file has been % loaded: switch plus_or_minus$ % off by default %% The integrator does not handle integrands containing the %% `plus_or_minus' operator at all well. This may require some %% re-writing of the integrator (!). Temporarily, just turn off the %% use of the `plus_or_minus' operator when necessary. % This switch controls some `solve' calls locally within `ODESolve': switch odesolve_plus_or_minus$ % TEMPORARY % off by default -- it's to odangerous at present! % !*odesolve_plus_or_minus := t$ % TEMPORARY symbolic operator AlgSolve$ symbolic procedure AlgSolve(u, v); %% Return either a list of EXPLICIT solutions of a single scalar %% expression `u' for variable `v' or nil. begin scalar soln, tail, !*plus_or_minus; !*plus_or_minus := t; tail := cdr(soln := solveeval1{u, v}); while tail do if car tail eq v then tail := cdr tail else tail := soln := nil; return soln end$ algebraic procedure SolvePM(u, v); %% Solve a single scalar expression `u' for variable `v', using the %% `plus_or_minus' operator in the solution of quadratics. %% *** NB: This messes up `root_multiplicities'. *** symbolic(solveeval1{u, v} where !*plus_or_minus := !*odesolve_plus_or_minus)$ %% This is a modified version of a routine in solve/quartic symbolic procedure solvequadratic(a2,a1,a0); % A2, a1 and a0 are standard quotients. % Solves a2*x**2+a1*x+a0=0 for x. % Returns a list of standard quotient solutions. % Modified to use root_val to compute numeric roots. SLK. if !*rounded and numcoef a0 and numcoef a1 and numcoef a2 then for each z in cdr root_val list mkpolyexp2(a2,a1,a0) collect simp!* z else begin scalar d; d := sqrtq subtrsq(quotsqf(exptsq(a1,2),4),multsq(a2,a0)); a1 := quotsqf(negsq a1,2); return if !*plus_or_minus then % FJW list(subs2!* quotsq(addsq(a1, multsq(!*kk2q newplus_or_minus(),d)),a2)) %% Must uniquefy here until newplus_or_minus does it %% for itself! else list(subs2!* quotsq(addsq(a1,d),a2), subs2!* quotsq(subtrsq(a1,d),a2)) end$ endmodule$ end$