share/contrib/rand/Ima

   1 Suggested MACSYMA Homework No.3,    Richard Rand
   2
   3 This homework involves using perturbations to treat the nonlinear
   4 boundary value problem (in the small e limit) :
   5
   6 (1)               v'' - e x v' (1+v) = 1     on  0 < x < 1
   7
   8 with the boundary conditions:
   9
  10 (2)               v(0) = v(1) = 0
  11
  12 1. Solve to order 2 in e by expanding v in a truncated power series in
  13 e:
  14
  15                   v = f[0](x) + f[1](x) e + f[2](x) e^2
  16
  17 Use MACSYMA in the interactive (non-programming) mode to find equations
  18 on the subscripted variables f[0],f[1] and f[2].  Use DEPENDS(F,X).
  19 Then solve these equations recursively using ODE2.  Evaluate the
  20 arbitrary constants %K1 and %K2 at each step by using the boundary
  21 conditions (2).
  22
  23 Plot your resulting expression for v(x) for e=0 and e=10 on the same
  24 graph.
  25
  26 2.  Write a program in the form of a BATCH file to accomplish the above
  27 task to arbitrary order n in e.
  28
  29 The program should:
  30 i)   Use READ to obtain the truncation order n from the keyboard.
  31 ii)  Expand v in a truncated power series:
  32
  33                    v = f[0] + f[1] e + ... + f[n] e^n
  34
  35 Use subscripted variables f[i] to SUM the series and DEPENDS(F,X).
  36 iii) Substitute v into (1) and TAYLOR(...,e,0,n) in order to facilitate
  37 collecting terms.  Use a FOR I:0 THRU N DO loop to store the
  38 coefficient of e^i in EQ[I].
  39 iv)  Prepare a list (now empty) to hold your intermediate results,
  40 i.e., your soon-to-be-derived values for f[0], f[1], ..., by writing:
  41      RESULTS:[];
  42 v)   Now set up another DO loop, the main loop, which does the
  43 following for each step of the recursive process:
  44      a) Plug RESULTS into EQ[I]. (Of course this will be lame for I=0,
  45 the first time thru.)
  46      b) Use ODE2 to solve EQ[I] for f[i].
  47      c) Use the boundary conditions (2) to evaluate the arbitrary
  48 constants %K1 and %K2.
  49      d) APPEND your value of f[i] to the list of RESULTS.
  50 vi)  Finally plug RESULTS into your expansion for the variable v (cf.
  51 step ii).
  52
  53 Run your program using BATCH("/ima/yourname/filename.filetype");
  54 Enter n=2 and cf. with your previously obtained expression for v (in
  55 part 1.)  Run it again with n=4 and plot n=2 and n=4 for e=10 on the
  56 same graph.
  57 -----------------------------------------------------------
  58
  59 /* nonlinear bvp        */
  60 /* ima hw3 28 june 89   */
  61
  62 n:read("enter truncation order");
  63 depends(f,x);
  64 v:sum(f[i]*e^i,i,0,n);
  65 de:diff(v,x,2)-e*x*diff(v,x)*(1+v)=1;
  66 de1:taylor(de,e,0,n);
  67 for i:0 thru n do eq[i]:coeff(de1,e,i);
  68 results:[];
  69 /* main loop */
  70 for i:0 thru n do (
  71    eq1[i]:ev(eq[i],results,diff),
  72    sol:ode2(eq1[i],f[i],x),
  73    rhs:rhs(sol),
  74    bc:[ev(rhs,x=0),ev(rhs,x=1)],
  75    const:solve(bc,[%k1,%k2]),
  76    sol1:ev(sol,const),
  77    print(sol1),
  78    results:append(results,[sol1]));
  79 v1:ev(v,results);
  80 -----------------------------------------------------------
  81 Suggested MACSYMA Homework No.4        R.Rand
  82
  83 This homework investigates the use of computer algebra to
  84 facilitate a coordinate transformation (i.e. a change of
  85 dependent variables) in a system of ode's.
  86
  87 1. Given a system of 2 autonomous ode's,
  88
  89 (1)       x1' = f1(x1,x2),   x2' = f2(x1,x2)
  90
  91 and a transformation of variables
  92
  93 (2)       x1 = g1(y1,y2),    x2 = g2(y1,y2)
  94
  95 it is desired to find the resulting system of 2 ode's
  96 on y1 and y2.
  97
  98 Write a program to accomplish this.  Your program should
  99 consist of a single function, TRANSFORM(), which should be
 100 written outside of MACSYMA using an editor and BATCHed in.
 101 It should READ the functions fi and gi from the keyboard.
 102
 103 Hint:  Plug (2) into (1) and SOLVE for [diff(y1,t),diff(y2,t)].
 104
 105 As a check on your program, try it on:
 106
 107 (3)       u' = -v + u^3 + u v^2,   v' = u + u^2 v + v^3
 108
 109 with the change of variables to polar coordinates:
 110
 111 (4)       u = r cos h,          v = r sin h
 112
 113 which should yield the result:
 114
 115 (5)       r' = r^3,             h' = 1
 116
 117 2. Now you will apply your program to the problem of determining
 118 the stabilty of the origin in the system:
 119
 120 (6)       x' = -y + x^2 y,      y' = x + x^2 y
 121
 122 Use your program to perform the near-identity transformation:
 123
 124           x = u + a30 u^3 + a21 u^2 v + a12 u v^2 + a03 v^3
 125 (7)
 126           y = v + b30 u^3 + b21 u^2 v + b12 u v^2 + b03 v^3
 127
 128 where aij and bij are as yet undetermined constants.  After obtaining
 129 the transformed ode's, eliminate terms of degree 4 and higher by using
 130
 131           TAYLOR( ... ,[u,v],0,3)
 132
 133 Next use your program again to transform to polar coordinates (4).
 134 Accomplish trigonometric simplification by using TRIGSIMP then
 135 TRIGREDUCE and finally EXPAND.
 136 Then select the coefficients aij and bij so that the resulting ode's
 137 are in the form:
 138
 139 (8)        r' = k1 r^3 + ...,   h' = 1 + k2 r^2 + ...
 140
 141 i.e., by requiring the coefficients of sin(2h),cos(2h),sin(4h),cos(4h)
 142 to vanish.
 143
 144 The sign of k1 will determine the stability of the origin.
 145
 146 --------------------------------------------------------------------
 147 /* hw4 ima 29 june 89 */
 148 /* transform 2 ode's  */
 149
 150 transform():=(
 151 x1:read("enter symbol for original variable 1"),
 152 x2:read("enter symbol for original variable 2"),
 153 y1:read("enter symbol for new variable 1"),
 154 y2:read("enter symbol for new variable 2"),
 155 print("the original eqs are of the form ",
 156 x1,"= f1(",x1,",",x2,"), ",
 157 x2,"= f2(",x1,",",x2,")"),
 158 f1:read("enter f1(",x1,",",x2,")"),
 159 f2:read("enter f2(",x1,",",x2,")"),
 160 print("the transformation is of the form ",
 161 x1,"= g1(",y1,",",y2,"), ",
 162 x2,"= g2(",y1,",",y2,")"),
 163 g1:read("enter g1(",y1,",",y2,")"),
 164 g2:read("enter g2(",y1,",",y2,")"),
 165 depends([x1,x2,y1,y2],t),
 166 eqs:ev([diff(x1,t)=f1,diff(x2,t)=f2],ev([x1=g1,x2=g2]),diff),
 167 unk:[diff(y1,t),diff(y2,t)],
 168 neweqs:solve(eqs,unk))$