#
# A 2nd order De is solved completely with application to a circuit problem.
#
# Art Belmonte
# Mon, 27/May/96
# Math 308-509 [Maple V Release 4]
# 
# Text, Section 3.9: Forced Vibrations
# 
# Ex 2-0 (S3.9, 190/18): Electric circuit
> 
> 
> 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):
> 
> L:=1; R:=5000; C:=1/4 / 1000000;

                                L := 1


                              R := 5000


                            C := 1/4000000

> unassign('Q'); setup:=LZ([L, R, 1/C]);

setup :=

            (2)
    [y -> (D   )(y) + 5000 D(y) + 4000000 y,

          2
    r -> y  + 5000 y + 4000000]

> LDO:=setup[1]; Z:=setup[2];

                          (2)
            LDO := y -> (D   )(y) + 5000 D(y) + 4000000 y


                              2
                   Z := r -> r  + 5000 r + 4000000

> homog_deq:=convert(LDO(Q)(t)=0, diff);

                  / 2      \
                  |d       |        /d      \
     homog_deq := |--- Q(t)| + 5000 |-- Q(t)| + 4000000 Q(t) = 0
                  |  2     |        \dt     /
                  \dt      /

> nonhomog_deq:=convert(LDO(Q)(t)=12, diff);

                   / 2      \
                   |d       |        /d      \
   nonhomog_deq := |--- Q(t)| + 5000 |-- Q(t)| + 4000000 Q(t) = 12
                   |  2     |        \dt     /
                   \dt      /

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

                    IC := {Q(0) = 0, D(Q)(0) = 0}


        / 2      \
        |d       |        /d      \
IVP := {|--- Q(t)| + 5000 |-- Q(t)| + 4000000 Q(t) = 12, Q(0) = 0,
        |  2     |        \dt     /
        \dt      /

    D(Q)(0) = 0}

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

                          rt := -4000, -1000

> Qc:=unapply(c1*exp(-1000*t) + c2*exp(-4000*t), t);

             Qc := t -> c1 exp(-1000 t) + c2 exp(-4000 t)

> Q:=Qc;

                               Q := Qc

> check:=simplify(homog_deq);

                            check := 0 = 0

# Ex 2-2: Ensure that (the terms of) g(t), 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);

                                  12

# Ex 2-3: Set up the subproblems, one for each term in g(t).
> unassign('Q');
> nh1:=convert(LDO(Q)(t)=12, diff);

               / 2      \
               |d       |        /d      \
        nh1 := |--- Q(t)| + 5000 |-- Q(t)| + 4000000 Q(t) = 12
               |  2     |        \dt     /
               \dt      /

# Ex 2-4: For each subproblem, assume the form of the particular
# solution to the corresponding nonhomogeneous subproblems.
> Qp1:=unapply(a, t);

                            Qp1 := t -> a

# Ex 2-5: Find the particular solution to each subproblem. Then sum them
# to form a particular solution to the full nonhomogeneous problem.
> simplify(subs(Q(t)=Qp1(t), nh1));

                            4000000 a = 12

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

                        sol := {a = 3/1000000}

> Qp1:=unapply(subs(sol, Qp1(t)), t);

                           Qp1 := 3/1000000

> Qp:=Qp1;

                           Qp := 3/1000000

# Ex 2-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.
> Q:=unapply(Qc(t)+Qp(t), t);

       Q := t -> c1 exp(-1000 t) + c2 exp(-4000 t) + 3/1000000

> check:=simplify(nonhomog_deq);

                           check := 12 = 12

# Ex 2-7: Resolve constants using ICs.
> IC; c_sols:=solve(", {c1, c2});

          {-1000 c1 - 4000 c2 = 0, c1 + c2 + 3/1000000 = 0}


                                                 -1
               c_sols := {c2 = 1/1000000, c1 = ------}
                                               250000

> final_solution; Q:=unapply(subs(c_sols, Q(t)), t);

                            final_solution


Q :=

    t -> - 1/250000 exp(-1000 t) + 1/1000000 exp(-4000 t) + 3/1000000

> final_check:=simplify(IVP);

                   final_check := {0 = 0, 12 = 12}

> #plot(u(t), t=0..0.6);
# Ex 2-8: Other requests...
> charge_after_one_thousandth_of_a_second:=Q(0.001);

                                                              -5
     charge_after_one_thousandth_of_a_second := .1546797874 10

> charge_after_one_hundredth_of_a_second:=Q(0.01);

                                                              -5
      charge_after_one_hundredth_of_a_second := .2999818400 10

> general_charge:=Q(t);

general_charge :=

    - 1/250000 exp(-1000 t) + 1/1000000 exp(-4000 t) + 3/1000000

> Limit(Q(t), t=infinity); limiting_charge:=value(");

     lim
t -> infinity

    - 1/250000 exp(-1000 t) + 1/1000000 exp(-4000 t) + 3/1000000


                     limiting_charge := 3/1000000

>