share/contrib/odes/odes.mac

   1
   2 /*
   3 odes  package for solving ordinary diferential equations
   4 version 2.01,  2013.11
   5 Copyright (C)  A.Domarkas 2013
   6 odes package is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License,
   8 */
   9
  10 /*
  11 Examples see in odes-doc.pdf
  12 */
  13
  14 /*
  15 1
  16 */
  17
  18 odecv(tr,eq,y,x):=block([_n,_t,tr1,itr,eq1,r,solvetrigwarn],
  19 solvetrigwarn:false,_n:derivdegree(eq,y,x),
  20 tr1:(solve(tr,x))[1],
  21 _t:last(listofvars(rhs(tr1))),
  22 itr:(solve(tr,_t))[1],
  23 depends(y,_t,_t,x),
  24 eq1:ev(eq,nouns),
  25 makelist(diff(itr,x,k),k,1,_n),
  26 subst(%%,eq1),
  27 ev(%%,nouns),
  28 r:subst(tr1,%%),
  29 remove([y,_t],dependency),
  30 radcan(r))$
  31
  32 /*
  33 2
  34 */
  35
  36 dchange(tr,eq,fun,var,nfun,nvar):=block(local(T,keit),
  37 diffpf(f,t,k):=block(
  38 if k = 0 then f[2] elseif k = 1
  39 then ratsimp(diff(f[2],t,1)/diff(f[1],t,1))
  40 else ratsimp(diff(diffpf(f,t,k-1),t,1)/diff(f[1],t,1))),
  41 T:makelist('diff(fun,var,k)=
  42 diffpf([rhs(tr),nfun],nvar,k),k,0,derivdegree(eq,fun,var)),
  43 keit:append(reverse(T),[tr]),
  44 ev(eq,keit))$
  45
  46 /*
  47 3
  48 */
  49
  50 odeC(eq,r,x):=block([derivsubst,_u,tr,itr,L],
  51 depends(_u,x),assume(_u>0),
  52 derivsubst:true,
  53 L:listofvars(r),
  54 tr:_u=r, itr:solve(tr,last(L))[1],
  55 subst(itr,eq),
  56 ev(%%, nouns), expand(%%),
  57 ode2(%%,_u,x),
  58 subst(tr,expand(%%))
  59 )$
  60
  61 /*
  62 4
  63 */
  64
  65 solvet(eq,x):=block([spr,k,maperror,mapprint],
  66 maperror:false,mapprint:false,
  67 spr:solve(eq,x),
  68 rectform(%%),
  69 if freeof(sin,%%) then return(sort(%%))
  70 else makelist(x=map(polarform,rhs(spr[k])),k,1,length(spr)),
  71 rectform(%%),
  72 trigsimp(%%),
  73 sort(%%))$
  74
  75 /*
  76 5
  77 */
  78
  79 ode1_ic(eq,y,x,_ic):=block([sol,gsol, spr,ats,p],
  80 solvetrigwarn:false,
  81 sol:ode2(eq,y,x),
  82 if %%=false then return(false),
  83 trigsimp(sol),
  84 radcan(%%),
  85 trigsimp(%%),
  86 trigreduce(%%),
  87 trigexpand(%%),
  88 gsol:rootscontract(%%),
  89 subst([x=_ic[1],y=_ic[2]],trigsimp(gsol)),
  90 solve(%%,%c),
  91 subst(%%,gsol),
  92 sol:ratsimp(logcontract(%%)),
  93 spr:solve(%%,y),
  94 if not freeof(y,map(rhs,%%)) then return(sol),
  95 if length(spr)=1 and freeof(y,rhs(spr[1])) then return(spr[1]),
  96 sublist(spr,lambda([e],is(radcan(subst([x=_ic[1], y=_ic[2]],e))))),
  97 trigsimp(%%[1]),
  98 radcan(%%),
  99 if not freeof(%e,%%) then
 100 radcan(%%) else %%,
 101 if not freeof(%i,%%) then
 102 ats:rootscontract(subst(%i=sqrt(-1),%%))
 103 else ats:%%,
 104 return(ats)
 105 )$
 106
 107 /*
 108 6
 109 */
 110
 111 ode2_ic(eq,y,x,_ic):=block([p,sol,ats,spr],
 112 solvetrigwarn:false,
 113 maperror:false,
 114 mapprint:false,
 115 assume(x>_ic[1]),
 116 sol:ode2(eq,y,x),
 117 solve(eq,'diff(y,x,2)),
 118 rhs(%%[1]),
 119 if listofvars(%%)=[y] and listp(sol) then
 120          (
 121           subst('diff(y,x,2)=p*'diff(p,y),eq),
 122           ode1_ic(%%,p,y,[_ic[2],_ic[3]]),
 123           subst(p='diff(y,x),%%),
 124           ats:ode1_ic(%%,y,x,[_ic[1],_ic[2]])
 125          )
 126 else
 127         (
 128          ic2(sol, x=_ic[1], y=_ic[2],'diff(y,x)=_ic[3]),
 129          ats:lhs(%%)=map(factor,rhs(%%)),
 130          if lhs(%%)#y then
 131               (
 132                spr:solve(ats,y),
 133                if spr#[] then
 134                     (
 135                      sublist(spr,lambda([e],
 136                      _ic[2]=ratsimp(ev(rhs(e),x=_ic[1])) and
 137                      _ic[3]=ratsimp(at(diff(rhs(e),x),x=_ic[1])))),
 138                      ats:%%[1]
 139                     )
 140               ),
 141         forget(facts(x)),
 142         return(ats)
 143         )
 144 )$
 145
 146 /*
 147 6a
 148 */
 149
 150 ode_ic(eq,y,x,_ic):=block([],
 151 if  derivdegree(eq,y,x)=1 then ode1_ic(eq,y,x,_ic)
 152 elseif derivdegree(eq,y,x)=2 then ode2_ic(eq,y,x,_ic)
 153 else error("equation must be the first or second order"),
 154 rootscontract(%%)
 155 )$
 156
 157 /*
 158 7
 159 */
 160
 161 P_iter(eq,x,y,x0,y0,n):=block([f,t],
 162 solve(eq,'diff(y,x,1))[1],
 163 define(f(x,y),rhs(%%)),
 164 if n=0 then y0 else
 165 (assume(x>x0),
 166 subst(x=t,f(t,P_iter(eq,x,y,x0,y0,n-1))),
 167 y0+expand(integrate(%%,t,x0,x))))$
 168
 169 /*
 170 8
 171 */
 172
 173 ode1taylor(eq,x0,y0,n):=block([spr,nezin],
 174 atvalue(y(x),x=x0,y0),
 175 spr:taylor(y(x),x,x0,n),
 176 nezin:at(makelist(diff(y(x),x,k),k,1,n),x=x0),
 177 makelist(diff(eq,x,k),k,0,n-1),
 178 at(%%,x=x0),solve(%%,nezin),
 179 ratexpand(ev(spr,first(%%))),
 180 taylor(%%,x,x0,n)
 181 )$
 182
 183 /*
 184 9
 185 */
 186
 187 ode2taylor(eq,x0,y0,y1,n):=block([spr,nezin],
 188 atvalue(y(x),x=x0,y0),
 189 atvalue(diff(y(x),x),x=x0,y1),
 190 spr:taylor(y(x),x,x0,n),
 191 nezin:at(makelist(diff(y(x),x,k),k,2,n),x=x0),
 192 makelist(diff(eq,x,k),k,0,n-2),
 193 at(%%,x=x0),solve(%%,nezin),
 194 ratexpand(ev(spr,first(%%))),
 195 taylor(%%,x,x0,n)
 196 )$
 197
 198 /*
 199 10
 200 */
 201
 202 ode1exact(eq):=block([P,Q,F],
 203 P:coeff(expand(lhs(eq)-rhs(eq)),dx),
 204 Q:coeff(expand(lhs(eq)-rhs(eq)),dy),
 205 if trigsimp(trigexpand(diff(P,y)-diff(Q,x)))=0 then
 206 (F:integrate(P,x)+g(y),diff(F,y)=Q,
 207 ode2(%%,g(y),y),
 208 subst([%%,%c=0],F)=C,
 209 trigexpand(%%),
 210 trigsimp(%%),
 211 expand(%%))
 212 else false)$
 213
 214 /*
 215 11
 216 */
 217
 218 intfactor1(eq,omega):=block([P,Q,d,t,ans],
 219 P:coeff(expand(lhs(eq)-rhs(eq)),dx),
 220 Q:coeff(expand(lhs(eq)-rhs(eq)),dy),
 221 (diff(P,y)-diff(Q,x))/(Q*diff(omega,x)-P*diff(omega,y)),
 222 ratsimp(ev(%%,diff)),
 223 trigexpand(%%),
 224 ratsubst(t,omega,%%),
 225 d:map(radcan,expand(%%)),
 226 if listofvars(%%)=[t] or listofvars(%%)=[] then (
 227 exp(integrate(d,t)),
 228 ans:subst(t=omega,%%))
 229 else return(false),
 230 if freeof(integrate, ans) then return(ans)
 231 else return(false)
 232 )$
 233
 234 /*
 235 12
 236 */
 237
 238 odeL(eq,y,x):=block([n,_yp,gamma_expand,i,_Y,_C],
 239 gamma_expand:true,
 240 _Y:fs(eq,y,x),
 241 n:derivdegree(eq,y,x),
 242 _yp:partsol(eq,y,x),
 243 _C:makelist(concat(C,i),i,1,n),
 244 y=_C._Y+_yp)$
 245
 246 /*
 247 13
 248 */
 249
 250 odeL_ic(eq,y,x,ic):=block(odeL(eq,y,x), icn(%%,y,x,ic))$
 251
 252 icn(sol,y,x,ic):=block([_spr,n],
 253 _spr:rhs(sol),
 254 n:length(ic)-1,
 255 makelist(at(diff(_spr,x,k),x=ic[1])=ic[k+2],k,0,n-1),
 256 solve(%%,_C:makelist(concat(C,i),i,1,n)),
 257 y=subst(%%,_spr)
 258 )$
 259
 260 /*
 261 14
 262 */
 263
 264 fs(eq,y,x):=block([_f,n,j,cr,_k,ce],
 265 _f(r):=block([a,b,i],
 266        a:radcan(realpart(r[1])),
 267        b:radcan(imagpart(r[1])),
 268        if b=0 then
 269        return(makelist(x^(i-1)*exp(a*x),i,1,r[2]))
 270        else
 271        join(makelist(x^(i-1)*exp(a*x)*cos(b*x),i,1,r[2]),
 272        makelist(x^(i-1)*exp(a*x)*sin(b*x),i,1,r[2]))),
 273 n:derivdegree(eq,y,x),
 274 makelist('diff(y,x,j)=_k^j,j,0,n),
 275 reverse(%%),
 276 ce:subst(%%,lhs(eq))=0,
 277 solve(ce,_k),
 278 map(rhs,%%),
 279 map(rootssimp,%%),
 280 makelist([%%[j],multiplicities[j]],j,1,length(%%)),
 281 cr:sublist(%%,lambda([e],is(radcan(imagpart(e[1]))>=0))),
 282 map(_f,cr),
 283 flatten(%%),
 284 trigreduce(%%),
 285 sort(%%))$
 286
 287 /*
 288 15
 289 */
 290
 291 partsolL_UC(eq,y,x):=block([spr,u],
 292 if rhs(eq)=0 then return(0),
 293 subst(y=u(x),eq),
 294 rhs(desolve([%%],[u(x)])),
 295 args(%%),
 296 sublist(%%,lambda([e],freeof(u,ilt,e))),
 297 spr:apply("+",%%),
 298 if not atom(spr) and op(spr)="+" then
 299 (args(expand(spr)),
 300 sublist(%%,lambda([e],ode_check(lhs(eq),y=e)#0) ),
 301 apply("+",%%),return(ratsimp(%%)))
 302 else
 303 return(ratsimp(spr))
 304 )$
 305
 306 partsolL_VP(eq,y,x):=block([W,Wi,_n,_dp],
 307 if rhs(eq)=0 then return(0),
 308 FS:fs(eq,y,x),
 309 F:rhs(eq),
 310 W:wronskian(FS,x),
 311 _n:length(FS),
 312 makelist(0,i,1,_n-1),
 313 _dp:append(%%,[F]),
 314 Wi:radcan(trigsimp(invert(W))),
 315 integrate(Wi._dp,x),
 316 ratsimp(%%),
 317 trigsimp(%%),
 318 logcontract(%%),
 319 FS.%%,
 320 ratsimp(%%),
 321 trigsimp(%%),
 322 trigreduce(%%),
 323 ratsimp(%%)
 324 )$
 325
 326 partsol(eq,y,x):=block([s,eq1],
 327 eq1:subst(y=y(x),eq),
 328 if freeof(laplace,laplace(rhs(eq),x,s))
 329 and freeof(log,tan,atan,sec,csc,rhs(eq))
 330 then partsolL_UC(eq,y,x)
 331 else partsolL_VP(eq,y,x)
 332 )$
 333
 334 /*
 335 16
 336 */
 337
 338 odeM(A,F,t):=block([W,gamma_expand,n,Ya],
 339 gamma_expand:true,
 340 n:length(F),
 341 W:matrix_exp(A,t),
 342 W.integrate(trigsimp(invert(W).F),t),
 343 trigrat(expand(%%)),
 344 Ya:expand(%%),
 345 W.transpose(makelist(concat(C,i),i,1,n))+Ya)$
 346
 347 /*
 348 17
 349 */
 350
 351 odeM_ic(A,B,t,t0,Y0):=block([bspr],
 352 bspr:odeM(A,B,t),
 353 subst(t=t0,%%)-Y0,
 354 list_matrix_entries(%%),
 355 solve(%%,makelist(concat(C,i),i,1,length(Y0))),
 356 subst(%%[1],bspr),
 357 expand(%%)
 358 )$
 359
 360 /*
 361 18
 362 */
 363
 364 matrix_exp(A,r):=block([n,B,s,t,Lap,f],
 365 n:length(A),
 366 B:invert(s*ident(n)-A),
 367 Lap(f):=ilt(f, s, t),
 368 matrixmap(Lap,B),
 369 subst(t=r,%%))$
 370
 371 /*
 372 19
 373 */
 374
 375 odelinsys(A,F,t,t0,Y0):=block([s,g,expA],
 376 assume(t>t0),
 377 define(g(t),F),
 378 mat_function(exp,A*t),ratsimp(%%),
 379 define(Aexp(t),%%),
 380 Aexp(t-t0).Y0+integrate(Aexp(t-s).g(s),s,t0,t),
 381 expand(%%)
 382 )$
 383
 384 /*
 385 20
 386 */
 387
 388 wronskian(functlist,var):=block([end],end:length(functlist)-1,functlist:[
 389 functlist],thru end do functlist:
 390 endcons(map(lambda([x],diff(x,var)),last(functlist)),functlist),
 391 apply('matrix,functlist))$
 392
 393 /*
 394 etc
 395 */
 396
 397 itr2can(eq,y,x):=block(coeff(lhs(eq),'diff(y,x))/coeff(lhs(eq),'diff(y,x,2)),
 398 t=radcan(integrate(exp(-integrate(%%,x)),x)))$
 399
 400 rootssimp(x):=block([maperror,mapprint],
 401 maperror:false,
 402 mapprint:false,
 403 map(polarform,x),rectform(%%),trigsimp(%%))$
 404
 405 trigvalue(r):=block(
 406 [f,x,sol,spr,spr1,solvetrigwarn,algebraic],
 407 solvetrigwarn:false,
 408 algebraic:true,
 409 if freeof(%pi,r) then return(r),
 410 f:part(r,0),
 411 if part(r,0)="-" then f:part(r,1,0),
 412 if f=cot then f:tan,
 413 sol:solve(x=r,%pi)[1],
 414 sol*denom(rhs(sol)),
 415 map(f,%%),
 416 trigexpand(%%),
 417 factor(%%),
 418 spr:solve(%%,x),
 419 if (length(spr)<=2 or not freeof(%i,%%)) then return(r),
 420 spr1:sublist(spr,lambda([e],is(abs(rhs(e)-r)<ratepsilon))),
 421 if %%=[] then return(r),
 422 rhs(spr1[1]),
 423 sqrtdenest(%%),
 424 factor(%%)
 425 )$
 426