> 
# A fourth order DE is completely solved and checked using 
# Variation of Parameters.
--------------------------------------------------------------------------------
# 
# Art Belmonte
# Thu, 15/Feb/96
# Math 308-509 (Worksheet C)
# 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):
--------------------------------------------------------------------------------
# a. Homebrew: a 4th order problem (T-205/5 with a different right 
# hand side, namely 1/x.)
# 
# 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, 0, -4, 0, 0]);\
L:=setup[1]; Z:=setup[2];\
h_deq:={convert(L(y)(x)=0, diff)};\
n_deq:={convert(L(y)(x)=1/x, diff)};\
rt:=solve(Z(r)=0, r);\
y1:=x->1; y2:=x->x; y3:=unapply(exp(2*x), x); y4:=unapply(exp(-2*x), x);
--------------------------------------------------------------------------------
# b:
# 
# 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));\
y4_satisfies_h_deq:=simplify(subs(y(x)=y4(x), h_deq));\
hsols:=[y1(x), y2(x), y3(x), y4(x)];\
WW:=Wronskian(hsols, x); WW_det:=simplify(det(WW));
--------------------------------------------------------------------------------
# c:
# 
# 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(4, 1, [0, 0, 0, 1/x]);
--------------------------------------------------------------------------------
# d:
# 
# 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!).
> uprimes:=linsolve(WW, G);\
map(int, uprimes, x);\
convert(", vector); u:=convert(", list);
--------------------------------------------------------------------------------
# e:
# 
# 5. To obtain a particular solution yp to the nonhomogeneous DE, 
# take the dot product of u and hsols. Simplify to taste... (NOTE: 
# Type ?Ei to learn about the Exponential Integral.)
> dotprod(u, hsols); yp:=simplify(");
--------------------------------------------------------------------------------
# f:
# 
# 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, c4], hsols) + yp;\
y:=unapply(", x);\
check:=simplify(n_deq);
--------------------------------------------------------------------------------
# VICTORY!!!
>