> 
# A third order DE is completely solved and checked using 
# Variation of Parameters.
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Thu, 15/Feb/96
# Math 308-509 (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!]
# 
> with(linalg):
--------------------------------------------------------------------------------
# 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.
> unassign('y'); setup:=LZ([1, -2, -1, 2]);\
L:=setup[1]; Z:=setup[2];\
h_deq:={convert(L(y)(x)=0, diff)};\
n_deq:={convert(L(y)(x)=exp(4*x), diff)};\
rt:=solve(Z(r)=0, r);\
y1:=x->exp(x); y2:=x->exp(2*x); y3:=unapply(exp(-x), 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));\
y2_satisfies_h_deq:=simplify(subs(y(x)=y2(x), h_deq));\
y3_satisfies_h_deq:=simplify(subs(y(x)=y3(x), h_deq));\
hsols:=[y1(x), y2(x), y3(x)];\
WW:=Wronskian(hsols, x); WW_det:=simplify(det(WW));
--------------------------------------------------------------------------------
# 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)]);
--------------------------------------------------------------------------------
# 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, ");\
map(int, uprimes, x);\
convert(", vector); u:=convert(", list);
--------------------------------------------------------------------------------
# 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(");
--------------------------------------------------------------------------------
# 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;\
y:=unapply(", x);\
check:=simplify(n_deq);
--------------------------------------------------------------------------------
# VICTORY!!!
>