share/diffequations/linde1.mc

   1 /*  THIS LITTLE PACKAGE SOLVES FIRST ORDER LINEAR
   2      ORDINARY DIFFERENTIAL EQUATIONS SUBJECT TO A
   3      BOUNDARY CONDITION (B.C.) AT AN INITIAL POINT
   4        THE CALLING PROCEDURE IS
   5           IVPSOL(DIFFEQ,Y,X,A,BCEQ);
   6      WHERE DIFFEQ IS THE DIFFERENTIAL EQUATION TO BE SOLVED
   7             (E.G.:  X*'DIFF(Y,X)+2*Y=1)
   8            Y IS THE DEPENDENT VARIABLE
   9            X IS THE INDEPENDENT VARIABLE
  10            A IS THE POINT AT WHICH THE B.C. IS APPLIED
  11            BCEQ IS THE B.C. EQUATION (E.G.: Y=2)
  12 */
  13 soldiff(eq1,y,x,a):=
  14      block([eq2,a1,a2,a3,a4,a5],
  15         eq2: lhs(eq1)-rhs(eq1),
  16         a1:  ratcoef(eq2,'diff(y,x)),
  17         a2:  ratcoef(eq2,y),
  18         a3:  eq2-a1*'diff(y,x)-a2*y,
  19         a4:  integrate(a2/a1,x),
  20         a5:  integrate(%e**a4*a3/a1,x),
  21            'const1/%e**a4-a5/%e**a4)$
  22 ivpsol(diffeq,y,x,a,bceq):=
  23      block([z,dery,dya,ya,co1],
  24         z:   soldiff(diffeq,y,x,a),
  25         dery:diff(z,x,1),
  26         dya: subst(a,x,dery),
  27         ya:  subst(a,x,z),
  28         co1: solve(subst(ya,y,subst(dya,'diff(y,x),bceq)),'const1),
  29         if length(co1)>1 /* LISTP(CO1) */ then
  30            [print("SOLUTIONS ARE NOT UNIQUE, POSSIBLE SOLUTIONS ARE:
  31              "),
  32             map(lambda([l1],apply('display,[l1])),co1), return(z)]
  33            else   if co1=[] then return(print("No solution for given initial
  34 condition"),z) else subst(rhs(co1[1]),'const1,z))$
  35 /* FIND RIGHT EIGENVECTORS WITH OTHER COMPONENT = 1
  36           I.E. X IS SOLUTION OF M*X=MU*X    */
  37 evec(m,mu,modes):=block([equations,unknowns,x],
  38 (for l:1 thru modes do[
  39      equations:[],
  40      unknowns:[],
  41 /*     KILL(X), */
  42      x[l]:1,
  43      for i:1 thru modes do if not (i=l) then [
  44          unknowns:cons(x[i],unknowns),
  45          equations:cons(sum(m[i,j]*x[j],j,1,modes)
  46               =mu[l]*x[i],equations)],
  47      print(solve(equations,unknowns))]))$
  48 /* CONSTRUCT 3 X 3 MATRIX WITH DESIRED EIGENVALUES */
  49
  50 consmat3(l1,l2,l3):=block([a,b,mat3,cp,cpm,sol1,l],
  51 mat3:matrix([1,2,3],[b,3,a],[1,1,l1+l2+l3-4]),
  52 cp:expand(-(l-l1)*(l-l2)*(l-l3)),
  53 cpm:charpoly(mat3,l),globalsolve:true,
  54 sol1:solve([ratcoef(cp,l)=ratcoef(cpm,l),
  55 ratcoef(cp,l,0)=ratcoef(cpm,l,0)],[a,b]),
  56 return(apply('ev,[mat3])))$