share/contrib/rand/Laplace.bk1

   1 /* Filename <name>.mac
   2
   3    ***************************************************************
   4    *                                                             *
   5    *                     <package name>                          *
   6    *                <functionality description>                  *
   7    *                                                             *
   8    *          from: Computer Algebra in Applied Math.            *
   9    *                   by Rand (Pitman,1984)                     *
  10    *                Programmed by Richard Rand                   *
  11    *      These files are released to the public domain          *
  12    *                                                             *
  13    ***************************************************************
  14 */ /*
  15 \r(d4) This program solves a system of N first order
  16
  17
  18 linear nonhomogeneous constant coefficient O.D.E's
  19
  20
  21 of the form         X[I]' = F(X[1],...,X[N],T)
  22
  23
  24 by using Laplace transforms.  To call it, type:
  25
  26
  27                     LAPTR()
  28
  29 */
  30 \rlaptr():=(input1(),setup1(),input2(),setup2(),step1(),step2(),step3(),step4(),
  31       step5(),output())$
  32 input1():=n:read("ENTER NUMBER OF 1ST ORDER ODE'S")$
  33 setup1():=(remarray(x),for i thru n do x[i]:x[i](t))$
  34 input2():=(
  35        print("ENTER RHS'S OF ODE'S\
  36 USE X[I] FOR ITH DEPENDENT VARIABLE\
  37 AND T FOR INDEPENDENT VARIABLE"),
  38        for i thru n do (f[i]:read("RHS OF EQ",i,":"),print(f[i])),
  39        print("ENTER INITIAL CONDITIONS"),
  40        for i thru n do z[i]:read("VALUE OF X[",i,"] AT T=0:"))$
  41 setup2():=for i thru n do eq[i]:diff(x[i],t) = f[i]$
  42 step1():=for i thru n do transf[i]:laplace(eq[i],t,s)$
  43 step2():=for i thru n do
  44        transf[i]:subst(-z[i],part(lhs(transf[i]),2),transf[i])$
  45 step3():=for i thru n do
  46        (for j thru n do transf[j]:subst(l[i],laplace(x[i],t,s),transf[j]))$
  47 step4():=lt:solve(makelist(transf[i],i,1,n),makelist(l[i],i,1,n))$
  48 step5():=(remarray(x),for i thru n do x[i]:ilt(rhs(part(lt,1,i)),s,t))$
  49 output():=for i thru n do (print("x[",i,"] =",x[i]),print(" "))$