# # A third order DE is completely solved and checked using Variation of # Parameters. # # Art Belmonte # Mon, 27/May/96 # Math 308-509 [Maple V Release 4] (Worksheet B) # Sections 3.7, 4.4: The Method of Variation of Parameters # # ("The government taught me, and they taught me well..." # --- The Smuggler's Blues) [P.S.: Thanks, Maury!] # # M-103/9a: # # 0. Make sure the DE is in standard form (coefficient of highest order # derivative must be 1). # # 1. If necessary, construct a fundamental set of solutions {y1, ..., # yn} to the corresponding homogeneous differential equation. > with(linalg): > unassign('y'); setup:=LZ([1, -2, -1, 2]); setup := [ (3) (2) 3 2 y -> (D )(y) - 2 (D )(y) - D(y) + 2 y, r -> y - 2 y - y + 2 ] > L:=setup[1]; Z:=setup[2]; (3) (2) L := y -> (D )(y) - 2 (D )(y) - D(y) + 2 y 3 2 Z := r -> r - 2 r - r + 2 > h_deq:={convert(L(y)(x)=0, diff)}; / 3 \ / 2 \ |d | |d | /d \ h_deq := {|--- y(x)| - 2 |--- y(x)| - |-- y(x)| + 2 y(x) = 0} | 3 | | 2 | \dx / \dx / \dx / > n_deq:={convert(L(y)(x)=exp(4*x), diff)}; / 3 \ / 2 \ |d | |d | /d \ n_deq := {|--- y(x)| - 2 |--- y(x)| - |-- y(x)| + 2 y(x) = exp(4 x)} | 3 | | 2 | \dx / \dx / \dx / > rt:=solve(Z(r)=0, r); rt := 1, 2, -1 > y1:=x->exp(x); y2:=x->exp(2*x); y3:=unapply(exp(-x), x); y1 := exp y2 := x -> exp(2 x) y3 := x -> exp(-x) # M-103/9b: # # 2. Next, verify that the functions y1, ..., yn are linearly # independent solutions to the corresponding homogeneous differential # equation. To do this, show that: # # a. They each give an identity when substituted into the homogeneous # DE. Also collect the functional EXPRESSIONS into a list [y1, ..., yn] # and label it hsols. # # b. Their Wronskian determinant is nonzero. Maple's linalg routine # Wronskian produces the Wronskian matrix (which we'll also need and # hence will label WW since W is reserved--Lambert's omega function), # then hitting it with the linalg routine det produces the determinant # `o` Wronski. > y1_satisfies_h_deq:=simplify(subs(y(x)=y1(x), h_deq)); y1_satisfies_h_deq := {0 = 0} > y2_satisfies_h_deq:=simplify(subs(y(x)=y2(x), h_deq)); y2_satisfies_h_deq := {0 = 0} > y3_satisfies_h_deq:=simplify(subs(y(x)=y3(x), h_deq)); y3_satisfies_h_deq := {0 = 0} > hsols:=[y1(x), y2(x), y3(x)]; hsols := [exp(x), exp(2 x), exp(-x)] > WW:=Wronskian(hsols, x); WW_det:=simplify(det(WW)); [exp(x) exp(2 x) exp(-x) ] [ ] WW := [exp(x) 2 exp(2 x) -exp(-x)] [ ] [exp(x) 4 exp(2 x) exp(-x) ] WW_det := 6 exp(2 x) # M-103/9c: # # 3. Form the n x 1 column matrix G, each of whose entries is zero # EXCEPT the last one, which is the functional expression corresponding # to g, the function on the RHS of the nonhomogeneous DE. > G:=matrix(3, 1, [0, 0, exp(4*x)]); [ 0 ] [ ] G := [ 0 ] [ ] [exp(4 x)] # M-103/9d: # # 4. Solve the matrix equation (WW)(uprimes) = G for uprimes, then # integrate each component of the result via the map command. ("Why do # them one at a time, when you can do them ALL at once?!") Convert the # result to a vector, then to a list, naming the result u (since you # integrated u' to obtain it!). > linsolve(WW, G); uprimes:=map(simplify, "); [ exp(4 x)] [- 1/2 --------] [ exp(x) ] [ ] [ exp(4 x) ] [ 1/3 -------- ] [ exp(2 x) ] [ ] [ exp(4 x) ] [ 1/6 -------- ] [ exp(-x) ] [- 1/2 exp(3 x)] [ ] uprimes := [ 1/3 exp(2 x) ] [ ] [ 1/6 exp(5 x) ] > map(int, uprimes, x); [- 1/6 exp(3 x)] [ ] [ 1/6 exp(2 x) ] [ ] [1/30 exp(5 x) ] > convert(", vector); u:=convert(", list); [- 1/6 exp(3 x), 1/6 exp(2 x), 1/30 exp(5 x)] u := [- 1/6 exp(3 x), 1/6 exp(2 x), 1/30 exp(5 x)] # M-103/9e: # # 5. To obtain a particular solution yp to the nonhomogeneous DE, take # the dot product of u and hsols. Simplify to taste... > dotprod(u, hsols); yp:=simplify("); 2 - 1/6 exp(3 x) exp(x) + 1/6 exp(2 x) + 1/30 exp(5 x) exp(-x) yp := 1/30 exp(4 x) # M-103/9f: # # 6. Finally, obtain the general solution y to the nonhomogeneous DE by # adding the particular solution yp to the general solution yc of the # homogeneous DE, collecting like terms and renaming constants as # necessary. A check shows that we have indeed produced the general # solution! > dotprod([c1, c2, c3], hsols) + yp; c1 exp(x) + c2 exp(2 x) + c3 exp(-x) + 1/30 exp(4 x) > y:=unapply(", x); y := x -> c1 exp(x) + c2 exp(2 x) + c3 exp(-x) + 1/30 exp(4 x) > check:=simplify(n_deq); check := {exp(4 x) = exp(4 x)} # VICTORY!!! >