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