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