1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530253125322533253425352536253725382539254025412542254325442545254625472548254925502551255225532554255525562557255825592560256125622563256425652566256725682569257025712572257325742575257625772578257925802581258225832584258525862587258825892590259125922593259425952596259725982599260026012602260326042605260626072608260926102611261226132614261526162617261826192620262126222623262426252626262726282629263026312632263326342635263626372638263926402641264226432644264526462647264826492650265126522653265426552656265726582659266026612662266326642665266626672668266926702671267226732674267526762677267826792680268126822683268426852686268726882689269026912692269326942695269626972698269927002701270227032704270527062707270827092710271127122713271427152716271727182719272027212722272327242725272627272728272927302731273227332734273527362737273827392740274127422743274427452746274727482749275027512752275327542755275627572758275927602761276227632764276527662767276827692770277127722773277427752776277727782779278027812782278327842785278627872788278927902791279227932794279527962797279827992800280128022803280428052806280728082809281028112812281328142815281628172818281928202821282228232824282528262827282828292830283128322833283428352836283728382839284028412842284328442845284628472848284928502851285228532854285528562857285828592860286128622863286428652866286728682869287028712872287328742875287628772878287928802881288228832884288528862887288828892890289128922893289428952896289728982899290029012902290329042905290629072908290929102911291229132914291529162917291829192920292129222923292429252926292729282929293029312932293329342935293629372938293929402941294229432944294529462947294829492950295129522953295429552956295729582959296029612962296329642965296629672968296929702971297229732974297529762977297829792980298129822983298429852986298729882989299029912992299329942995299629972998299930003001300230033004300530063007300830093010301130123013301430153016301730183019302030213022302330243025302630273028302930303031303230333034303530363037303830393040304130423043304430453046304730483049305030513052305330543055305630573058305930603061306230633064306530663067306830693070307130723073307430753076307730783079308030813082308330843085308630873088308930903091309230933094309530963097309830993100310131023103310431053106310731083109311031113112311331143115311631173118311931203121312231233124312531263127312831293130313131323133313431353136313731383139314031413142314331443145314631473148314931503151315231533154315531563157315831593160316131623163316431653166316731683169317031713172317331743175317631773178317931803181318231833184318531863187318831893190319131923193319431953196319731983199 |
- comment
- Title: CRACK
- Authors:
- Andreas Brand
- Institut fuer Angewandte und Technische Informatik
- Friedrich Schiller Universitaet Jena
- 07740 Jena, Germany
- email: maa@hpux.rz.uni-jena.de
- tel.: + 49 3641 630790 ,
- Thomas Wolf
- School of Mathematical Sciences
- Queen Mary and Westfield College
- University of London
- London E1 4NS
- email: T.Wolf@maths.qmw.ac.uk
- tel.: + 44 71 975 5493
- Date of release: 25 August 93
- Abstract:
- CRACK is a package for solving overdetermined systems of partial or
- ordinary differential equations (PDEs, ODEs). Examples of programs
- which make use of CRACK for investigating ODEs (finding symmetries,
- first integrals, an equivalent Lagrangian or a "differential
- factorization") are added. The manual CRACKMAN gives further details.
- REDUCE version: 3.4.
- CRACK uses the package ODESOLVE of Malcolm MacCallum which is included
- in REDUCE 3.4.
- T. Wolf, An Analytic Algorithm for Decoupling and Integrating
- systems of Nonlinear Partial Differential Equations, J. Comp.
- Phys., no. 3, 60 (1985) 437-446.
- T. Wolf, The Symbolic Integration of Exact PDEs, preprint.
- M.A.H. MacCallum, An Ordinary Differential Equation Solver for REDUCE,
- Proc. ISAAC'88, Springer Lect. Notes in Comp Sci. 358,
- 196--205.
- Keywords: partial differential equations, computer analytic
- $
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- % %
- % CRACK Version 25 Aug. 93 %
- % %
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %**********************************************************************
- module crackstar$
- %**********************************************************************
- % Main program
- % Author: Andreas Brand
- % 1991,1993
- symbolic fluid '(REDUCEFUNCTIONS_ cont_ odesolve_ fcteval_ print_
- logoprint_ sp_cases independence_ tr_gensep decouple_
- tr_genint contradiction_ poly_only fname_ factorize_ nfct_
- fnew_ special_cases genint_ ineq_ time_)$
- symbolic operator setcrackflags$
- symbolic procedure setcrackflags$
- <<REDUCEFUNCTIONS_:= '(ABS ACOS ACOSD ACOSH ACOT ACOTD ACOTH ACSC ACSCD
- ACSCH ASEC ASECD ASECH ASIN ASIND ASINH ATAN
- ATAND ATANH ATAN2 ATAN2D CBRT COS COSD COSH COT
- COTD COTH CSC CSCD CSCH EXP FACTORIAL HYPOT LN
- LOG LOGB LOG10 SEC SECD SECH SIN SIND SINH SQRT
- TAN TAND TANH PLUS DIFFERENCE DF MINUS TIMES
- QUOTIENT EXPT INT)$
- if not fixp nfct_ then
- nfct_:=1$ % index of new functions and constants initialized
- cont_:=nil$ % if t then the user will be asked if a de
- % is too long for integration or substitution
- sp_cases:=nil$ % special cases for arbitrary functions or constants
- % are not of interest
- independence_:=nil$% if t then the user will be asked whether or not
- % expr. are considered to be linear independent
- genint_:=nil$ % generalized integration disabled/enabled
- factorize_:=5$ % recursion depth for factorizing equations
- odesolve_:=100$ % maximal length of a de (number of terms) to be
- % integrated as ode
- fcteval_:=100$ % maximal length of an expression to be substituted
- % for a function
- decouple_:=2$ % max. number of decoupling processes for a function
- tr_gensep:=nil$ % Trace generalized separation
- tr_genint:=nil$ % Trace generalized integration
- fname_:='c$ % name of new functions and constants (integration)
- logoprint_:=t$ % print logo for crack call
- poly_only:=t$ % all equations are polynomials only
- time_:=t$ % print the time needed for running crack
- print_:=8>>$ % maximal length of an expression to be printed
- symbolic operator crackhelp$
- symbolic procedure crackhelp$
- <<
- terpri()$
- write "The call: CRACK({expr1,...,exprn},{fct1,...,fctm},",
- "{var1,...vark})"$
- terpri()$terpri()$
- write "expr1,...,exprn are the equations expr1=0,...,exprn=0."$
- write "fct1,...,fctm are the unknown functions depending on ",
- "arguments"$
- terpri()$
- write " declared by DEPEND fcti,...$ ."$terpri()$
- write "var1,...,vark are further independent free variables which"$
- terpri()$
- write " none of the fcti depends on."$terpri()$terpri()$
- write "The Return: {res1,...,resr}, each resi representing a",
- " solution of the form"$terpri()$
- write "resi = {list_remain_ness_and_suff_conditions, ",
- "list_solved_functions,"$
- terpri()$
- write " list_free_or_unsolved_functions}"$terpri()$
- write "i.e. CRACK returns {} if it found that no solution exists"$
- terpri()$terpri()$
- write "The performance of CRACK can be controlled with special flags."$
- terpri()$
- write "The procedure SETCRACKFLAGS() sets the defaults for all flags"$
- terpri()$
- write "and the procedure LISTCRACKFLAGS() shows the flags and their ",
- "defaults."$
- terpri()$
- write "All flags are GLOBAL or FLUID symbolic mode variables."$
- terpri()$
- write "They can be changed with: SYMBOLIC (flagname:=value)$"$
- terpri()$terpri()$
- write "This helptext can be recalled by: CRACKHP()"$terpri()$
- write "For more detailed information see the manual."
- >>$
- symbolic operator listcrackflags$
- symbolic procedure listcrackflags$
- <<
- write
- "The values in parenthesis are the defaults. ",
- "For more details see the manual."$
- terpri()$terpri()$write
- "CONT_ (nil): ",
- "if t then the user will be asked if a de is too long"$
- terpri()$write
- " (wrt. odesolve_, fcteval_) ",
- "for integration or substitution"$
- terpri()$write
- "DECOUPLE_ (2): maximal number of decoupling attempts ",
- "for each function"$
- terpri()$write
- "FACTORIZE_ (5): recursion depth for factorizing equations"$
- terpri()$write
- "FCTEVAL_ (100): maximal length of an expression to be ",
- "substituted"$
- terpri()$write
- " for a function"$
- terpri()$write
- "FNAME_ ('c): name of new functions and constants (integr.)"$
- terpri()$write
- "GENINT_ (nil): generalized integration disabled/enabled"$
- terpri()$write
- "INDEPENDENCE_ (nil): if t then the user will be asked whether or not"$
- terpri()$write
- " expressions are considered to be linear ",
- "independent"$
- terpri()$write
- "ODESOLVE_ (100): maximal length of a de (number of terms) to be"$
- terpri()$write
- " integrated as ode"$
- terpri()$write
- "PRINT_ (8): maximal length of an expression to be printed"$
- terpri()$write
- "SP_CASES (nil): if nil then dropping also such factors after ",
- "factorization"$
- terpri()$write
- " which contain unknown functions or constants ",
- "which are"$
- terpri()$write
- "TIME_ (T): print the time needed for running crack"$
- terpri()$write
- " parametric and not to calculate"$
- terpri()$write
- "TR_GENINT (nil): Trace generalized integration"$
- terpri()$write
- "TR_GENSEP (nil): Trace generalized separation"$
- terpri()$
- >>$
- symbolic operator crack$
- symbolic procedure crack(el,il,fl,vl)$
- begin scalar l,l1,l2,n,m,ezgcdold,mcdold,gcdold,expold,ratold,ratargold$
- if print_ and logoprint_ then
- <<
- "**************************************************************",
- "************"$
- terpri()$
- write
- "This is CRACK - a solver for overdetermined partial differential ",
- "equations"$
- terpri()$
- write "Version 25-08-1993"$terpri()$
- write
- "***************************************************************",
- "***********"$
- terpri()$
- >>$
- ezgcdold:=!*ezgcd$
- gcdold:=!*gcd$
- expold:=!*exp$
- mcdold:=!*mcd$
- ratold:=!*rational$
- ratargold:=!*ratarg$
- !*ezgcd:=t$
- !*gcd:=nil$
- !*exp:=t$
- !*mcd:=t$
- %!*rational:=nil$
- !*rational:=t$
- !*ratarg:=t$
- fnew_:=nil$
- special_cases:=nil$
- vl:=union(reverse argset cdr fl,cdr vl)$
- % if print_ then
- % <<terpri()$write "CRACK is called with :"$terpri()$
- % terpri()$write "equations : "$deprint cdr el$terpri()$
- % terpri()$write "functions : "$fctprint cdr fl$terpri()$
- % terpri()$write "variables : "$fctprint vl$terpri()>>$
- n:=time()$
- m:=gctime()$
- l:=union(crack1(cdr el,nil,cdr il,cdr fl,vl,factorize_,nil,nil),nil)$
- if !*time or time_ then
- <<terpri()$write "CRACK needed : ",time()-n," ms GC time : ",
- gctime()-m," ms">>$
- !*ezgcd:=ezgcdold$
- !*gcd:=gcdold$
- !*exp:=expold$
- !*mcd:=mcdold$
- !*rational:=ratold$
- !*ratarg:=ratargold$
- l:=for each a in l collect
- <<l1:=nil$
- l2:=caddr a$
- for each b in cadr a do
- if (pairp b) and (car b = 'EQUAL) then l1:=cons(b,l1)
- else l2:=cons(b,l2)$
- list(car a,l1,l2)>>$
- return if l then
- cons('LIST,for each a in l collect list('LIST,cons('LIST,car a),
- cons('LIST,cadr a),
- cons('LIST,caddr a)))
- else list('LIST)
- end$
- symbolic procedure crack1(ss,sf,ineq,forg,vl,fact,level,ass)$
- % Main program
- begin scalar q,p,s,ftem,fges,f,ls,lf,l,l1,se,ndecouple,decfl$
- contradiction_:=nil$
- ineq_:=ineq$
- ndecouple:=0$
- ss:=desort ss$ % sort the de's w.r.t. its length
- sf:=desort sf$
- fges:=union(smemberl(fnew_,append(forg,append(ss,sf))),forg)$
- % add new functions
- fnew_:=nil$
- rep: ftem:=smemberl(fges,union(ss,sf))$ % only occuring functions
- if print_ and (ss or sf) then
- <<terpri()$write "equations: "$
- deprint(ss)$
- deprint(sf)>>$
- if print_ and ineq_ then
- <<terpri()$write "non-vanishing expressions: "$
- for each aa in ineq_ do eqprint aa>>$
- if print_ then
- <<terpri()$write "functions: "$
- fctprint(forg)>>$
- if print_ then fctprint(setdiff(ftem,forg))$
- for each p in union(sf,ss) do
- if member(p,ftem) then l:=cons(cons(0,p),l)$
- while l or ss do
- <<for each p in l do
- <<l1:=ls:=lf:=nil$
- if print_ then
- <<terpri()$write "substitution : "$
- substprint(cdr p,car p)>>$
- if ineq_ and null (ineq_:=ineqsubst(p,ftem)) then l:=ss:=nil$
- while not contradiction_ and ss do
- <<s:=car ss$
- ss:=cdr ss$
- ls:=cons(reval subst(car p,cdr p,s),ls)$
- contradiction(s,fges)>>$
- while not contradiction_ and sf do
- <<s:=car sf$
- sf:=cdr sf$
- if freeof(s,cdr p) then lf:=cons(s,lf)
- else
- <<ls:=cons(reval subst(car p,cdr p,s),ls)$
- if contradiction(s,fges) then contradiction_:=t>> >>$
- if not contradiction_ then
- <<for each s in se do
- if freeof(s,cdr p) then l1:=cons(s,l1)$
- se:=reverse l1$
- ss:=desort delete(0,union(ls,nil))$
- sf:=desort delete(0,union(lf,nil))>>$
- forg:=fctsubst(car p,cdr p,forg)$
- ftem:=delete(cdr p,ftem)$
- fges:=delete(cdr p,fges)$
- >>$ % alle Fkt. subst.
- ls:=ss$
- l:=nil$
- while not contradiction_ and ls do
- <<q:=car ls$
- ls:=cdr ls$
- ss:=delete(q,ss)$
- if pairp q then
- <<l1:=desort gensep(q,smemberl(union(fnew_,ftem),q),vl)$
- sf:=union(reverse l1,sf)$
- ftem:=union(fnew_,smemberl(ftem,append(sf,ss)))$
- fges:=union(fnew_,smemberl(fges,append(sf,ss)))$
- for each p in l1 do
- if member(p,ftem) then l:=union(list cons(0,p),l)
- else se:=cons(p,se)$
- if l then ls:=nil>>
- else se:=cons(q,se)$
- se:=desort union(se,nil)$
- if not contradiction_ then
- while se and not l do
- <<if q:=fcteval(car se,ftem,vl) then
- <<l:=union(list q,l)$
- ls:=nil>>$
- se:=cdr se>>
- >>$
- while se and not l do
- <<if q:=fcteval(car se,ftem,vl) then l:=union(list q,l)$
- se:=cdr se>>
- >>$
- if contradiction_ then RETURN NIL$
- if fact and (fact>0) then
- q:=splitandcrack(sf,ineq_,forg,ftem,vl,fact,level,ass)
- else q:=nil$
- decfl:=makedecfl(fges,decfl)$
- if q then p:=nil
- else % decoupling process starts
- if (ndecouple<decouple_) and
- (p:=decouple(sf,forg,ftem,fges,vl,decfl,fact)) then
- if contradiction_ then RETURN NIL
- else
- <<ndecouple:=add1 ndecouple$
- if setdiff(car p,sf) or setdiff(sf,car p) or (cadr p) then
- <<sf:=car p$
- ss:=cadr p$
- forg:=caddr p$
- fges:=cadddr p$
- decfl:=car cddddr p$
- goto rep>> >>$
- if not q then
- <<q:=list list(
- sf,
- forg,
- smemberl(setdiff(fges,forg),
- append(sf,for each l in equallist forg collect cddr l)))$
- if print_ then resout(caar q,cadar q,caddar q)>>$
- ineq_:=ineq$
- return q
- end$ % of crack1
- symbolic procedure ineqsubst(p,ftem)$
- % tests all q's in ineq_ for subst(car p, cdr p,q)=0
- % result: nil, if 0 occurs
- % otherwise list of the subst(car p,...)
- begin scalar l,a$
- while ineq_ do
- if not freeof(car ineq_,cdr p) then
- <<a:=reval subst(car p, cdr p,car ineq_)$
- if zerop a then
- <<if print_ then
- <<terpri()$write "contradiction :"$
- eqprint list('EQUAL,cdr p,car p)$
- write "because of non-vanishing expression:"$
- deprint list car ineq_>>$
- contradiction_:=t$
- l:=nil$
- ineq_:=nil>>
- else
- <<if pairp a and (car a='QUOTIENT) then a:=cadr a$
- a:=cdr reval list('FACTORIZE,a)$
- for each s in a do
- if smemberl(ftem,s) and not member(s,l) then l:=cons(s,l)$
- ineq_:=cdr ineq_>> >>
- else
- <<l:=cons(car ineq_,l)$
- ineq_:=cdr ineq_>>$
- return reverse l$
- end$
- symbolic procedure contradiction(s,fges)$
- % tests expression s for contradiction an non-zero expressions
- if (s=0) or smemberl(fges,s) then
- if member(s,ineq_) then
- <<if print_ then
- <<terpri()$write "non-zero expression vanishes : "$
- deprint list s>>$
- t>>
- else nil
- else
- <<if print_ then
- <<terpri()$write "Contradiction found : "$
- deprint list s>>$
- t>>$
- symbolic procedure makedecfl(fges,decfl)$
- % makes decouple flag list of functions in fges
- begin scalar s$
- while decfl do
- <<if member(caar decfl,fges) then s:=cons(car decfl,s)$
- decfl:=cdr decfl>>$
- for each f in fges do
- if not assoc(f,s) then s:=cons(cons(f,decouple_),s)$
- return s$
- end$
- symbolic procedure resout(sorg,forg,fnew)$
- begin scalar l$
- if print_ then
- <<for each a in forg do
- if (pairp a) and (car a='EQUAL) then l:=cons(a,l)
- else fnew:=cons(a,fnew)$
- terpri()$write "End of this CRACK run"$terpri()$terpri()$
- write "The solution : "$
- if l then
- fctprint l$
- if sorg then
- <<write "Remaining conditions : "$
- deprint sorg$
- if fnew then
- <<write "for the functions : "$fctprint fnew>> >>
- else if fnew then
- <<write "Free functions or constants : "$fctprint fnew>>$
- terpri()$
- write
- "*************************************************************",
- "*************"$terpri()$terpri()>>$
- end$
- symbolic procedure splitandcrack(sorg,ineq,forg,ftem,vl,fact,level,ass)$
- % If possible one de form sorg is factoirized and
- % crack1 is called with each factor instead
- begin scalar l,p,q,a,n$
- l:=sorg$
- while l do
- if (p:=splitde(car l,ftem)) then % factorizing a de
- <<a:=car l$
- if print_ then
- <<terpri()$
- write "factorizing "$deprint(list a)$
- write "we get the alternative equations : "$deprint(p)>>$
- l:=nil>>
- else l:=cdr l$
- if p then
- <<n:=0$
- for each d in p do % for each factor
- if not member(d,ineq_) then
- <<ass:=cons(d,ass)$
- n:=n+1$
- level:=cons(n,level)$
- if print_ then
- <<terpri()$
- write "Level : "$
- for each m in reverse level do write m,"."$terpri()$
- write "CRACK is now called with the assumtion(s) : "$
- deprint(reverse ass)>>$
- % sorg:=cons(d,delete(a,sorg))$
- l:=crack1(list(d),delete(a,sorg),ineq,
- union(forg,smemberl(ftem,sorg)),
- vl, if fact<2 then nil else fact-1,level,ass)$
- ineq:=union(list d,ineq)$
- ass:=cdr ass$
- level:=cdr level$
- if l then q:=union(clearphantomfct(l,forg),q)>> >>$
- % crack1 is called
- return q
- end$
- symbolic procedure clearphantomfct(l,forg)$
- % functions from higher level ftem beeing removed from forg
- for each a in l collect clrphfct(a,forg)$
- symbolic procedure clrphfct(a,forg)$
- begin scalar l,l1,fl$
- if a then
- <<for each f in forg do
- if (pairp f) and (car f='EQUAL) then fl:=cons(cadr f,fl)
- else fl:=cons(f,fl)$ % List of functions in forg
- for each f in cadr a do
- if (pairp f) and (car f='EQUAL) then
- (if member(cadr f,fl) then l:=cons(f,l)) % Member of forg
- else l1:=cons(f,l1)$ % free or unknown function
- l1:=smemberl(l1,append(car a,append(fl,l)))$
- % occuring in forg or sorg
- return list(car a,l,union(l1,caddr a))>>
- end$
- endmodule$
- %********************************************************************
- module decoupling$
- %********************************************************************
- % Routines for decoupling de's
- % Author: Andreas Brand
- % August 1991
- symbolic procedure eqderiv(p,f,vl)$
- % yields a list of de's by differentiating p w.r.t all those
- % variables from vl on which f not depends
- begin scalar l$
- vl:=setdiff(vl,fctargs f)$ % variables on which f not depends
- for each v in vl do
- if not freeof(p,v) then l:=cons(reval list('DF,p,v),l)$
- return l
- end$
- symbolic procedure decouplelist(liste,f)$
- % produces a list of de's for decoupling process w.r.t function f
- % the elements are of the type: ((expr (deriv.pot)) flag),
- % f Funktion
- begin scalar p,q$
- p:=nexteq(liste,t)$ % find first de marked with t
- q:=nexteq(liste,nil)$ % find first de marked with nil
- return if minausdp(p,q) then
- append(derivlist(list(list(q,nil),list(p,t)),f,
- difdiff(cadr q,cadr p),fctargs f,t),liste)
- % simplier de is marked with t
- % the other with nil
- % than the de's are diff.
- else append(derivlist(list(list(p,t),list(q,nil)),f,
- difdiff(cadr p,cadr q),fctargs f,nil),liste)
- end$
- symbolic procedure nexteq(liste,flag)$
- % finds the first element in liste marked with flag
- begin scalar stop$
- while (liste and (not stop)) do
- if flag=cadar liste then stop:=t
- else liste:=cdr liste$
- return if stop then caar liste
- else nil$
- end$
- symbolic procedure derivlist(dl,f,l,v,flag)$
- % differentiates de's from dl, dl:list of de's
- % f:function, l:list of numbers, v:list of variables,
- % flag:de to be substituted by a simplier one
- begin scalar m,p,q,li$
- for each x in v do
- <<m:=car l$
- l:=cdr l$
- if m<0 then % second de to be diff.
- <<if not q then q:=car nexteq(dl,not flag)$
- q:=reval list('DF,q,x,list('MINUS,m))$
- li:=cons(list(cons(q,ldiff(q,f)),not flag),li)>>
- else if 0<m then % first de to be differentiated
- <<if not p then p:=car nexteq(dl,flag)$
- p:=reval list('DF,p,x,m)$
- li:=cons(list(cons(p,ldiff(p,f)),flag),li)>> >>$
- return li$
- end$
- symbolic procedure fctchoose(stem,fl)$
- % finds function from fl with the maximal number of arguments
- if null (fl:=smemberl(fl,stem)) then nil
- else begin
- scalar s,n$
- s:=car fl$
- n:=fctlength s$
- for each f in cdr fl do
- if length(fctargs f)>n then
- <<s:=f$n:=fctlength s>>$
- return s
- end$
- symbolic procedure dechoose(gl,f)$
- % chooses two de's from list gl witch are optimal to be decoupled
- % w.r.t function f
- begin scalar d,c$
- d:=shortdesearch(gl,f)$
- c:=secondde(delete(d,gl),d,f)$
- return if null c then list(list(d,nil),list(nil,t))
- else
- if minausdp(d,c) then list(list(c,nil),list(d,t))
- else list(list(d,nil),list(c,t))
- end$
- symbolic procedure shortdesearch(gl,f)$
- % finds shortest de from list gl in which function f occurs
- if null gl then nil
- else if freeof(caar gl,f) then shortdesearch(cdr gl,f)
- else shortde(car gl,shortdesearch(cdr gl,f))$
- symbolic procedure shortde(g1,g2)$
- % finds the shorter of two de's
- if null g1 then g2
- else if null g2 then g1
- else begin scalar n,m$
- n:=delength(car g1)$
- m:=delength(car g2)$
- return if n<m then g1
- else if n=m then minausd(g1,g2)
- % the simplier de
- else g2$
- end$
- symbolic procedure secondde(gl,g,f)$
- % Auswahl der 2. Dgl. aus gl,g 1.Dgl.,f Funktion
- begin scalar c$
- c:=secondde1(gl,g,f)$
- if null c then c:=shortdesearch(gl,f)$
- return c
- end$
- symbolic procedure secondde1(gl,g,f)$
- % Auswahl der kuerzesten Gl. aus gl, die die Funktion f enthaelt und
- % einfacher als die Gleichung g ist
- if null gl then nil
- else
- if freeof(caar gl,f) or minausdp(car gl,g) then secondde1(cdr gl,g,f)
- % Gl., die f nicht enthalten oder
- % schwieriger als g sind, nicht
- else shortde(car gl,secondde1(cdr gl,g,f))$
- % die kuerzere aus der 1. Gl. und
- % der kuerzesten des Rests
- symbolic procedure decouple(sorg,forg,ftem,fges,vl,decfl,fact)$
- begin scalar repe,a,f,fl,l,l1,s,p,q,r,ss,s13,stem,ftem1,ftem2,stop2,h,
- fctflag$
- fl:=t$
- repe:=1$
- ftem1:=ftem$
- stem:=sorg$
- rep:
- stop2:=nil$
- ftem1:=smemberl(ftem1,stem)$
- ftem2:=nil$
- for each a in ftem1 do
- if assoc(a,decfl) then ftem2:=cons(a,ftem2)$
- if null ftem2 then
- <<l:=union(setdiff(stem,sorg),setdiff(s13,sorg))$
- % new equations
- stop2:=t$
- if null l or repe<=0 then
- <<for each a in l do
- sorg:=union(list reval aeval a,sorg)$
- sorg:=desort sorg$
- fl:=nil>>
- else <<repe:=sub1 repe$
- for each a in l do sorg:=union(list reval aeval a,sorg)$
- sorg:=desort sorg>> >>
- else if (f:=fctchoose(stem,ftem2)) then
- <<fctflag:=t$
- f:=assoc(f,decfl)$
- decfl:=cons(cons(car f,sub1 cdr f),delete(f,decfl))$
- f:=car f$
- if print_ then
- <<terpri()$write "decoupling: "$eqprint(f)>>$
- stem:=powappend(stem,f)$
- if print_ then <<terpri()$write "new equations: ">>$
- begin scalar stop3$
- rep1: stop3:=nil$
- l1:=dechoose(stem,f)$ % optimale Gl. suchen
- p:=caar l1$ % p schwierigere Gl.
- if (null caadr l1) and fctflag then
- if length varslist(p,smemberl(ftem,p),vl)>fctlength f then
- % f kommt nur in einer
- % Gl. vor und haengt nicht
- % von allen Variablen ab
- <<fctflag:=nil$
- stem:=union(powappend(eqderiv(car p,f,vl),f),stem)$
- % Hinzufuegen neuer Gl.
- l1:=dechoose(stem,f)$
- p:=caar l1>>$
- if not null caadr l1 then
- <<repeat
- <<l1:=decouplelist(l1,f)$ %list of de's to be decoupled
- h:=mkldiff(f,cadaar l1)$ % h leading derivative
- s:=car nexteq(l1,t)$
- if pairp s and (car s='QUOTIENT) then
- <<if (a:=casecheck(caddr s,vl)) then
- ineq_:=union(list a,ineq_)$
- s:=cadr s$>>$
- r:=car nexteq(l1,nil)$
- if pairp r and (car r='QUOTIENT) then
- <<if (a:=casecheck(caddr r,vl)) then
- ineq_:=union(list a,ineq_)$
- r:=cadr r>>$
- s:=reval !*q2a simpresultant list(s,r,h)$ % alg. decoupl.
- s:=car simplifyde(s,ftem1,vl,nil)$
- if contradiction(s,ftem1) then
- <<stem:=nil$
- stop2:=t$
- stop3:=t$
- contradiction_:=t$
- fl:=nil>>
- else
- <<ftem1:=union(fnew_,ftem1)$
- fges:=union(fnew_,fges)$
- if freeof(s,f) and
- stardep(s,smemberl(ftem1,s),varslist(s,ftem1,vl)) then
- <<ss:=list s$
- sorg:=desort sorg$
- fges:=union(fnew_,fges)$
- stem:=nil$
- stop2:=t$
- stop3:=t>> >>$
- if not stop3 then
- if (q:=fcteval(s,ftem1,vl)) then
- <<if print_ then
- <<terpri()$write "substitution : "$
- substprint(cdr q,car q)>>$
- if member(cdr q,ftem) then
- <<l:=nil$
- while sorg do
- <<a:=car sorg$
- sorg:=cdr sorg$
- if freeof(a,cdr q) then l:=cons(a,l)
- else
- <<ss:=cons(reval subst(car q,cdr q,a),ss)$
- if contradiction(car ss,ftem) then sorg:=nil
- >>
- >>$
- if not contradiction_ then
- <<sorg:=desort reverse l$
- ss:=desort delete(0,union(ss,nil))$
- forg:=fctsubst(car q,cdr q,forg)$
- stop3:=t$
- stop2:=t>>
- >>
- else
- <<stem:=powdelete(stem)$
- stem:=substandsep(q,stem,ftem1,vl)$
- if not contradiction_ then
- <<stem:=powappend(stem,f)$
- s13:=substandsep(q,s13,ftem1,vl)>>$
- stop3:=t$
- if contradiction_ then
- <<stop2:=t$
- fl:=t>>
- >>$
- ftem1:=delete(cdr q,ftem1)>>$
- if not stop3 then
- <<s:=cons(s,ldiff(s,f))$ % fuehr. Abl. + Potenz
- l1:=reverse l1$ % Reihenfolge umkehren
- l:=list(car l1)$ % einfach. Ausgangsgl.
- l1:=cdr l1$ % Glgn., von denen eine
- % ersetzt wird
- while not null l1 and minausdsp(caar l1,s) do
- <<l:=cons(car l1,l)$ % erste einfachere Gl.
- l1:=cdr l1>>$
- l1:=cons(list(s,not cadar l),l)$ % durch s ersetzen
- l:=nil$
- stop3:=nil$
- if print_ and (length l1=2) then
- <<eqprint(list('EQUAL,0,car s))$
- write "with leading derivative "$
- write reval list('EXPT,mkldiff(f,cadr s),cddr s)$
- write " replaces a de with "$
- write reval list('EXPT,mkldiff(f,cadr p),cddr p)$
- terpri()>> >> >>
- until ((length(l1)=2) and not q) or stop3$
- if not stop3 then
- if null car s or zerop car s then
- <<stem:=delete(p,stem)$
- p:=nil$
- stop2:=nil>>
- else
- <<
- stem:=cons(s,delete(p,stem))$
- % stem:=cons(cons(algsimple(car s,f,ftem1,vl),cdr s),
- % delete(p,stem))$
- p:=s:=nil$
- stop2:=nil>> >>
- else
- <<s13:=cons(car p,s13)$
- stem:=delete(p,stem)$
- p:=nil$
- stop3:=t>>$
- if not stop3 then goto rep1$
- end$
- stem:=powdelete(stem)>>$
- if not stop2 then goto rep$
- return if fl then list(sorg,ss,forg,fges,decfl)
- else nil
- end$
- endmodule$
- %*********************************************************************
- module separation$
- %*********************************************************************
- % Routines for separation of de's
- % Author: Andreas Brand
- % June 1990
- symbolic procedure termsep(a,x)$
- % Separieren eines Produktes
- % a Ausdr. in LISP-Notation, x Var.
- % Ergebnis: nil, falls keine Sep. moegl.,sonst Liste von Produktpaaren
- % ((abh. von x)(unabh. von x))
- begin scalar l,p,q,sig,l1,l2$
- if atom a then if x=a then l:=list(a,1) % Variable
- else l:=list(1,a) % Konstante
- else if freeof(a,x) then l:=list(1,a) % a unabh. von x
- else
- <<if car a='MINUS then
- <<a:=cadr a$
- sig:=not sig>>$
- if pairp a and (car a='TIMES) then l:=cdr a
- else l:=list a$ % l Liste der Faktoren
- p:=nil$ % Liste der Faktoren,
- % die x enthalten
- q:=nil$ % Liste der Faktoren,
- % die x nicht enth.
- while l do
- <<if freeof(car l,x) then q:=cons(car l,q)
- % Faktor enth. x nicht
- else
- if pairp car l and (caar l='EXPT) and freeof(cadar l,x)
- and (pairp caddar l) and (car caddar l='PLUS) then
- <<for each s in cdr caddar l do
- if freeof(s,x) then l1:=cons(s,l1)
- else l2:=cons(s,l2)$
- if l1 then
- <<if cdr l1 then l1:=cons('PLUS,l1)
- else l1:=car l1$
- q:=cons(list('EXPT,cadar l,l1),q)>>$
- if l2 and cdr l2 then l2:=cons('PLUS,l2)
- else l2:=car l2$
- p:=cons(list('EXPT,cadar l,l2),p)>>
- else p:=cons(car l,p)$
- l:=cdr l>>$
- if p then
- if length p>1 then p:=cons('TIMES,p)
- else p:=car p$
- if q then
- if length q>1 then q:=cons('TIMES,q)
- else q:=car q$
- if p or q then % war Sep. moegl.?
- if null p then if sig then l:=list(1,list('MINUS,q))
- else l:=list(1,q)
- else
- if q then if sig then l:=list(p,list('MINUS,q))
- else l:=list(p,q)
- else if sig then l:=list(p,list('MINUS,1))
- else l:=list(p,1)>>$
- return l
- end$
- symbolic procedure sumsep(l,x)$
- % Separieren einer Liste von Summanden
- % l Liste von Ausdr. in LISP-Notation, x Var.
- % Ergebnis: nil, falls keine Sep. moegl.,
- % sonst Liste von Listen von Summanden
- begin scalar cl,p,q,s,qflag$
- while l do
- if pairp car l and (caar l='QUOTIENT) then
- <<p:=termsep(cadar l,x)$
- if not q then q:=termsep(caddar l,x)$ % Quotient immer gleich
- if p and q then
- <<l:=cdr l$
- if car q=1 then s:=car p
- else s:=list('QUOTIENT,car p,car q)$
- if cadr q=1 then p:=list(s,cadr p)
- else p:=list(s,list('QUOTIENT,cadr p,cadr q))$
- cl:=termsort(cl,p)>>
- else
- <<l:=nil$
- cl:=nil>> >>
- else
- <<p:=termsep(car l,x)$ % Separieren eines Summanden
- if p then % erfolgreich
- <<l:=cdr l$
- cl:=termsort(cl,p)>> % Terme einordnen
- else
- <<l:=nil$ % Abbruch
- cl:=nil>> >>$ % keine Separ. moegl.
- if cl and print_ and (length cl>1) then % output for test
- <<terpri()$write "separation w.r.t. "$fctprint list x$
- % of linear independence
- if pairp caar cl and (caaar cl='QUOTIENT) then
- l:=for each s in cl collect cadar s
- else l:=for each s in cl collect car s$
- if not linearindeptest(l,list x) then cl:=nil>>$
- return cl % Liste der Terme, die von x
- % unabh. sind
- end$
- symbolic procedure linearindeptest(l,vl)$
- begin scalar l1,flag$
- l1:=l$flag:=t$
- while flag and pairp l1 do
- if atom car l1 then l1:=cdr l1
- else if member(car l1,vl) then l1:=cdr l1
- else if caar l1='EXPT then
- if (numberp caddar l1) and member(cadar l1,vl)
- then l1:=cdr l1
- else flag:=nil
- else flag:=nil$
- if not flag then
- <<terpri()$write "linear independent expressions : "$
- for each x in l do eqprint(x)$
- if independence_ then
- if yesp "Are the expressions linear independent?"
- then flag:=t
- else flag:=nil
- else flag:=t>>$
- return flag$
- end$
- symbolic procedure termsort(cl,p)$
- % Zusammenfassen der Terme
- % cl Liste von Paaren,p Paar
- % Sind bei einem Element von cl und bei p die ersten Teile gleich,
- % wird der zweite Teil von p an den des Elements von cl angehaengt
- if null cl then list p
- else if caar cl=car p then cons(cons(car p,cons(cadr p,cdar cl)),cdr cl)
- else cons(car cl,termsort(cdr cl,p))$
- symbolic procedure eqsep(eql,ftem)$
- % Separieren einer Liste von: Gl. als Liste von Koeffizient + Summ.
- % + Liste der Var. nach denen schon erfolglos sep. wurde
- % + Liste der Var. nach denen noch nicht separiert wurde
- begin scalar vlist1,vlist2,a,x,l,eql1$
- while eql do
- <<a:=car eql$ % erste Gl. +Listen
- vlist1:=cadr a$ % Var. nach d. erfolglos sep. wurde
- vlist2:=caddr a$ % Var. nach denen nicht sep. wurde
- eql:=cdr eql$
- if vlist2 then % noch Var. zu sep.
- <<x:=car vlist2$
- vlist2:=cdr vlist2$
- if freeof(cdar a,x) then
- eql:=cons(list(car a,vlist1,vlist2),eql)
- else if fctdepend(list cdar a,x,ftem)
- then eql:=cons(list(car a,cons(x,vlist1),vlist2),eql)
- else <<l:=sumsep(cdar a,x)$
- % Liste der Gl. die sich durch Sep.
- % nach x ergeben
- if l then
- eql:=append(varapp(l,caar a,nil,append(vlist1,vlist2)),eql)
- % nach erfolgr. Sep. wird nach bisher
- % erfolglosen Var. sep.
- else
- eql:=cons(list(car a,cons(x,vlist1),vlist2),eql)>> >>
- % erfolgloses Sep.,x wird als
- % erfolglose Var. registriert
- else eql1:=cons(a,eql1)>>$
- return eql1$
- end$
- symbolic procedure separ(p,ftem,varl)$
- % Die Gl. p (in LISP-Notation) wird nach den Var. aus varl separiert
- % varl Liste der Variabl.
- begin scalar eql,eqlist,a,q,l,s$
- if p and not zerop p then
- if not (pairp p and (car p='QUOTIENT) and
- intersection(argset smemberl(ftem,cadr p),varl)) then
- <<if pairp p and (car p='QUOTIENT) then
- <<q:=cdr reval list('FACTORIZE,caddr p)$
- if length q>1 then q:=cons('TIMES,q)
- else q:=car q$
- p:=cadr p>>$
- if pairp p and (car p='PLUS) then
- a:=cons(nil,if not q then cdr p
- else for each b in cdr p collect list('QUOTIENT,b,q))
- else if not q then a:=list(nil,p)
- else a:=list(nil,List('QUOTIENT,p,q))$
- % Gl. als Liste von Summanden
- eql:=list(list(a,nil,varl))$
- % Listen der Var. anhaengen
- eql:=eqsep(eql,ftem)$
- while eql do
- <<a:=caar eql$ % Listen der Var. streichen
- if cddr a then a:=cons(car a,cons('PLUS,cdr a))
- else a:=cons(car a,cadr a)$ % PLUS eintragen
- if car a then
- if cdar a then a:=cons(cons('TIMES, car a),cdr a)
- else a:=cons(caar a,cdr a)
- else a:=cons(1,cdr a)$
- eqlist:=cons(a,eqlist)$
- eql:=cdr eql>>$
- if eqlist then eql:=union(cdr eqlist,list car eqlist)$
- eqlist:=nil$
- while eql do
- <<a:=car eql$
- l:=eql:=cdr eql$
- for each b in l do
- <<s:=reval list('QUOTIENT,cdr b,cdr a)$
- if not smemberl(append(varl,ftem),s) then
- <<eql:=delete(b,eql)$
- a:=cons(reval list('PLUS,car a,list('TIMES,s,car b)),cdr a)>>
- >>$
- eqlist:=cons(a,eqlist)>>
- >>
- else eqlist:=list cons(1,p) % FTEM functions in the DENR of p
- else eqlist:=list cons(0,0)$
- return eqlist
- end$
- symbolic procedure separate(p,ftem,vl)$
- % Separieren der Gleichung p
- % ftem Liste der Fkt., vl Liste der Variablen
- begin scalar l$
- l:=separ(p,ftem,vl)$
- return union(for each a in l collect cdr a,nil)$
- end$
- symbolic procedure stardep(p,ftem,vl)$
- % yields: nil, if a function (from ftem) in p depends
- % on all variables (from vl)
- % cons(v,n) otherweise, with v beeing the list of variables
- % on which depend a minimal number n of functions
- begin scalar ft,b,v,n,m$
- vl:=varslist(p,ftem,vl)$
- ft:=ftem$
- n:=length vl$ % number of all variables
- while not b and ft do % searching a fct of all vars
- if fctlength car ft=n then b:=t
- else ft:=cdr ft$
- if not b then
- <<n:=sub1 length ftem$
- while vl do % searching var.s on which depend a
- % minimal number of functions
- <<if n> (b:=for each h in ftem sum
- if member(car vl,fctargs h) then 1
- else 0)
- then
- <<n:=b$ % a new minimum
- v:=list car vl>>
- else if b=n then v:=cons(car vl,v)$
- vl:=cdr vl>> >>$
- return if v then cons(v,n) % on each varible from v depend n
- % functions
- else nil
- end$
- symbolic procedure gensep(p,ftem,vl)$
- begin scalar a,aa,q,l,l1,l2,ft,pl,pl1,x,vl1,deno$
- if pairp p and (car p = 'QUOTIENT) then
- <<casecheck(caddr p,vl)$
- p:=cadr p>>$
- ftem:=smemberl(ftem,p)$
- vl:=varslist(p,ftem,vl)$
- if not (a:=stardep(p,ftem,vl)) then pl:=list p
- else
- if zerop cdr a then pl:=separate(p,ftem,vl)
- else % e.g. a=((x y z).2)
- <<pl:=list p$
- if print_ then <<terpri()$write "generalized separation ">>$
- if tr_gensep then
- <<terpri()$write "de to be separated : "$
- deprint pl$
- terpri()$write "variable list for separation : ",car a$
- terpri()$write "for each of these variables ",cdr a,
- " function(s) do(es) not depend on it ">>$
- a:=car a$
- while pairp a do % i.e. being able to make a separation for each x
- <<x:=car a$
- q:=p$ ft:=l2:=nil$
- for each f in ftem do
- if member(x,fctargs f) and not freeof(q,f) then ft:=cons(f,ft)$
- if tr_gensep then
- <<terpri()$write "to separate directly w.r.t. ",x$
- write " the expression : "$deprint list q$
- write "will be differentiated">>$
- ft:=reverse fctsort ft$ % sorting w.r.t. number of args
- while ft do % throwing out all functions ft
- if (l1:=felim(q,car ft,ftem,vl,l2)) then
- <<q:=car l1$
- l2:=cadr l1$
- aa:=stardep(q,smemberl(ftem,q),vl)$
- if not aa or zerop cdr aa then ft:=nil
- else ft:=smemberl(cdr ft,q)>>
- else ft:=nil$
- if l1 then
- <<if (pairp q) and (car q='QUOTIENT) then
- <<deno:=caddr q;
- q:=cadr q;
- ft:=smemberl(ftem,q)>> else ft:=nil$
- if ft then
- <<vl1:=nil$
- for each y in vl do if freeof(ft,y) then vl1:=cons(y,vl1)
- >> else vl1:=vl$
- if aa and zerop cdr aa and not (q=0) then
- <<if tr_gensep then
- <<terpri()$write "trying direct separation of "$
- deprint list q$
- write "Remaining variables: ",vl1>>;
- l1:=separate(q,ftem,vl1)$
- if tr_gensep then
- <<write "The result of direct separation: "$deprint l1>>$
- if l1 and cdr l1 and tr_gensep then
- <<terpri()$
- write "direct separation yields ",length l1," equations">>
- >>
- else l1:=list q$
- l1:=desort for each s in l1 collect
- if deno then
- backint(list('QUOTIENT,s,deno),l2,union(fnew_,ftem),vl)
- else
- backint(s,l2,union(fnew_,ftem),vl)$
- ftem:=union(fnew_,ftem)$
- pl:=union(l1,pl)$
- while l1 do
- <<if fctevalstrict(car l1,ftem,vl) then
- <<pl1:=cons(car l1,pl1)$
- a:=list nil$
- ftem:=union(fnew_,ftem)>>$
- l1:=cdr l1
- >>
- >>$
- a:=cdr a$
- >>$
- >>$
- if pl1 then
- pl:=union(simplifyde(reval p,union(fnew_,ftem),vl,genint_),pl1)
- else
- <<if null pl then pl:=list p$
- for each p in pl do
- pl1:=union(simplifyde(reval p,union(fnew_,ftem),vl,genint_),pl1)$
- pl:=pl1$
- if cdr pl and print_ then
- <<terpri()$write "separation yields ",length pl," equations"$
- if tr_gensep then deprint pl>>
- >>$
- return pl$
- end$
- symbolic procedure felim(q,f,ftem,vl,l2)$
- begin scalar a,b,l,l1,ft1,v,prflag$
- v:=car setdiff(vl,fctargs f)$ % getting rid of f through diff. wrt. v
- ft1:=nil$
- for each f in ftem do if freeof(f,v) then ft1:=cons(f,ft1)$
- if not (pairp q and (car q='QUOTIENT) and smemberl(ft1,caddr q)) then
- <<prflag:=print_$print_:=nil$
- l:=desort separ(q,ft1,list v)$ % det. all lin. ind. factors with v
- if tr_gensep then
- <<terpri()$write "To get rid of ",f,
- " we differentiate w.r.t. : ",v>>$
- print_:=prflag$
- l1:=nil$
- while l do
- <<a:=car l$
- b:=nil$
- if not freeof(cdr a,f) and (not zerop car a) then
- if (pairp cdr a) and (cadr a='PLUS) then
- <<for each c in cddr a do if not freeof(c,f) then b:=cons(c,b)$
- if b then b:=cons('PLUS,b)>> else b:=cdr a$
- if b then
- <<l1:=cons(car a,l1)$
- % q:=reval list('DIFFERENCE,q,list('times,b,car a))
- >>$
- l:=cdr l
- >>$
- if tr_gensep then
- <<terpri()$
- write v," - depending coefficients of terms containing ",f," : "$
- for each ss in l1 do eqprint ss>>$
- while l1 do
- <<b:=reval aeval car l1$
- l1:=cdr l1$
- if not zerop b then
- <<a:=reval aeval list('QUOTIENT,q,b)$
- % if not freeof(a,v) then
- <<l:=cons(b,l)$ % for later backward integrations
- l1:=for each c in l1 collect
- reval list('DF,list('QUOTIENT,c,b),v)$
- casecheck(b,vl)$
- q:=reval list('DF,a,v)$
- if tr_gensep then
- <<write "The new equation: "$
- eqprint q$
- write "The remaining factors:"$
- for each ss in l1 do eqprint ss$>>
- >>
- % else <<q:=b$ l1:=nil>>
- >>
- >>$
- if tr_gensep then
- <<terpri()$write "new expression (should not depend on ",f,") : "$
- eqprint q$>>$
- if l then l2:=cons(list(v,l),l2)$
- if tr_gensep and l2 then
- <<write "The list of multiplications and integrations ",
- "to go backwards after direct separation:"$
- for each ss in l2 do <<write "integr. wrt. ",car ss$ terpri()$
- write "multiply with "$
- for each aa in cadr ss do
- eqprint aa>>
- >>$
- l1:=list(q,l2)
- >>$
- return l1$
- end$
- symbolic procedure backint(s,l2,ftem,vl)$
- begin scalar fl,ft,q,l,v,v1,v2,vf,s2,p,f2,fnew1$
- fnew1:=fnew_$
- fl:=q:=t$
- p:=s$
- while l2 and fl do
- <<l:=car l2$
- l2:=cdr l2$
- v:=car l$
- if tr_gensep then
- <<terpri()$
- write "backward integration w.r.t. ",v," of the expression : "$
- eqprint p>>$
- l:=cadr l$
- while l and q do
- <<ft:=smemberl(ftem,p)$
- %terpri()$write "vor int: p= "$eqprint p$
- fnew_:=nil$
- q:=integratepde(p,ft,list v,nil)$
- fnew1:=append(fnew_,fnew1)$
- if q then
- <<fl:=t$
- p:=reval list('TIMES,car l,car q)$
- % Substituting the new functions of integration by derivatives
- % of them such that back-integration can be made
- % hat fnew_ nur ein element, d.h. wird nur eine Integration gemacht
- % oder mehrere?
- for each f1 in fnew_ do
- <<f2:=f1;
- vf:=setdiff(vl,fctargs f1);
- for each s1 in reverse l2 do
- <<v1:=car s1;
- if not freeof(f1,v1) then
- % only then integration makes difficulties
- <<s2:=reverse cadr s1;
- while s2 do
- <<if not smemberl(vf,car s2) then
- f2:=reval list('DF,list('QUOTIENT,f2,car s2),v1);
- % actually only dividing through those factors of (car s2)
- % which depend on v1 and which do not contain variables
- % which f2 does not depent on.
- s2:=cdr s2
- >>
- >>
- >>;
- if f1 neq f2 then
- <<if tr_gensep then <<write f1," is replaced by ";
- eqprint f2>>;
- p:=subst(f2,f1,p)$
- >>
- >>;
- ftem:=union(fnew_,ftem)
- >> else fl:=nil$
- l:=cdr l
- >>
- >>$
- if tr_gensep then if fl then <<terpri()$write "yields : "$
- eqprint p$>>
- else <<terpri()$
- write "was not successful.">>$
- fnew_:=union(fnew1,fnew_)$
- return if fl then p
- else s
- end$
- endmodule$
- %*********************************************************************
- module integration$
- %*********************************************************************
- % Routines for integration of pde's containing unnown functions
- % Author: Andreas Brand
- % June 1990
- symbolic procedure intminderiv(p,ftem,vl)$
- % yields a list l of variables from vl, such that in p for each v in l
- % there exists a derivative of a function f from ftem wrt. v
- begin scalar l,ft,a$
- for each v in vl do
- <<ft:=ftem$
- a:=t$
- while ft and a do
- if member(v,fctargs car ft) and highdiff(p,car ft,v)=0
- then a:=nil
- else ft:=cdr ft$
- if a then l:=cons(v,l)>>$
- return l
- end$
- symbolic procedure integratepde(p,ftem,varl,genflag)$
- % Generalized integration of the expression p
- % if not genflag then "normal integration"
- begin scalar l,l1,vl,n,l2,q,f,flag$
- varl:=if genflag then varl
- else intminderiv(p,ftem,varl)$
- q:=p$
- if tr_genint then
- <<terpri()$write "integration of the expression : "$
- eqprint q$
- >>$
- %ftem:=smemberl(ftem,p)$
- vl:=argset ftem$
- if zerop q then flag:='success$
- for each v in varl do
- if member(v,vl) then
- <<n:=0$
- if flag neq 'genint then
- <<p:=q$
- while (p and not zerop p) do
- <<p:=intpde(q,ftem,vl,v)$
- if p then
- if zerop cadr p then
- <<q:=car p$n:=add1 n$flag:='success>>
- else if not genflag then p:=nil
- else
- <<if print_ and not tr_genint then
- <<terpri()$write "generalized Integration: ">>$
- if (l2:=partint(cadr p,smemberl(ftem,cadr p),vl,v))
- then
- <<q:=reval list('PLUS,car p,car l2)$
- n:=add1 n$
- ftem:=union(fnew_,ftem)$
- l1:=append(cdr l2,l1)$ % additional de's
- flag:='genint>>$
- p:=nil$
- >>$
- >> >>$
- l:=cons(n,l)>>
- else
- <<flag:='success$
- l:=cons(1,l)$
- q:=reval list('INT,q,v)$
- l2:=1$
- >>$
- l:=reverse l$
- if not zerop q then
- <<ftem:=reverse fctsort ftem$
- while ftem do
- if fctlength car ftem<length vl then ftem:=nil
- else if fctlinear(q,f) then
- <<f:=car ftem$
- ftem:=nil>>
- else ftem:=cdr ftem$
- if f then
- <<l2:=lderiv2(q,f,fctargs f)$
- l2:=reval coeffn(q,reval car l2,cdr l2)>>
- else l2:=1>>$
- return if flag then
- <<if varl then
- for each v in varl do
- <<q:=list('PLUS,q,intconst(l2,vl,v,car l))$
- l:=cdr l>>
- else
- <<q:=newfct(fname_,vl,nfct_)$
- nfct_:=add1 nfct_$
- fnew_:=cons(q,fnew_)>>$
- cons(reval q,l1)>>
- else nil$
- end$
- symbolic procedure intpde(p,ftem,vl,x)$
- % integration of an polynomial expr. p w.r.t. x
- begin scalar f,ft,l,l1,vl,k,s,a,iflag,flag$
- ft:=ftem$
- vl:=cons(x,delete(x,vl))$
- while ftem and not flag do
- <<f:=car ftem$
- if member(x,fctargs f) then
- <<l1:=l:=lderiv2(p,f,vl)$
- while not (flag or (iflag:=intlintest(p,f,l,x))) do
- <<k:=reval coeffn(p,car l,cdr l)$
- if intcoefftest(lderiv1(k,f,vl),car l,vl) then
- <<a:=decderiv(car l,x)$
- k:=reval list('INT,subst('v_a_r_,a,k),'v_a_r_)$
- k:=reval subst(a,'v_a_r_,k)$
- s:=cons(k,s)$
- p:=reval aeval list('DIFFERENCE,p,list('DF,k,x))$
- if (l:=lderiv2(p,f,vl))=l1 then flag:='neverending
- else l1:=l>>
- else flag:='coeffld>>$
- if iflag='nofct then ftem:=smemberl(ftem,p)
- else if fctlength f<length vl then
- <<ftem:=cdr ftem$flag:=nil>>
- else flag:=(iflag or flag)>>
- else ftem:=cdr ftem>>$
- return if not flag then
- <<l:=explicitpart(p,ft,x)$
- l:=list('INT,l,x)$
- s:=reval aeval cons('PLUS,cons(l,s))$
- p:=list('DIFFERENCE,p,list('DF,l,x))$
- if poly_only then
- if ratexp(s,ft) then list(s,reval aeval p)
- else nil
- else list(s,reval aeval p)>>
- else nil$
- end$
- symbolic procedure explicitpart(p,ft,x)$
- % part of a sum p which only explicitly depends on a variable x
- begin scalar l$
- if not member(x,argset smemberl(ft,p)) then l:=p
- else if pairp p then
- <<if car p='MINUS then
- l:=reval list('MINUS,explicitpart(cadr p,ft,x))$
- if (car p='QUOTIENT) and not member(x,argset smemberl(ft,caddr p))
- then
- l:=reval list('QUOTIENT,explicitpart(cadr p,ft,x),caddr p)$
- if car p='PLUS then
- <<for each a in cdr p do
- if not member(x,argset smemberl(ft,a)) then l:=cons(a,l)$
- if pairp l then l:=reval cons('PLUS,l)>> >>$
- if not l then l:=0$
- return l$
- end$
- symbolic procedure intconst(co,vl,x,n)$
- begin scalar l,l2,f$
- while n>0 do
- <<f:=newfct(fname_,delete(x,vl),nfct_)$
- nfct_:=add1 nfct_$
- fnew_:=cons(f,fnew_)$
- l:=cons(list('TIMES,f,list('EXPT,x,sub1 n)),l)$
- n:=sub1 n>>$
- if length l>1 then l:=cons('PLUS,l)
- else if pairp l then l:=car l
- else l:=0$
- if length co>1 then
- <<if car co='TIMES then co:=cdr co
- else co:=list co$
- while co do
- <<if freeof(car co,x) then l2:=cons(car co,l2)$
- co:=cdr co>> >>
- else if not (co=x) then l2:=list co$
- return reval if l2 then cons('TIMES,cons(l,l2))
- else l
- end$
- symbolic procedure intcoefftest(l,l1,vl)$
- begin scalar s$
- return if not pairp l then t
- else if car l='DF then
- <<s:=decderiv(l1,car vl)$
- if pairp s and pairp cdr s then s:=cddr s
- else s:=nil$
- if diffrelp(cons(cddr l,1),cons(s,1),vl) then t
- else nil>>
- else t$
- end$
- symbolic procedure fctlinear(p,f)$
- <<p:=ldiffp(p,f)$
- (null car p) and (cdr p=1)>>$
- symbolic procedure intlintest(p,f,l,x)$
- % Test,ob in einem Paar (fuehrende Ableitung.Potenz)
- % eine x-Ableitung linear auftritt
- if freeof(l,f) then 'nofct
- else % Fkt. tritt auf
- if (car l) and (cdr l=1) then % fuehr. Abl. tritt linear auf
- if pairp car l then % Abl. der Fkt. tritt auf
- if caar l='DF then
- if member(x,cddar l) then nil
- % x-Abl. tritt auf
- else if member(x,fctargs cadar l) then 'noxdrv
- else 'noxdep
- else 'nodrv
- else 'nodrv
- else 'nonlin$
- symbolic procedure decderiv(l,x)$
- % in Diff.ausdr. l wird Ordn. d. Abl. nach x um 1 gesenkt
- begin scalar l1$
- return if l then if car l='DF then
- <<l1:=decderiv1(cddr l,x)$
- if l1 then cons('DF,cons(cadr l,l1))
- else cadr l>>
- else l
- else nil$
- end$
- symbolic procedure decderiv1(l,x)$
- if null l then nil
- else
- if car l=x then
- if cdr l then
- if numberp cadr l then
- if cadr l>2 then cons(car l,cons(sub1 cadr l,cddr l))
- else cons(car l,cddr l)
- else cdr l
- else nil
- else cons(car l,decderiv1(cdr l,x))$
- symbolic procedure integratede(q,ftem,genflag)$
- % Integration of a de
- % result: newde if successfull, nil otherwise
- begin scalar l,l1,fl$
- ftem:=smemberl(ftem,q)$
- if l1:=integrableode(q,ftem) then % looking for an integrable ode
- if l1:=integrateode(q,car l1,cadr l1,caddr l1,ftem) then
- % trying to integrate it
- <<q:=l1$
- fl:=t$
- ftem:=smemberl(union(fnew_,ftem),q)>>$
- if l1:=integratepde(q,ftem,argset ftem,genflag) then
- % trying to integrate a pde
- <<q:=car l1$
- l:=for each a in cdr l1 conc
- <<ftem:=union(fnew_,ftem)$
- if (l1:=integratede(a,ftem,nil)) then l1
- else list a>>$
- fl:=t$
- ftem:=union(fnew_,ftem)>>$
- if fl then
- <<l:=cons(q,l)$
- l:=for each a in l collect reval aeval a$
- l:=for each a in l collect
- if pairp a and (car a='QUOTIENT) then cadr a
- else a>>$
- return l$
- end$
- endmodule$
- %********************************************************************
- module generalized_integration$
- %********************************************************************
- % Routines for generalized integration of pde's containing unnown
- % functions
- % Author: Andreas Brand
- % December 1991
- symbolic procedure gintorder(p,vl,x)$
- % reorder equation p
- begin scalar l,l1,q,n,m,b,c$
- if pairp(p) and (car p='QUOTIENT) then <<q:=caddr p$p:=cadr p>>$
- if pairp(p) and (car p='PLUS) then p:=cdr p % list of summands
- else p:=list p$
- while p do
- <<l1:=gintorder1(car p,x)$
- if DepOnAllVars(car l1,x,vl) then l:=p:=nil
- else <<l:=termsort(l,l1)$p:=cdr p>> >>$
- if l then
- <<l:=for each a in l collect if cddr a then
- <<b:=car a$
- c:=cdr reval coeff1(cons('PLUS,cdr a),x,nil)$
- m:=0;
- while c and (car c=0) do <<c:=cdr c$m:=add1 m>>$
- if m>0 then b:=list('TIMES,list('EXPT,x,m),b)$
- cons(reval b,c)>>
- else
- cons(reval car a,cdr reval coeff1(cadr a,x,nil))$
- if q then
- if freeof(q,x) then
- l:=for each a in l collect
- cons(car a,for each s in cdr a collect
- reval list('QUOTIENT,s,q))
- else
- l:=for each a in l collect
- cons(reval list('QUOTIENT,car a,q),cdr a)>>$
- return l$
- end$
- symbolic procedure DepOnAllVars(c,x,vl)$
- % tests for occurence off vars from vl in factors of c depending on x
- begin scalar l$
- if pairp c and (car c='TIMES) then c:=cdr c
- else c:=list c$
- while c and vl do
- <<if not freeof(car c,x) then
- for each v in vl do if not freeof(car c,v) then l:=cons(v,l)$
- vl:=setdiff(vl,l)$
- c:=cdr c
- >>$
- return (null vl)$
- end$
- symbolic procedure gintorder1(p,x)$
- % reorder a term p
- begin scalar l1,l2,sig$ % l2:list of factors of p not depending
- % on x or beeing a power of x
- % l1:all other factors
- if pairp p and (car p='MINUS) then <<sig:=t$p:=cadr p>>$
- if pairp p and (car p='TIMES) then p:=cdr p
- else p:=list p$
- for each a in p do
- <<if freeof(a,x) then l2:=cons(a,l2)
- else if a=x then l2:=cons(a,l2)
- else if pairp a and (car a='EXPT) and (cadr a=x) and fixp caddr a
- then l2:=cons(a,l2)
- else l1:=cons(a,l1)>>$
- if pairp l1 then
- if cdr l1 then l1:=cons('TIMES,l1)
- else l1:=car l1$
- if pairp l2 then
- <<if cdr l2 then l2:=cons('TIMES,l2)
- else l2:=car l2$
- if sig then l2:=list('MINUS,l2)>>$
- return list(if l1 then l1 else 1,if l2 then l2 else 1)$
- end$
- symbolic procedure partint(p,ftem,vl,x)$
- begin scalar f,neg,l1,n,k,l$
- if tr_genint then
- <<terpri()$write "generalized integration of the unintegrated rest : "$
- eqprint p
- >>$
- l:=gintorder(p,vl,x)$
- l:=for each s in l collect
- <<f:=newfct(fname_,varslist(car s,ftem,vl),nfct_)$
- nfct_:=add1 nfct_$
- fnew_:=cons(f,fnew_)$
- neg:=t$
- n:=sub1 length cdr s$
- k:=-1$
- l1:=cons(reval list('DIFFERENCE,car s,list('DF,f,x,add1 n)),l1)$
- reval cons('PLUS, for each sk on cdr s collect
- <<neg:=not neg$
- k:=add1 k$
- reval
- list('TIMES,if neg then -1 else 1,
- list('DF,f,x,n-k),tailsum(sk,k,x))>>
- )>>$
- if l then l:=cons(reval cons('PLUS,l),l1)$
- if tr_gensep then
- <<terpri()$
- write "result (without constant or function of integration): "$
- eqprint car l$
- write "additional equations : "$
- deprint cdr l$
- >>$
- return l$
- end$
- symbolic procedure tailsum(sk,k,x)$
- begin scalar j$
- j:=-1$
- return reval cons('PLUS,
- for each a in sk collect
- <<j:=j+1$
- list('TIMES,a,prod(j+1,j+k),list('EXPT,x,j)) >> )
- end$
- symbolic procedure prod(m,n)$
- if m>n then 1
- else for i:=m:n product i$
- endmodule$
- %********************************************************************
- module odeintegration$
- %********************************************************************
- % Routines for integration of ode's containing unnown functions
- % Author: Thomas Wolf
- % August 1991
- symbolic procedure integrateode(de,fold,xnew,ordr,ftem)$
- begin scalar newde,newnewde,l$ % Integral einer Gleichung
- return
- if not xnew then begin % Integr. einer alg. Gl. fuer eine Abl.
- newde:=cadr solveeval list(de,fold)$
- newde:=reval list('PLUS,cadr newde,list('MINUS,caddr newde))$
- return
- if newde=de then nil
- else if (l:=integratede(newde,ftem,nil)) then car l
- else nil$
- end else % eine ode fuer ein f?
- if not pairp fold then % i.e. not df(...,...), i.e. fold=f
- odeconvert(de,fold,xnew,ordr)
- % --> ode fuer eine Abl. von f
- else begin
- newde:=odeconvert(de,fold,xnew,ordr)$
- % if print_ then
- % if newde then <<terpri()$write "One ode has been integrated.">>$
- return
- if not newde then nil
- else <<
- newnewde:=cadr solveeval list(newde,fold)$
- newnewde:=reval list('PLUS,cadr newnewde,list('MINUS,
- caddr newnewde))$
- ftem:=union(fnew_,ftem)$
- newnewde:=integratede(newnewde,ftem,nil)$
- if newnewde then car newnewde
- else newde >>
- end
- end$ % of integrateode
- symbolic operator polyans$
- symbolic procedure polyans(ordr,dgr,x,y,d_y,fn)$
- % generates a polynom
- % for i:=0:dgr sum fn"i"(x,y,d_y(1),..,d_y(ordr-1))*d_y(ordr)**i
- % with fn as the function names and d_y as names or derivatives of y
- % w.r.t. x
- begin scalar ll,fl,a,i,f$
- i:=sub1 ordr$
- while i>0 do
- <<ll:=cons(list(d_y,i),ll)$
- i:=sub1 i>>$
- ll:=cons(y,ll)$
- ll:=reverse cons(x,ll)$
- fl:=nil$
- i:=0$
- while i<=dgr do
- <<f:=newfct(fn,ll,i)$
- fl:=(f . fl)$
- a:=list('PLUS,list('TIMES,f,list('EXPT,list(d_y,ordr),i)),a)$
- i:=add1 i>>$
- return list('list,reval a,cons('list,fl))
- end$ % of polyans
- symbolic operator sepans$
- symbolic procedure sepans(kind,v1,v2,fn)$
- % Generates a separation ansatz
- % v1,v2 = lists of variables, fn = new function name + index added
- % The first variable of v1 occurs only in one sort of the two sorts of
- % functions and the remaining variables of v1 in the other sort of
- % functios.
- % The variables of v2 occur in all functions.
- % Returned is a sum of products of each one function of both sorts.
- % form: fn1(v11;v21,v22,v23,..)*fn2(v12,..,v1n;v21,v22,v23,..)+...
- % the higher "kind", the more general and difficult the ansatz is
- % kind = 0 is the full case
- begin scalar n,vl1,vl2,h1,h2,h3,h4,fl$
- if cdr v1 = nil then <<vl1:=cdr v2;vl2:=cdr v2>>
- else <<vl1:=cons(cadr v1,cdr v2)$
- vl2:=append(cddr v1,cdr v2)>>$
- return
- if kind = 0 then <<vl1:=append(cdr v1,cdr v2)$
- h1:=newfct(fn,vl1,'_)$
- list('LIST,h1,list('LIST,h1))>>
- else
- if kind = 1 then <<h1:=newfct(fn,vl1,1)$
- list('LIST,h1,list('LIST,h1))>>
- else
- if kind = 2 then <<h1:=newfct(fn,vl2,1)$
- list('LIST,h1,list('LIST,h1))>>
- else
- if kind = 3 then <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- list('LIST,reval list('PLUS,h1,h2),
- list('LIST,h1,h2))>>
- else
- if kind = 4 then <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- list('LIST,reval list('TIMES,h1,h2),
- list('LIST,h1,h2))>>
- else
- if kind = 5 then <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- h3:=newfct(fn,vl1,3)$
- list('LIST,reval list('PLUS,list('TIMES,h1,h2),h3),
- list('LIST,h1,h2,h3))>>
- else
- if kind = 6 then <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- h3:=newfct(fn,vl2,3)$
- list('LIST,reval list('PLUS,list('TIMES,h1,h2),h3),
- list('LIST,h1,h2,h3))>>
- else
- if kind = 7 then <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- h3:=newfct(fn,vl1,3)$
- h4:=newfct(fn,vl2,4)$
- list('LIST,reval list('PLUS,
- list('TIMES,h1,h2),h3,h4),
- list('LIST,h1,h2,h3,h4))>>
- else
- % ansatz of the form FN = FN1(v11,v2) + FN2(v12,v2) + ... + FNi(v1i,v2)
- if kind = 8 then <<n:=1$ vl1:=cdr v1$ vl2:=cdr v2$
- fl:=();
- while vl1 neq () do <<
- h1:=newfct(fn,cons(car vl1,vl2),n)$
- vl1:=cdr vl1;
- fl:=cons(h1, fl);
- n:=n+1
- >>;
- list('LIST, cons('PLUS,fl), cons('LIST,fl))>>
- else
- <<h1:=newfct(fn,vl1,1)$
- h2:=newfct(fn,vl2,2)$
- h3:=newfct(fn,vl1,3)$
- h4:=newfct(fn,vl2,4)$
- list('LIST,reval list('PLUS,list('TIMES,h1,h2),
- list('TIMES,h3,h4)),
- list('LIST,h1,h2,h3,h4))>>
- end$ % of sepans
- symbolic procedure odecheck(ex,fint,ftem)$
- % assumes an revaled expression ex
- % Does wrong if car ex is a list!
- begin scalar a,b,op,ex1$
- %***** ex is a ftem-function *****
- if ex=fint then % list(ex,0,0,..)
- <<a:=list ex$
- ex:=fctargs ex$
- while ex do
- <<a:=append(list(0,0),a)$
- ex:=cdr ex>>$
- % not checked if it is a function of an expression of x
- op:=reverse a>>
- else if pairp ex then
- %***** car ex is 'df *****
- if (car ex)='df then
- <<a:=odecheck(cadr ex,fint,ftem)$
- if not pairp a then op:=a
- else % a is list(fctname,0,0,..,0,0)
- <<op:=list(car a)$
- a:=fctargs car a$ % a is list(variables), not checked
- ex:=cddr ex$ % ex is list(derivatives)
- while a do
- <<ex1:=ex$
- while ex1 and ((car a) neq (car ex1)) do ex1:=cdr ex1$
- if null ex1 then op:=cons(0,cons(0,op))
- else
- <<if not cdr ex1 then b:=1 % b is number of deriv. of that var.
- else
- <<b:=cadr ex1$
- if not numberp b then b:=1>>$
- op:=cons(b,cons(b,op))>>$
- a:=cdr a>>$
- op:=reverse op>> >>
- else
- %***** car ex is a standard or other function *****
- <<a:=car ex$ % for linearity check
- ex:=cdr ex$
- if a='INT then ex:=reval car ex$
- while (op neq '!_abb) and ex do
- <<b:=odecheck(car ex,fint,ftem)$
- if b then % function found
- if b eq '!_abb then op:='!_abb else % occures properly
- if not poly_only then op:=odetest(op,b) else % linearity check
- if (a = 'PLUS) or (a = 'MINUS) then op:=odetest(op,b) else
- if a = 'TIMES then
- if smemberl(ftem,cdr ex) then op:='!_abb
- %car b is f.-name->freeof works
- else op:=odetest(op,b)
- else op:='!_abb$
- ex:=cdr ex>> >>$
- return op
- end$
- symbolic procedure integrableode(p,ftem)$
- if delength p>(if odesolve_ then odesolve_ else 0) then
- (if cont_ then
- if yesp("expression to be integrated ? ") then
- integrableode1(p,ftem))
- else integrableode1(p,ftem)$
- symbolic procedure integrableode1(p,ftem)$
- begin scalar a,b,u,vl,le,uvar,
- fint,fivar,% the function to be integrated and its variables
- fold, % the new function of the ode
- xnew, % the independ. variable of the ode
- ordr1, % order of the ode
- ordr2, % order of the derivative for which it is an ode
- intlist$ % list of ode's
- ftem:=smemberl(ftem,p)$
- vl:=argset ftem$
- % p muss genau eine Funktion aus ftem von allen Variablen enthalten.
- % Die Integrationsvariable darf nicht Argument anderer in p enthaltener
- % ftem-Funktionen sein.
- a:=ftem$
- b:=nil$
- le:=length vl$
- while a and vl do
- <<u:=car a$
- uvar:=fctargs u$
- if (length uvar) = le then
- if b then
- <<vl:=nil$a:=list(nil)>>
- else
- <<b:=t$
- fint:=u$
- fivar:=uvar>>
- else vl:=setdiff(vl,uvar)$
- a:=cdr a>>$
- if not b then vl:=nil$
- le:=length p$
- if ((1<le) and vl) then
- <<a:=odecheck(p,fint,ftem)$
- if not atom a then % The equation is an ode
- <<ordr1:=0$
- ordr2:=0$
- xnew:=nil$
- a:=cdr a$
- b:=fivar$
- while b do
- <<if (car a) neq 0 then
- <<fold:=cons(car b,fold)$
- if (car a) > 1 then fold:=cons(car a,fold)>>$
- ordr2:=ordr2+car a$
- if (car a) neq (cadr a) then
- <<xnew:=car b$
- if not member(xnew,vl) then
- <<b:=list(nil)$vl:=nil>>$
- ordr1:=(cadr a) - (car a)>>$
- b:=cdr b$
- a:=cddr a>>$
- fold:=reverse fold$
- %fold is the list of diff.variables + number of diff.
- if fold then fold:=cons('df,cons(fint,fold))
- else fold:=fint$
- if vl and ((ordr1 neq 0) or (ordr2 neq 0)) then
- intlist:=list(fold,xnew,ordr1,ordr2)
- >> % of variable found
- >>$ % of if
- return intlist
- end$ % of integrableode1
- symbolic procedure odetest(op,b)$
- if not op then b
- else % op=nil --> first function found
- if (car op) neq (car b) then '!_abb else % f occurs in differ. fct.s
- begin scalar dif,a$
- dif:=nil$ % dif=t --> different derivatives
- a:=list(car op)$ % in one variable already found
- op:=cdr op$
- b:=cdr b$
- while op do
- <<a:=cons(max(cadr op,cadr b),cons(min(car op,car b),a))$
- if (car a) neq ( cadr a) then
- if dif then
- <<a:='!_abb$
- op:=list(1,1)>>
- else dif:=t$
- op:=cddr op$
- b:=cddr b>>$
- if not atom a then a:=reverse a$
- return a % i.e. '!_abb or (fctname,min x1-der.,max x1-der.,...)
- end$
- symbolic procedure odeconvert(de,ford,xnew,ordr)$
- begin scalar j,ford_,newco,oldde,newde,newvl,null_$
- ford_:='y_$%gensym$
- depl!*:=delete(assoc(ford_,depl!*),depl!*)$
- depend1(ford_,xnew,t)$
- oldde:=reval subst(ford_,ford,de)$
- if pairp ford then % i.e. (car ford) eq 'DF then
- << for j:=1:ordr do
- oldde:= subst( reval list('DF,ford_,xnew,j),
- reval list('DF,ford,xnew,j), oldde)>>$
- algebraic !!arbconst:=0$
- algebraic ( explog:= {
- cot(~x) => 1/tan(x),
- e**(~x+~y) => e**x*e**y,
- sqrt(e)**(~x+~y) => sqrt(e)**x*sqrt(e)**y,
- e**((~x+~y)/~z) => e**(x/z)*e**(y/z),
- sqrt(e)**((~x+~y)/~z) => sqrt(e)**(x/z)*sqrt(e)**(y/z),
- sqrt(e)**(log(~y)/~x) => y**(1/x/2),
- sqrt(e)**(-log(~y)/~x) => y**(-1/x/2),
- sqrt(e)**(~x*log(~y)/~z) => y**(x/z/2),
- sqrt(e)**(-~x*log(~y)/~z) => y**(-x/z/2),
- sqrt(e)**((~x*log(~y))/~z) => y**(x/z/2),
- e**(log(~y)/~x) => y**(1/x),
- e**(~x*log(~y)/~z) => y**(x/z),
- e**((~x*log(~y))/~z) => y**(x/z),
- int(df(~y,~x)/~y,~x) => log(y) })$
- algebraic( let explog )$
- newde:=algebraic first
- odesolve(symbolic oldde,symbolic ford_,symbolic xnew)$
- if (cadr newde neq oldde) then begin % solution found
- % Test der Loesung klappt nur, wenn Loesung explizit gegeben
- if cadr newde neq ford_ then
- (if print_ then
- <<write "Solution of the ode is not explicitly given."$
- algebraic write "Equation is: ",algebraic symbolic oldde$
- algebraic write "Solution is: ",algebraic symbolic newde>>)
- else begin
- null_:=reval reval aeval subst(caddr newde, ford_, oldde)$
- % reval reval because of a REDUCE bug for special data,
- % to be droped as soon as possible
- if (null_ neq 0) then begin
- % newde:=nil;
- if print_ then begin
- write "odesolve solves : "$
- deprint list oldde$
- write "to"$
- deprint list newde$
- Write "which inserted in the equation yields"$
- deprint list null_$
- end
- end
- end
- end$
- algebraic ( clear explog )$
- depl!*:=delete(assoc(ford_,depl!*),depl!*)$
- if newde then
- <<newde:=list('PLUS,cadr newde,list('MINUS,caddr newde))$
- if 0 = reval list('PLUS,newde,list('MINUS,oldde)) then newde:=nil>>$
- return
- if not newde then nil
- else
- <<newde:=subst(ford,ford_,newde)$
- newvl:=delete(xnew,if (pairp ford and (car ford='DF))
- then fctargs cadr ford
- else fctargs ford)$
- % if pairp ford then newvl:=delete(xnew,cdr assoc(cadr ford,depl!*))
- % else newvl:=delete(xnew,cdr assoc(ford,depl!*))$
- for j:=1:ordr do begin
- newco:=newfct(fname_,newvl,nfct_)$
- nfct_:=add1 nfct_$
- fnew_:=cons(newco,fnew_)$
- newde:=subst(newco,list('arbconst,j),newde)
- % newde:=subst(newco, prepf !*kk2f list('arbconst,j),newde)
- % newde:=reval subst(newco,list('arbconst,j),newde)
- % newde:=reval subst(newco, prepf !*kk2f list('arbconst,j),newde)
- end$
- newde>>
- end$
- endmodule$
- %********************************************************************
- module simplifications$
- %********************************************************************
- % Routines for simplifications and substitution of functions
- % Author: Andreas Brand
- % August 1991
- symbolic procedure clearfactors(sorg,ftem,vl)$
- % liefert Liste der Gl. die entstehen, wenn in allen Gl. aus sorg
- % gemeinsame Faktoren, die keine der Funktionen aus ftem sowie deren
- % Ableitungen enthalten, gekuerzt werden
- if sorg and ftem then
- begin scalar l,ft$
- while sorg do
- <<if not zerop car sorg then
- if not(ft:=smemberl(ftem,car sorg)) then l:=cons(car sorg,l)
- else l:=cons(reval algsimple(reval car sorg,car ft,ft,vl),l)$
- sorg:=cdr sorg>>$
- return l
- end
- else sorg$
- symbolic procedure algsimple(g,f,ftem,vl)$
- % Gl. die entsteht, wenn in g gemeinsame Faktoren, die keine der Fkt.
- % aus ftem sowie deren Ableitungen enthalten, gekuerzt werden
- if pairp g then
- begin scalar h,l$
- if car g='MINUS then g:=algsimple(cadr g,f,ftem,vl)
- else if car g='QUOTIENT then
- g:=reval aeval list('QUOTIENT,algsimple(cadr g,f,ftem,vl),
- algsimple(caddr g,f,ftem,vl))
- else
- <<h:=ftem$
- if pairp g and (car g='PLUS) and (pairp cadr g) and (caadr g='MINUS)
- then g:=reval list('MINUS,g)$
- if freeof(g,f) then
- <<f:=nil$
- while h do
- if not freeof(g,car h) then <<f:=car h$h:=nil>>
- else h:=cdr h>>$
- if f then
- <<l:=lderiv2(g,f,fctargs f)$
- h:=car reverse cdr reval list('COEFF,g,reval car l)$
- %highest coefficient
- if lowpow!*>1 then
- g:=reval list('QUOTIENT,g,list('EXPT,car l,lowpow!*-1))$
- if not zerop h then
- if pairp h then
- <<h:=independpart(reval list('GCD,g,h),ftem,vl)$
- if not zerop h then g:=reval list('QUOTIENT,g,h)>>
- else g:=reval list('QUOTIENT,g,list('GCD,g,h))>>
- >>$
- return g
- end
- else if zerop g then 0
- else if freeof(ftem,g) then 1
- else g$
- symbolic procedure lcasepossible(l,vl)$
- begin scalar l1$
- for each a in l do l1:=union(list casepossible(a,vl),l1)$
- return delete('NIL,l1)
- end$
- symbolic procedure casepossible(a,vl)$
- if pairp a then
- if (car a='QUOTIENT) or (car a='EXPT) or (car a='MINUS) then
- casepossible(cadr a,vl)
- else
- if car a='TIMES then begin scalar l$
- return if (l:=lcasepossible(cdr a,vl)) then reval cons('TIMES,l)
- end
- else if setdiff(idlist a,append(idlist vl,REDUCEFUNCTIONS_)) then a
- else if setdiff(idlist a,append(idlist vl,REDUCEFUNCTIONS_)) then a$
- symbolic procedure casecheck(p,vl)$
- begin scalar a$
- if (a:=casepossible(p,vl)) and not member(a,special_cases) then
- <<special_cases:=cons(a,special_cases)$
- ineq_:=union(list a,ineq_)$
- if print_ then
- <<terpri()$write "Special case to be considered : "$
- deprint list a>> >>$
- return a
- end$
- symbolic procedure independpart(g,ftem,vl)$
- % determines factors not depending on functions in ftem
- begin scalar l,l1$
- g:=cdr reval list('FACTORIZE,g)$
- while g do
- <<if member(car g,cdr g) or member(car g,ineq_)
- then l:=cons(car g,l)
- else l1:=cons(car g,l1)$
- g:=cdr g>>$
- for each x in l1 do
- if not smemberl(ftem,x) then
- if not casecheck(x,vl) or not sp_cases then l:=cons(x,l)$
- return if l then if cdr l then cons('TIMES,l)
- else car l
- else 1
- end$
- symbolic procedure splitde(p,ftem)$
- begin scalar l$
- if pairp p and (car p='QUOTIENT) then p:=cadr p$
- for each x in cdr reval list('FACTORIZE,p) do
- if smemberl(ftem,x) then l:=cons(reval x,l)$
- return if length l>1 then union(l,nil)
- else nil
- end$
- symbolic procedure simplifyde(a,ftem,vl,genflag)$
- begin scalar l,f,p$
- ftem:=smemberl(ftem,a)$
- if f:=fctchoose(list a,ftem) then
- <<p:=reval algsimple(reval a,f,ftem,vl)$
- if not p or zerop p then
- <<write "***** Error in ALGSIMPLE : "$
- mathprint a>>
- else a:=p$
- if not stardep(a,ftem,vl) then
- if (l:=integratede(a,ftem,genflag)) and print_ then
- if null cdr l then
- <<terpri()$write "integrated equation : "$deprint(l)>>
- else
- <<terpri()$write "generalized integration yields ",
- length l," equations : "$deprint(l)>> >>$
- return if l then l
- else list a
- end$
- symbolic procedure fctsubst(ex,fo,forg)$
- % substitution of a function fo in the list forg by an expression ex
- <<%if not member(fo,forg) then depl!*:=delete(assoc(fo,depl!*),depl!*)$
- for each f in forg collect
- if f=fo then list('EQUAL,fo,ex)
- else if freeof(f,fo) then f
- else reval subst(ex,fo,f)>>$
- symbolic procedure substandsep(g,sorg,ftem,vlist)$
- % in allen Gl. aus sorg wird die Fkt. car g durch cdr g ersetzt
- begin scalar l,p,q,a$
- while sorg do
- <<p:=car sorg$
- sorg:=cdr sorg$
- if freeof(p,car g) then l:=union(list p,l)
- else
- <<q:=reval subst(cdr g,car g,p)$
- if (ineq_:=ineqsubst(g,ftem)) then
- if pairp q then
- <<a:=gensep(q,smemberl(union(fnew_,ftem),q),vlist)$
- if contradiction_ then l:=sorg:=nil
- else l:=union(car a,l)>>
- else
- if q neq 0 then l:=union(list q,l)>>
- >>$
- return l
- end$
- symbolic procedure fcteval(p,ftem,vl)$
- begin
- if pairp p and (car p='QUOTIENT) then
- <<casecheck(caddr p,vl)$
- p:=cadr p>>$
- return
- if delength p>(if fcteval_ then fcteval_ else 0) then
- (if cont_ then
- if yesp("function to be evaluated ? ") then
- fcteval1(p,ftem,vl))
- else fcteval1(p,ftem,vl)
- end$
- symbolic procedure fcteval1(p,ftem,vl)$
- begin scalar l,f,n,ft,a$
- ft:=smemberl(ftem,p)$
- n:=length varslist(p,ft,vl)$
- for each x in ft do
- if fctlength x>=n then
- <<l:=ldiffp(p,x)$
- if (null car l) and (cdr l=1) then
- if freeofzero(coeffn(p,x,1),delete(x,ft)) then f:=x>>$
- if f then
- p:=reval list('DIFFERENCE,f,list('QUOTIENT,p,coeffn(p,f,1)))$
- return if f then cons(p,f)
- else nil$
- end$
- symbolic procedure freeofzero(p,ft)$
- % liefert t, falls p nicht 0 wird
- if null ft then t
- else if freeof(p,car ft) then freeofzero(p,cdr ft)
- else
- begin scalar a$
- p:=cdr reval list('FACTORIZE,p)$
- while p do
- if null smemberl(ft,car p) or member(car p,ineq_) then p:=cdr p
- else <<p:=nil$
- a:=t>>$
- return (null a)
- end$
- symbolic procedure fctevalstrict(a,ftem,vl)$
- % fcteval with additional condition
- % substituted expr. contains only functions of less arguments
- % than the function
- begin scalar l,n,ft$
- l:=fcteval(a,ftem,vl)$
- if l then
- <<ft:=smemberl(ftem,car l)$
- n:=for each f in ft collect length fctargs f$
- if not (reval cons('MAX,n)<length fctargs cdr l) then l:=nil>>$
- return l$
- end$
- endmodule$
- %********************************************************************
- module utilities$
- %********************************************************************
- % Routines for finding leading derivatives and others
- % Author: Andreas Brand
- % June 1990
- symbolic procedure diffrel(p,q,v)$
- % liefert komplizierteren Differentialausdruck$
- if diffrelp(p,q,v) then q
- else p$
- symbolic procedure diffrelp(p,q,v)$
- % liefert q, falls p einfacherer Differentialausdruck, sonst nil
- % p, q Paare (liste.power), v Liste der Variablen
- % liste Liste aus Var. und Ordn. der Ableit. in Diff.ausdr.,
- % power Potenz des Differentialausdrucks
- if null v then % alle Variable untersucht ?
- if cdr p>cdr q then nil
- else t
- else begin
- scalar a,b$
- a:=diffdeg(car p, car v)$ % Ordnung der Ableitung nach
- b:=diffdeg(car q,car v)$ % der ersten Variablen
- return if a<b then t
- else if b<a then nil
- else diffrelp(p,q,cdr v) % falls Ableitungen
- % erste Variable gleich, dann
- % restliche Variablen
- end$
- symbolic procedure diffdeg(p,v)$
- % liefert Ordnung der Ableitung von p nach v$
- % p Liste Varible+Ordnung der Ableitung, v Variable (Atom)
- if null p then 0 % alle Variable bearbeitet ?
- else if car p=v then % v naechste Variable ?
- if cdr p then
- if numberp(cadr p) then cadr p % folgt eine Zahl ?
- else 1
- else 1
- else diffdeg(cdr p,v)$ % weitersuchen
- symbolic procedure ldiff(p,f)$
- % Suchen der fuehrenden Ableitung der Fkt. f(arg1.,...,argn) in p
- % p Ausdruck in Listenform,f Funktion in listenform
- % Ergebnis: (p,(liste.power))$ liste: Liste der Ordn. der Ableitungen
- begin scalar l$
- l:=ldiffp(p,f)$ % fuerende Liste mit Potenz
- return cons(ldiff1(car l,fctargs f),cdr l)$
- % aus Liste Variablen + Ordnung
- % wird Liste der Ordnungen
- end$
- symbolic operator fargs$
- symbolic procedure fargs f$
- cons('LIST,fctargs f)$
- symbolic procedure fctargs f$
- % arguments of a function
- %if pairp f then cdr f$
- if (f:=assoc(f,depl!*)) then cdr f$
- symbolic procedure fctlength f$
- % number of arguments
- length fctargs f$
- symbolic procedure argset(ftem)$
- % List of arguments of all functions in ftem
- if ftem then union(reverse fctargs car ftem,argset(cdr ftem))
- else nil$
- symbolic procedure ldiff1(l,v)$
- % liefert Liste der Ordnungen der Ableitungen nach den Variablen aus v
- % l Liste (Variable + Ordnung)$ v Liste der Variablen
- if null v then nil % alle Variable abgearbeitet ?
- else cons(diffdeg(l,car v),ldiff1(l,cdr v))$
- % Ordnung der Ableitung nach
- % erster Variable anhaengen
- symbolic procedure ldiffp(p,f)$
- % liefert Liste der Variablen + Ordnungen mit Potenz
- % p Ausdruck in LISP - Notation, f Funktion
- ldiffp1(p,f,fctargs f)$
- symbolic procedure ldiffp1(p,f,vl)$
- % liefert Liste der Variablen + Ordnungen mit Potenz
- % p Ausdruck in LISP - Notation, f Funktion, lv Variablenliste
- begin scalar a$
- a:=cons(nil,0)$
- if not atom p then
- if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,'QUOTIENT,'DF,'EQUAL))
- then
- % erlaubte Funktionen
- <<if (car p='PLUS) or (car p='TIMES) or (car p='QUOTIENT)
- or (car p='EQUAL) then
- <<p:=cdr p$
- while p do
- <<a:=diffrel(ldiffp1(car p,f,vl),a,vl)$
- p:=cdr p>> >>
- else if car p='MINUS then
- a:=ldiffp1(cadr p,f,vl)
- else if car p='EXPT then % Exponent
- if numberp caddr p then
- <<a:=ldiffp1(cadr p,f,vl)$
- % fuehrende Abl. aus der Basis
- a:=cons(car a,times(caddr p,cdr a))>>
- else a:=cons(nil,0)
- % Poetenz aus Basis wird mit
- % Potenz multipliziert
- else if car p='DF then % Ableitung
- if cadr p=f then a:=cons(cddr p,1)
- % f wird differenziert?
- else a:=cons(nil,0)>> % sonst Konstante bzgl. f
- else if p=f then a:=cons(nil,1)
- % Funktion selbst
- else a:=cons(nil,0) % alle uebrigen Fkt. werden wie
- else if p=f then a:=cons(nil,1)$ % Konstante behandelt
- return a
- end$
- symbolic procedure lderiv(p,f)$
- % fuehrende Ableitung in LISP-Notation ohne Potenz
- lderiv1(p,f,fctargs f)$
- symbolic procedure lderiv1(p,f,vl)$
- % fuerende Ableitung in LISP-Notation ohne Potenz
- % mit Angabe der Rangfolge der Variablen
- begin scalar l$
- l:=ldiffp1(p,f,vl)$
- return if car l then cons('DF,cons(f,car l))
- else if zerop cdr l then nil
- else f
- end$
- symbolic procedure lderiv2(p,f,vl)$
- % fuehrende Ableitung in LISP-Notation mit Potenz (als dotted pair)
- begin scalar l$
- l:=ldiffp1(p,f,vl)$
- return cons(if car l then cons('DF,cons(f,car l))
- else if zerop cdr l then nil
- else f
- ,cdr l)
- end$
- symbolic procedure highdiff(p,f,x)$
- % hoechste Ableitung einer Funktion f nach der Variable x
- % in einem Ausdruck p
- % p,f in LISP-Notation, x Atom
- diffdeg(car ldiffp1(p,f,list x),x)$
- symbolic procedure minausd(p,q)$
- % minimalen Differentialausdruck suchen$
- % p,q, Ausdr. in Listenform mit fuehrender Ableitung + Potenz
- % zur Berechnung wird nur fuehrender Ableitung + Potenz benutzt
- if minausdp1(cadr p,cddr p,cadr q,cddr q) then p
- else q$
- symbolic procedure minausdsp(p,q)$
- % falls p minimaler Diff.ausdr. und p nicht laenger als q, p sonst nil
- % p,q, Ausdr. in Listenform mit fuerender Ableitung und Potenz
- % benutzt wird nur fuerende Ableitung und Potenz
- begin scalar s$
- s:=minausdp1(cadr p,cddr p,cadr q,cddr q)$
- return if s then if s='EQUAL!= then % falls die Ausdruecke gleich
- % schwierig sind,
- if shorter(car p,car q) then p
- else nil
- % geht auch die Laenge ein
- else p
- else nil
- end$
- symbolic procedure minausdp(p,q)$
- % falls p minimaler Differentialausdruck, p sonst nil
- % p,q, Ausdr. in Listenform mit fuehrender Ableitung und Potenz
- % benutzt wird nur fuerende Ableitung und Potenz
- begin scalar j,k$
- return if minausdp1(cadr p,cddr p,cadr q,cddr q) then p
- else nil$
- end$
- symbolic procedure minausdp1(l1,p1,l2,p2)$
- % liefert t, wenn der erste von zwei Differentialausdr. echt einfacher,
- % 'EQUAL!=, wenn sie gleich schwierig sind, sonst nil
- % l1,l2 Listen der Ordnungen der Ableit., p1,p2 Potenzen
- if null l1 then % Ordnungen gleich ?
- if p1=p2 then 'EQUAL!=
- % Potenzen gleich ?
- else p1<p2
- else if null l2 then nil % darf eigentlich nicht sein
- else if car l1<car l2 then t % Ordn. der 1. Var. kleiner ?
- else if car l1>car l2 then nil
- % Ordn. der 1. Var. groeaer ?
- else minausdp1(cdr l1,p1,cdr l2,p2)$
- % restliche Variable testen
- symbolic procedure difdiff(r,s)$
- % liefert Liste der Differenzen der Elemente zweier Listen
- if null r then nil
- else cons (car r-car s,difdiff(cdr r,cdr s))$
- symbolic procedure mkldiff(f,l)$
- % erzeugt die Ableitung der Funktion f nach
- % den Variablen (einschliesslich Ordnungen) aus der Liste l
- cons('DF,cons(f,difflist(fctargs f,l)))$
- symbolic procedure difflist(v,l)$
- % erzeugt aus der Liste der Variablen und der Liste der Ordnungen der
- % Ableitungen eine Liste von Variablen + Ordnungen
- if null v then nil
- else if car l=0 then difflist(cdr v,cdr l)
- % Ordnung 0
- else if car l=1 then cons(car v,difflist(cdr v,cdr l))
- % bei Ordnung 1 wird Variable
- % ohne Ordnung angehaengt
- else cons(car v,cons(car l,difflist(cdr v,cdr l)))$
- symbolic procedure delength(d)$
- % Laenge eines Polynoms in LISP - Notation
- if not pairp d then if d then 1
- else 0
- else
- if car d='PLUS then length(d)-1 % Laenge in LISP ohne PLUS
- else if car d='EQUAL then delength(cadr d)+delength(caddr d)
- else if car d='QUOTIENT then delength(cadr d)
- else 1$
- symbolic procedure shorter(a,b)$
- delength a<=delength b$
- symbolic procedure desort(l)$
- % list sorting
- begin scalar l1,l2,l3,m,n$
- return
- if null l then nil
- else
- <<n:=delength car l$
- l2:=list car l$
- l:=cdr l$
- while l do
- <<m:=delength car l$
- if m<n then l1:=cons(car l,l1)
- else if m>n then l3:=cons(car l,l3)
- else l2:=cons(car l,l2)$
- l:=cdr l>>$
- append(desort(l1),append(l2,desort(l3)))>>
- end$
- symbolic procedure smemberl(fl,ex)$
- if fl and ex then
- if smember(car fl,ex) then cons(car fl,smemberl(cdr fl,ex))
- else smemberl(cdr fl,ex)$
- symbolic procedure fctsort(l)$
- % list sorting
- begin scalar l1,l2,l3,m,n$
- return
- if null l then nil
- else
- <<n:=fctlength car l$
- l2:=list car l$
- l:=cdr l$
- while l do
- <<m:=fctlength car l$
- if m<n then l1:=cons(car l,l1)
- else if m>n then l3:=cons(car l,l3)
- else l2:=cons(car l,l2)$
- l:=cdr l>>$
- append(fctsort(l1),append(l2,fctsort(l3)))>>
- end$
- symbolic procedure fctprint(fl)$
- % Ausdrucken der Funktionen aus fl
- begin scalar n,l$
- n:=0$
- for each f in fl do
- if pairp f then
- if car f='EQUAL then
- <<n:=delength reval f$
- if n>print_ then
- <<terpri()$write cadr f,"= expr. with ",n," terms"$terpri()>>
- else mathprint f$
- n:=0>>
- else
- <<if n eq 5 then
- <<terpri()$n:=0>>$
- write car f$
- if pairp cdr f then
- <<write "("$write cadr f$
- for each v in cddr f do
- <<write ","$write v>>$
- write ") ">>
- else
- <<write car f$
- write " ">>$
- n:=add1 n>>
- else
- <<if n eq 5 then
- <<terpri()$n:=0>>$
- write f$
- if (l:=fctargs f) then
- <<write "("$write car l$
- for each v in cdr l do write ",",v$
- write ")">>$
- write " "$
- n:=add1 n>>$
- %if n neq 0 then terpri()
- end$
- symbolic procedure eqprint(e)$
- % Ausdrucken der Gl. e
- if print_ then
- begin scalar n$
- n:=delength reval e$
- if n>print_ then
- <<write "expr. with "$write n$write " terms"$terpri()>>
- else
- mathprint e$
- end$
- symbolic procedure substprint(a,b)$
- if print_ then
- begin scalar n$
- n:=delength reval b$
- if n>print_ then <<fctprint list a$write "= expr. with ",n," terms"$
- terpri()>>
- else mathprint list('EQUAL,a,b)
- end$
- symbolic procedure deprint(l)$
- % Ausdrucken der Gl. aus der Liste l
- if l and print_ then for each x in l do eqprint(list('EQUAL,0,x))$
- symbolic procedure powappend(l,f)$
- % Anhaengen der fuehrenden Ableitung
- if null l then nil
- else cons(cons(car l,ldiff(car l,f)),powappend(cdr l,f))$
- symbolic procedure powdelete(l)$
- % Loeschen der fuehrenden Ableitung
- if null l then nil
- else cons(caar l,powdelete(cdr l))$
- symbolic procedure fctdepend(e,v,fctset)$
- % Test,ob variable v im Ausdruck e nicht in Abhaeng. von
- % Fkt. aus fctset vorkommt
- if member(v,argset smemberl(fctset,e)) then t
- else freeof(e,v) $
- symbolic procedure subset(a,b)$
- % test,ob a Teilmenge von b ist
- if null a then t
- else if member(car a,b) then subset(cdr a,b)
- else nil$
- symbolic procedure idlist(l)$
- if pairp l then union(idlist car l,idlist cdr l)
- else if l and idp l then list l
- else nil$
- symbolic procedure varapp(l,a,v1,v2)$
- % an jede Gl. aus l werden v1 und v2 angehaengt
- if null l then nil
- else
- cons(list(cons(cons(caar l,a),cdar l),v1,v2),varapp(cdr l,a,v1,v2))$
- symbolic procedure varslist(p,ftem,vl)$
- begin scalar l$
- ftem:=argset smemberl(ftem,p)$
- for each v in vl do
- if not freeof(p,v) or member(v,ftem) then l:=cons(v,l)$
- return reverse l$
- end$
- symbolic procedure newfct(id,l,nfct)$
- begin scalar f$
- f:=mkid(id,nfct)$
- depl!*:=delete(assoc(f,depl!*),depl!*)$
- %put(f,'simpfn,'simpiden)$
- %if pairp l then f:=cons(f,l)$
- if pairp l then depl!*:=cons(cons(f,l),depl!*)$
- if print_ then
- <<terpri()$
- if pairp l then
- <<write "new function: "$
- fctprint list f>>
- else
- write "new constant: ",f>>$
- return f$
- end$
- symbolic procedure equallist(forg)$
- % List of all members of forg which are pairs
- if forg then
- if pairp car forg then cons(car forg,equallist cdr forg)
- else equallist cdr forg$
- symbolic procedure spmin(a,b)$
- if null a then b
- else if null b then a
- else if cadr a<cadr b then a
- else b$
- symbolic procedure specons(a,l)$
- if freeof(l,car a) then cons(car a,cons(cadr a,l))
- else specons1(a,l)$
- symbolic procedure specons1(a,l)$
- if null l then nil
- else if car l=car a then
- if cadr l<cadr a then cons(car a,cons(cadr a,cddr l))
- else l
- else cons(car l,cons(cadr l,specons1(a,cddr l)))$
- symbolic procedure polyp(p,ftem)$
- begin scalar a$
- a:=t$
- while ftem do
- if polyp1(p,car ftem) then ftem:=cdr ftem
- else
- <<terpri()$
- write "***** This equation is not a polynomial w.r.t. "$
- write car ftem$write " and its derivatives : "$
- deprint(list p)$
- ftem:=nil$
- a:=nil>>$
- return a
- end$
- symbolic procedure polyp1(p,f)$
- % prueft, ob p f nur in polynomialer Form enthaelt
- if freeof(p,f) then t
- else
- begin scalar a$
- if atom p then a:=t
- else
- if member(car p,list('EXPT,'PLUS,'MINUS,'TIMES,'QUOTIENT,'DF)) then
- % erlaubte Funktionen
- <<if (car p='PLUS) or (car p='TIMES) then
- <<p:=cdr p$
- while p do
- if a:=polyp1(car p,f) then p:=cdr p
- else p:=nil>>
- else if (car p='MINUS) then
- a:=polyp1(cadr p,f)
- else if (car p='QUOTIENT) then
- <<if freeof(caddr p,f) then a:=polyp1(cadr p,f)>>
- else if car p='EXPT then % Exponent
- <<if (fixp caddr p) then
- if caddr p>0 then a:=polyp1(cadr p,f)>>
- else if car p='DF then % Ableitung
- if (cadr p=f) or freeof(cadr p,f) then a:=t>>
- else a:=(p=f)$
- return a
- end$
- symbolic operator rationalexp$
- symbolic procedure rationalexp(p,ftem)$
- ratexp(p,cdr ftem)$
- lisp flag('(rationalexp),'BOOLEAN)$
- symbolic procedure ratexp(p,ftem)$
- if null ftem then t
- else if ratexp1(p,car ftem) then ratexp(p,cdr ftem)
- else nil$
- symbolic procedure ratexp1(p,f)$
- % prueft, ob p f nur in rationaler Form enthaelt
- if not pairp p or freeof(p,f) then t
- else if car p='QUOTIENT then polyp1(cadr p,f) and polyp1(caddr p,f)
- else polyp1(p,f)$
- %symbolic procedure freeof(u,v)$
- %not(smember(v,u)) and freeofdepl(depl!*,u,v)$
- symbolic procedure freeofdepl(de,u,v)$
- if null de then t
- else if smember(v,cdar de) and smember(caar de,u) then nil
- else freeofdepl(cdr de,u,v)$
- symbolic procedure freeoflist(l,m)$
- % liefert t, falls kein Element aus m in l auftritt
- if null m then t
- else if freeof(l,car m) then freeoflist(l,cdr m)
- else nil$
- endmodule$
- %********************************************************************
- module dfint; % Patch to improve symbolic differentiation,
- % mainly of integrals.
- % Author: Francis J. Wright <F.J.Wright@QMW.ac.uk>
- % QMW, London and CBPF, Rio.
- % Date: Thu Nov 19 21:11:12 1992
- % REDUCE version: 3.4.1 (ONLY); PSL
- % Beta test version
- % 19 June 1992: Preliminary alpha test version (for CRACK)
- % 28 Oct 1992: !*k2q -> !*kk2q to avoid non-unique kernels
- % But maybe it should be done as in diffp instead?
- % 7 Nov 1992: Allowdfint mode added; completely revised REDUCE 3.4.1
- % version
- % DOCUMENTATION
- % =============
- % This patch has several functions: for examples see the test and
- % demonstration file DFINT.TST.
- % Let A, B, C, D be arbitrary SEQUENCES of kernels.
- % 1. REDUCE commutes the order of differentiations so as to collect
- % differentiation variables, but it misses the commutation that would
- % simplify any expression equivalent to df(y, A, y, B) to zero (unless
- % A is absent). A one-line patch to procedure diffp (below) fixes
- % this. This simplification must surely ALWAYS be correct and
- % desirable (since REDUCE commutes differentiations), and is very
- % cheap to include, so it is ALWAYS performed.
- % 2. However, the main purpose of this patch is to differentiate
- % integrals (hence the name DFINT) in a controlled manner. This is
- % implemented for single derivatives and differentiation wrt the first
- % variable in multiple derivatives via the new main procedure
- % dfform_int, together with a patch to the existing procedure diffp
- % (below) to deal with differentiation wrt later variables in multiple
- % derivatives.
- % It is not possible in general in REDUCE to perform any of these
- % simplifications by pattern matching because of the arbitrary number
- % of arguments. In particular, this patch makes redundant the
- % following specific rule defined in INT.RED:
- % put('df,'opmtch,'(((int !=y !=x) !=x) (nil . t) (evl!* !=y) nil) .
- % get('df,'opmtch));
- % This is equivalent (but NOT identical) to the algebraic command
- % for all y,x let df(int(y,x), x) = y;
- % which could therefore be cleared - see DFINT.TST. However, leaving
- % it effective should not matter. The code here to replace this rule
- % is not NECESSARY, but seems desirable since the code has to check
- % for differentiation of an integral anyway.
- % The precise operation of the new facility to differentiate integrals
- % is as follows:
- % (a). Differentiation wrt the integration variable:
- % df(int(F(A,x,B),x),C,x,D), or any equivalent composition of
- % derivatives applied to one integral, is ALWAYS simplified to
- % df(F(A,x,B),C,D).
- % (b). Differentiation wrt a parameter of the integrand:
- % If the switch ALLOWDFINT is turned on (it is ON by default) then "df
- % int" will be commuted into "int df" IF THIS SIMPLIFIES THE
- % INTEGRAND, meaning currently that the new integrand is not a
- % symbolic derivative. Thus differentiation under the integral sign
- % will be performed so that df(int(F(A,x,B),x),C,y,D) simplifies to
- % df(int(F_y(A,x,B),x),C,D) if F_y = dF/dy simplifies. It think this
- % is normally desirable.
- % If the switch DFINT is turned on (it is OFF by default) then the
- % commutation is ALWAYS performed. Thus df(int(F(A,x,B),x),C,y,D)
- % ALWAYS evaluates to int(df(F(A,x,B),C,y,D),x). I think this will
- % not normally be desirable, although for special purposes it is.
- % If used together with either Herbert Melenk's patch or mine (called
- % intpatch) that allow integrands to be implicitly dependent on the
- % integration variable, this code is intended to also work with such
- % dependent integrands, rather than only with explicit operators as
- % illustrated above. For examples see the end of the test and
- % demonstration file DFINT.TST.
- % This file is intended to be compiled using faslout and then loaded
- % using load_package, but alternatively the source code can be input
- % when required and either compiled or interpreted. It could also be
- % edited into the standard source file POLY.RED and built permanently
- % into REDUCE.
- % ** WARNING ** To use this patch with REDUCE 3.4 the new procedures
- % called by diffp in 3.4.1 need to be supplied, or the end of diffp
- % changed back to its 3.4 form.
- switch allowdfint, dfint; % dfint OFF by default
- deflist('((allowdfint ((t (rmsubs)))) (dfint ((t (rmsubs))))), 'simpfg);
- % There is no code to reverse the df-int commutation,
- % so no reason to call rmsubs when the switch is turned off.
- % on allowdfint; % But would be done when INput rather than LOADed!
- !*allowdfint := t; rmsubs(); % allowdfint ON by default
- !*dfint:=t;
- put('int, 'dfform, 'dfform_int);
- % dfform is a new hook, otherwise used only by taylor (I think).
- % This code 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.
- 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)
- begin scalar result, x, y;
- y := simp!* cadr u; % SQ form integrand
- x := caddr u; % kernel
- result :=
- if v eq x then y
- % df(int(y,x), x) -> y replacing the let rule in INT.RED
- else if (!*allowdfint or !*dfint) and << y := diffsq(y, v);
- !*dfint or not_df_p y % it has simplified
- >> then simp{'int, mk!*sq y, x} % MUST re-simplify it!!!
- % i.e. differentiate under the integral sign
- % df(int(y, x), v) -> int(df(y, v), x).
- % (Perhaps I should use prepsq - kernels are normally
- % true prefix?)
- else !*kk2q{'df, u, v}; % remain unchanged
- if not(n eq 1) then
- result := multsq( (((u .** (n-1)) .* n) .+ nil) ./ 1, result);
- return result
- end;
- symbolic procedure not_df_p y;
- % True if the SQ form y is not a df kernel.
- not(denr y eq 1 and
- not domainp (y := numr y) and eqcar(mvar y, 'df));
- % Author: Anthony C. Hearn.
- % Copyright (c) 1991 RAND. All rights reserved.
- % Modifications by Francis J. Wright, QMW, London and CBPF, Rio
- % (clearly flagged -- search for FJW).
- % This procedure is from the REDUCE 3.4.1 module poly in file POLY.RED
- % and is a little different from the 3.4 version! It needs a minor
- % addition to improve the simplification of MULTIPLE derivatives.
- % Note that in PSL-REDUCE POLY is compiled into the base system, so
- % this procedure is always defined and can safely be re-defined.
- % Under other Lisps this may not be true, so BEWARE!!!
- 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) 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))
- then return nil ./ 1;
- w := list('df,u,v);
- go to j;
- h: % Final check for possible kernel deriv.
- if car u eq 'df
- then if depends(cadr u,v)
- %% BEGIN addition to 3.4.1 version by FJW
- and not(cadr u eq v) then
- % (df (df v A) v) ==> 0
- << 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)
- % 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 := aeval{'int, mk!*sq w, caddr cadr u} . cddr u;
- go to j >>; % derivative absorbed
- %% END addition to 3.4.1 version by FJW
- %% BEGIN difference from 3.4
- if (x := find_sub_df(w:= cadr u . derad(v,cddr u),
- 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
- %% END difference from 3.4
- >> %FJW
- else return nil ./ 1
- else if depends(u,v) then w := list('df,u,v)
- else return nil ./ 1;
- j: w := if x := opmtch w then simp x else mksq(w,1);
- go to e
- end;
- endmodule;
- module intpatch; % Integrate dependent variables & rational powers
- % Author: Francis J. Wright <F.J.Wright@QMW.ac.uk>
- % Date: 19 June 1992 (Very minor update 16/11/92)
- % REDUCE version: 3.4 or 3.4.1; PSL
- % This version (19 June 1992) fixes a bug that integrals that remained
- % symbolic were not returned as unique kernels - I now use !*kk2q.
- % This bug caused rather obscure symptoms, such as failures with on
- % factor.
- % DOCUMENTATION
- % =============
- % This patch has two separate functions:
- % 1: It allow integrals containing IMPLICITLY dependent variables,
- % created using the DEPEND command, to remain symbolic rather than
- % cause an error, whilst preserving other error handling as normal.
- % ON FAILHARD turns this facility off.
- % This facility was developed from a patch by Herbert Melenk,
- % which this patch is intended to replace.
- % 2: It integrates simple rational powers of the integration
- % variable that the integrator currently fails to integrate.
- % This file is intended to be compiled using faslout and then loaded
- % using load_package, but alternatively the source code can be input
- % when required and either compiled or interpreted. It could also be
- % edited into the standard source file INT.RED and built permanently
- % into REDUCE.
- % Because of the autoload definitions for int (see ENTRY.RED) it is
- % safe to load this module before the integrator has been loaded, and
- % it does not cause the integrator to load. The integrator can safely
- % be AUTOLOADED by calling INT afterwards, either from algebraic mode
- % or from symbolic mode using any of the methods described below, but
- % MUST NOT BE EXPLICITLY LOADED AFTER this patch, otherwise this patch
- % will not work because the simpfn property of int will get
- % overwritten. I cannot see any solution to this, other than avoiding
- % doing it!
- % WARNING: To call the integrator from symbolic mode when this patch
- % is loaded DO NOT CALL SIMPINT DIRECTLY. Instead, there are several
- % possibilities that should work, including any of the following:
- % 1. Use the form ????{'int, integrand, int_var}, where ???? is
- % reval, aeval or simp, etc, depending on the form of result you want,
- % and the arguments are in prefix (or pseudo-prefix) form. This is
- % equivalent to calling the integrator from algebraic mode, and I
- % currently think this is probably the best approach.
- % 2. Call simpintpatch directly, but this is not general and so not
- % recommended.
- % 3. Replace calls of "simpint" by (say) "integrator", and do
- % essentially the following:
- % global'(integration_function);
- % symbolic(integration_function := get('int, 'simpfn));
- %
- % symbolic procedure integrator u;
- % apply1(integration_function, u);
- % Note: apply1 is defined in RLISP.RED.
- fluid '(!*failhard);
- %% fluid '(soft!-rerror!-number);
- % Hope not necessary - it seems not to be.
- symbolic procedure SimpIntPatch u;
- % Driver for various patches:
- % 1: Catch errors from SimpInt, trap error number 7 only,
- % and pass on all other errors as normal hard REDUCE errors.
- % 2: Post-process unintegrated rational powers.
- begin scalar r, !*usermode,!*redefmsg, !*uncached; !*uncached := t;
- % !*redefmsg rebound to avoid PSL messages
- % about redefinition of rerror.
- %% integer soft!-rerror!-number; % defaults to 0, not nil
- put('int, 'simpfn, 'SimpInt); % assumed & reset by SimpInt
- copyd('rerror, 'intpatch!-rerror); % redefine rerror
- r := errorset!*({'SimpInt, mkquote u}, nil);
- copyd('rerror, 'rerror!%); % restore rerror
- put('int, 'simpfn, 'SimpIntPatch); % reset INT interface
- if pairp r then <<
- % First call of SimpInt succeeded -
- % try to reprocess any integrals left:
- put('int, 'simpfn, 'SimpIntRatPow);
- u := resimp car r; % this works ONLY with !*uncached := t;
- put('int, 'simpfn, 'SimpIntPatch);
- return u
- >>
- else if !*failhard or not(r eq 7) then
- rederr EMSG!* % Error failure
- else return !*kk2q('int . u) % Remain symbolic
- end;
- % Integrator error trap patch to allow controlled error handling
- % ==============================================================
- % The error numbers generated by SIMPINT and the corresponding
- % error message that would be output by INT are the following,
- % collected from the INT source code:
- % 1 = "Improper number of arguments to INT"
- % 2 = "Improper number of arguments to INT"
- % 3 = "Too many arguments to INT"
- % 4 = "FAILHARD switch set"
- % 5 = "Invalid polynomial in int-quadterm"
- % 6 = "Empty list to mapply"
- % 7 = "Can't integrate in the presence of side-relations" (TRAPPED)
- % 8 = "Invalid exponent"
- % 9 = "FAILHARD switch set"
- % If any other error number, such as 0, should occur then it
- % corresponds to some other non-specific error.
- symbolic procedure rerror!%(packagename,number,message);
- % This is precisely the definition of rerror in RLISP.RED,
- % but redefining it here makes sure it is loaded,
- % and also avoids the need to save it.
- % Precisely this procedure is also defined in SOFTSOLV.
- rederr message;
- symbolic procedure intpatch!-rerror(packagename,number,message);
- %% << soft!-rerror!-number := number; error1() >>;
- % The following will suffice provided errorsets
- % are not nested in the integrator.
- % It makes error message text available in EMSG!*.
- error(number, message);
- % Integrator postprocessor patch to integrate simple rational
- % powers that the integrator currently fails to integrate.
- % =======================================================
- symbolic procedure SimpIntRatPow u; % u = (integ var)
- % Integrate integrands of the form var**(m/n),
- % which the integrator leaves in a rather bizarre form -
- % hence the precise form of the following code.
- % Returns original integral if it has the wrong form.
- begin scalar integ, var, power;
- integ := car u; var := cadr u;
- % assumes true prefix forms, already evaluated by SimpInt.
- % power := errorset!*(
- % {'FindRatPow, mkquote integ, mkquote var}, nil);
- % errorset!*(u,v) == errorset(u,v,!*backtrace)
- % Backtrace from this is unlikely to be interesting, so ...
- power := errorset(
- {'FindRatPow, mkquote integ, mkquote var}, nil, nil);
- if errorp power then return !*kk2q('int . u);
- power := car power; % correct form of integrand found.
- % integrand = var**power, so return integral:
- power := reval {'plus, power, 1};
- return simp!* {'quotient, {'expt, var, power}, power}
- end;
- symbolic procedure FindRatPow(monom, var);
- % Return power of a monomial in var, as a
- % rational number in UNSIMPLIFIED prefix form
- % or cause error return to enclosing errorset.
- if eqcar(monom, 'quotient) then
- {'plus, FindRatPow(cadr monom, var),
- {'minus, FindRatPow(caddr monom, var)}}
- else if eqcar(monom, 'times) then
- 'plus . for each el in cdr monom collect FindRatPow1(el, var)
- else FindRatPow1(monom, var);
- symbolic procedure FindRatPow1(monom, var);
- if monom eq 1 then 0 % only possible constant by linearity
- else if monom = var then 1
- else if eqcar(monom, 'expt) and
- cadr monom = var then caddr monom
- else error1(); % wrong form
- endmodule;
- load_package int$
- % die folgende Aenderung verhindert das Erzeugen von int* ...
- remd('simpint!*)$
- symbolic procedure simpint!* u$
- begin scalar x;
- return if (x := opmtch('int . u)) then simp x
- else simpiden('int . u)
- % statt else simpiden('int!* . u)
- end$
- setcrackflags()$
- !*fullroots:=t;
- % Zum Integrieren
- put('int, 'simpfn, 'SimpIntPatch);
- share !*intflag!*; % True inside integrator.
- algebraic let {
- int(1/~x^(~n),~x) => -x/(x^n*(n-1)) when numberp n,
- ~x^(~m/~n)*~x => x**((m+n)/n) when (numberp n and numberp m),
- int(~z/~y,~x) => log(y) when z = df(y,x)};
- load_package crack2;
- end$
|