# A nonhomogeneous 2nd order ODE with constant coefficients is solved
# from 
# beginning to end using the Method of Undetermined Coefficients.
# 
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4] (Worksheet A)
# Text, Sections 3.6, 4.3: Lecture: The Method of Undetermined
# Coefficients
> 
> unprotect(LZ):

> LZ:=proc(c)

> local a, d, k, L, n, p, r, s, x, y, Z;

> n:=nops(c) - 1;

> a:=array(0..n, c);

> p:=0; s:=0;

> for k from 0 to n do

>   p:=p + a[k]*r^(n-k);

#       s:=s + a[k]*(D@@(n-k))(y)

> od;

> L:=unapply(s, y);

> Z:=unapply(p, r);

> [eval(L), eval(Z)]

> end:

> protect(LZ):

# 
# Ex 1-0 (S3.6, 162/15): Solution to 2nd order linear nonhomogenous IVP;
# setup.
> unassign('y'); setup:=LZ([1, -2, 1]);
> L:=setup[1]; Z:=setup[2];

                       (2)                         2
      setup := [y -> (D   )(y) - 2 D(y) + y, r -> y  - 2 y + 1]


                               (2)
                   L := y -> (D   )(y) - 2 D(y) + y


                                   2
                        Z := r -> r  - 2 r + 1

> homog_deq:=convert(L(y)(x)=0, diff);

                        / 2      \
                        |d       |     /d      \
           homog_deq := |--- y(x)| - 2 |-- y(x)| + y(x) = 0
                        |  2     |     \dx     /
                        \dx      /

> nonhomog_deq:=convert(L(y)(x)=x*exp(x)+4, diff);

                    / 2      \
                    |d       |     /d      \
    nonhomog_deq := |--- y(x)| - 2 |-- y(x)| + y(x) = x exp(x) + 4
                    |  2     |     \dx     /
                    \dx      /

> IC:={y(0)=1, D(y)(0)=1}; IVP:={nonhomog_deq} union IC;

                    IC := {y(0) = 1, D(y)(0) = 1}


IVP := {

    y(0) = 1, D(y)(0) = 1,

    / 2      \
    |d       |     /d      \
    |--- y(x)| - 2 |-- y(x)| + y(x) = x exp(x) + 4}
    |  2     |     \dx     /
    \dx      /

# Ex 1-1: Find general solution, yc, of the corresponding homogeneous
# problem.
> rt:=solve(Z(r)=0, r);

                              rt := 1, 1

> y:=unapply(c1*exp(x) + c2*x*exp(x), x);

                  y := x -> c1 exp(x) + c2 x exp(x)

> yc:=eval(y);

                  yc := x -> c1 exp(x) + c2 x exp(x)

> check:=simplify(homog_deq);

                            check := 0 = 0

# Ex 1-2: Ensure that (the terms of) g(x), the RHS in the nonhomogeneous
# DE, is of the form in Table 3.6.1 (T-159) or Table 4.3.1 (T-204).
# Otherwise, use the method of variation of parameters. It is in this
# case, the method of undetermined coefficients is the ticket!
> rhs(nonhomog_deq);

                             x exp(x) + 4

# Ex 1-3: Set up the subproblems, one for each term in g(x).
> unassign('y');
> nh1:=convert(L(y)(x)=x*exp(x), diff);

                 / 2      \
                 |d       |     /d      \
          nh1 := |--- y(x)| - 2 |-- y(x)| + y(x) = x exp(x)
                 |  2     |     \dx     /
                 \dx      /

> nh2:=convert(L(y)(x)=4, diff);

                     / 2      \
                     |d       |     /d      \
              nh2 := |--- y(x)| - 2 |-- y(x)| + y(x) = 4
                     |  2     |     \dx     /
                     \dx      /

# Ex 1-4: For each subproblem, assume the form of the particular
# solution to the corresponding nonhomogeneous subproblems. Looking at
# the Notes beneath Table 3.6.1, we see that s = 2 = #times alpha=1 is a
# root of the characteristic equation Z(r)=0.
> yp1:=unapply(x^2*(a*x+b)*exp(x), x);

                                2
                   yp1 := x -> x  (a x + b) exp(x)

> yp2:=x->c;

                            yp2 := x -> c

# Ex 1-5: Find the particular solution to each subproblem. Then sum them
# to form a particular solution to the full nonhomogeneous problem.
> simplify(subs(y(x)=yp1(x), nh1));

                 6 x a exp(x) + 2 exp(x) b = x exp(x)

> sol:=solve(identity(", x), {a, b});

                       sol := {b = 0, a = 1/6}

> yp1:=unapply(subs(sol, yp1(x)), x);

                                       3
                      yp1 := x -> 1/6 x  exp(x)

> 
> simplify(subs(y(x)=yp2(x), nh2));

                                c = 4

> sol:=solve(identity(", x), {c});

                            sol := {c = 4}

> yp2:=unapply(subs(sol, yp2(x)), x);

                               yp2 := 4

> 
> yp:=yp1+yp2;

                            yp := yp1 + 4

# Ex 1-6: Finally, obtain the general solution to the full
# nonhomogeneous problem by adding the general solution to the
# homogeneous problem and the particular solution to the full
# nonhomogeneous problem.
> y:=unapply(yc(x)+yp(x), x);

                                                 3
        y := x -> c1 exp(x) + c2 x exp(x) + 1/6 x  exp(x) + 4

> check:=simplify(nonhomog_deq);

                 check := x exp(x) + 4 = x exp(x) + 4

# Ex 1-7: Since this was an IVP, we now determine c1 and c2.
> IC; c_sols:=solve(", {c1, c2});

                      {c1 + 4 = 1, c1 + c2 = 1}


                     c_sols := {c1 = -3, c2 = 4}

> final_solution; y:=unapply(subs(c_sols, y(x)), x);

                            final_solution


                                                 3
         y := x -> -3 exp(x) + 4 x exp(x) + 1/6 x  exp(x) + 4

> final_check:=IVP;

         final_check := {x exp(x) + 4 = x exp(x) + 4, 1 = 1}

# "Victory at sea..."
>