1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
/* ************** beginning of main ******************************** */
lode2(oeq):=block([lhs,ay,ypx],tobegin(dum),
remvalue(y1,y2), swro:0, ys:f, nhstddeq(oeq), if rr#0 then
(print("first",str1,eq)), kt:0, ys:hlode2(eq),
if rr=0 then return(ys), print("the solution of the homog. eq. is",ys),
ldf0(8), ypx:nonhom(rr,y1,y2), return(ys+ypx))$
tobegin(dum):=block([],dyx2:'diff(y,x,2),dyx1:'diff(y,x),
dyt2:'diff(y,t,2),dyt1:'diff(y,t),dux2:'diff(u,x,2),dux1:'diff(u,x),
str1:"we solve",str2:"we use",str3:"the result is",str4:"type y or n",
str5:"positive,zero,or negative? type p,z,or n")$
nhstddeq(oeq):=block([], weq:expand(oeq), lhs:part(weq,1),
rhs:part(weq,2), lhs:lhs-rhs, deqcoeff(lhs), wlhs:lhs,
lhs:za*dyx2+zb*dyx1+zc*y, rhs:-ratsimp(wlhs-lhs),
weq:lhs=rhs, steq:stddeq(weq), lhs:part(steq,1), eq:lhs=0,
if lhs#wlhs then print(str1,steq))$
hlode2(eq):=block([exp,ro], aeq[kt]:eq, swsv:f,
ll1, leq:part(eq,1), tp:brchde1(leq),
if zc=0 then (ldf0( 9),svode1(rb),return(ys)),
if tp= 50 then (ldf0(10),skelp(rb,rc), if ys#f then return(ys)),
if swpt=y then (ldf0(11),gptm(rb,rc), if ys#f then return(ys)
else if rptm=t then (if freeof(y,vtr) then idvprep(1)
else idvprep(2),go(ll1))),
if tp>60 then (lv:listofvars(za), if length(lv)>1 then
(ldf0(16), if stddsp(za)=t then go(ll2)), tblspt(za,rb,rc,x,n),
if nspt>3 then (ldf0(12),if grmaps(dum)=y then go(ll1))),
ha:hipow(za*q1+zb*q2+zc*q3,x), netv:listofvars(ha),lg:length(netv),
if lg>0 then (ldf0(13),if ppowr(rb,rc,lg) then tp:70),
if tp<=70 then go(ll2), vecdeg(za,zb,zc,x), tp:brchde2(rb,rc),
if tp=110 then (ldf0(2),svconf(rb,rc),return(ys)) else
if tp=120 then (ldf0(3),svhypr(rb,rc),return(ys)) else
if tp=125 then (idvprep(2),go(ll1)),
if tp=126 then (idvprep(1),go(ll1)),
if tp=130 then return(ys), if swsv=tdv or tp=990 then go(ll3),
ll2,ldf0(5),neq:ivtrns(rb,rc,rr),if 10<tp and tp<14 then go(ll3),
if neq#f then (idvprep(1),go(ll1)), if tp=0 then go(ll4),
ll3, ldf0(6),neq:dvtrns(za,rb,rc,rr,ro,lam1),
if neq#f then (idvprep(2),go(ll1)),
ll4, nosol(oeq),return(oeq))$
deqcoeff(exp):=block([],za:coeff(exp,dyx2,1),zb:coeff(exp,dyx1,1),
zc:coeff(exp,y,1),rb:ratsimp(zb/za),rc:ratsimp(zc/za))$
brchde1(exp):=block([],
if not freeof(sin,cos,tan,cot,exp) then return(10),
if not freeof(%e,exp) then return(20),
if not freeof(log,exp) then return(30),
if not freeof(sinh,cosh,tanh,coth,exp) then return(40),
if not freeof(sn,cn,dn,pe,th,exp) then return(50),
if not freeof(sqrt(x),exp) then return(60),
return(990))$
brchde2(wb,wc):=block([ws,wi],
if nspt=1 then (if rk[1]=1 then return(110)),
if nspt=2 then (if rk[1]<=0 and denom(rk[2])<3 and rk[2]<=1 then
(if spt[1]#0 then (axb:x-spt[1],return(91)) else return(110))),
if nspt=3 and rk[1]=0 and rk[2]=0 then (if rk[3]=0 then
(if nf=1 and maxacf=2 then (ldf0(14),if powr(wb,wc,2) then return(80)),
if da=2 and db=1 and dc=0 then
(ldf0(15),tp:eqpt1(dum), if tp#f then return(tp)),
if da=3 and db=2 and dc=1 then
(ldf0(18),tp:eqpt12(wb,wc),if tp#f then return(tp))),
if da=4 and db=3 and dc<=4 then
(ldf0(17),tp:eqpt11(dum),if tp#f then return(tp)) else
if rk[3]=0 then return(120)),
ldf0(7),tp:brchde3(wb,wc),return(tp))$
tblspt(za,wb,wc,x,ksw):=block([k,n,wa,wf,wfi], k:0, n:hipow(za,x),
if n=0 then go(lf), pspt:1, maxacf:0, wf:facza(za), j:0,
for i:1 thru nf do (pwf[i]:part(wf,i),if pwf[i]=x then j:i),
if nf>1 and j>1 then (ws:pwf[1],pwf[1]:pwf[j],pwf[j]:ws),
for i:1 thru nf do (wfi:solve(pwf[i],x),l1:length(wfi),
lacf[i]:l1, l2:polord(za,pwf[i]), maxacf:max(maxacf,l2),
for j:1 thru l1 do (k:k+1,spt[k]:part(wfi[j],2),
if ksw#n then kindspt(wb,wc,k),
if l1=1 then pspt:pspt*spt[k] else pspt:0)),
lf,k:k+1,nspt:k,spt[k]:inf,if ksw#n then (kindspt(wb,wc,k),ro:rk[k]))$
facza(za):=block([w1,w2,w3,w4], w1:za,w2:diff(w1,x),w3:gcd(w1,w2),
w4:ratsimp(w1/w3), wf:factor(w4), if wf=w4 then nf:1 else nf:length(wf),
if nf=1 then wf:[wf], return(wf))$
vecdeg(za,zb,zc,x):=block([],da:hipow(za,x),db:hipow(zb,x),dc:hipow(zc,x))$
nosol(oeq):=block([], print("we cannot solve"), print(oeq))$
idvprep(i):=block([],print(str2,vtr),print(str3,neq),
if i=1 then eq:subst(x,t,neq) else eq:subst(y,u,neq), vprep(dum))$
vprep(dum):=block([],eq:stddeq(eq),print(str1,eq),kt:kt+1,aeq[kt]:eq,
avtr[kt]:vtr)$
stddeq(oeq):=block([wexp,wr,deq],wexp:part(oeq,1),
wr:part(oeq,2),deqcoeff(wexp),rr:ratsimp(wr/za),
za:calca(rb,rc,x), zb:ratsimp(rb*za), zc:ratsimp(rc*za),
deq:dyx2+rb*dyx1+rc*y=rr, return(deq))$
calca(wb,wc,x):=block([],w1:denom(wb),w2:denom(wc),za:lcm(w1,w2,x))$
lcm(p,q,x):=block([r,s],r:gcd(p,q,x),s:ratsimp(p*q/r),return(s))$
polord(exp,sexp):=block([we,n],we:exp,n:0,ll,we:ratsimp(we/sexp),
if denom(we)#1 then return(n), n:n+1,go(ll))$
kindspt(wb,wc,k):=block([bt,ct,b1,c1,dbt,dct],
if spt[k]=inf then alpha:0 else alpha:spt[k],
bt:ratsubst(1/t+alpha,x,wb), b1:ratsimp((2*t-bt)/t^2),
ct:ratsubst(1/t+alpha,x,wc), c1:ratsimp(ct/t^4),if spt[k]=inf then
(bt:subst(x,t,b1),ct:subst(x,t,c1),dbt:ratdeg(wb,x),dct:ratdeg(wc,x))
else (bt:ratsimp(wb),ct:ratsimp(wc), dbt:ratdeg(b1,t),dct:ratdeg(c1,t)),
rk[k]:max(dbt,dct/2)+1, if rk[k]=0 and freeof(%i,alpha)
then regsingpt(x-alpha,bt,ct,alpha))$
ratdeg(rf,x):=block([wrf,nrf,dgnrf,drf,dgdrf,dgrf],if rf=0 then
return(-99), wrf:ratsimp(rf), nrf:num(wrf), dgnrf:hipow(nrf,x),
drf:denom(wrf),dgdrf:hipow(drf,x), dgrf:ratsimp(dgnrf-dgdrf),
return(dgrf))$
regsingpt(xa,wb,wc,alpha):=block([f1,g1,flx], remvalue(l),
f1:ratsimp(xa*wb), g1:ratsimp(xa^2*wc), flx:l^2+(f1-1)*l+g1,
cheq[k]:subst(alpha,x,flx),cflx[k]:flx)$
reqans1(dum):=block([],ans:read("continue?",str4))$
ldf0(i):=block([],if ldsw[0,i]#y then (ldsw[0,i]:y,
if i= 2 then load(pconfl) else
if i= 3 then load(phypgm) else
if i= 5 then load(pivtr) else
if i= 6 then load(pdvtr) else
if i= 7 then load(pbrnch) else
if i= 8 then load(pnonh) else
if i= 9 then load(psode1) else
if i=10 then load(pskelp) else
if i=11 then load(pgptm) else
if i=12 then load(premap) else
if i=13 then load(ppow) else
if i=14 then load(pnpow) else
if i=15 then load(peqp1) else
if i=16 then load(pstdsp) else
if i=17 then load(peqp11) else
if i=18 then load(peqp12) ))$
lodesave(filename,[lode_functions]):=
if status(feature,its)=true
then apply('fassave,cons(filename,lode_functions))
else block([packagefile:true],
if status(feature,unix)=true
then apply('save,
cons(concat(first(filename),".l"),lode_functions))
else error("unknown system"))$
lodesave([pmain,fasl],lode2,tobegin,nhstddeq,hlode2,brchde1,vecdeg,
brchde2,deqcoeff,tblspt,facza,nosol,idvprep,vprep,stddeq,calca,
lcm,polord,kindspt,ratdeg,regsingpt,reqans1,ldf0);
lodesave([lodsav,fasl],lodesave)$
|
