Direction Field Plots

Note also that intuition gained from a direction field plot needs to be nailed down by careful analysis. Just giving a picture will not garner full credit.

To get a direction field plot of the equation in #9, above, with the solution through the point [0,-1] shown, we would use

> with(DEtools):
> de:=diff(y(x),x)=x-y(x);
                             d
                       de := -- y(x) = x - y(x)
                             dx

> DEplot(de,y(x),x=0..3,y=-1..2,[[0,-1]],linecolor=black,scaling=constrained);

(Export to Mapletext strips out graphs, so you must do this in Maple to see the picture.)

In the direction field plot, you can trace out the solution passing through [0,1] by following the arrows with your finger. Observe that you do in fact see a relative minimum appearing, as indicated in c) and that the solution does not cross the asymptote y = x - 1. Again note that pictures do not constitute a full credit solution to this problem.

The Euler Approximation Method

The curve that you see in the direction field plot was derived by the forward Euler (Oi-ler, not You-ler) Approximation Method. You are expected to be able to use this method both in the tabular (four column table) method and by Maple.

x      y      m     Delta y
0      1     -1     -0.1
0.1   0.9   -0.8    -0.08
0.2   0.82  -0.62   -0.062

The Maple version of the Euler approximation is

> sol:=dsolve({de,y(0)=1},numeric,method=classical[foreuler],stepsize=0.1);
               sol := proc(x_classical)  ...  end proc

> sol(0.0);sol(0.1);sol(0.2);

                         [x = 0., y(x) = 1.]


                [x = 0.1, y(x) = 0.900000000000000024]


                [x = 0.2, y(x) = 0.820000000000000064]

To find the approximate value of y at x = 2, one would have to write 20 lines in the tabular method, since we start at x = 0 and go to x = 2 in steps of size 0.1. With the Maple command, we would simply say sol(2.0).

To plot the Euler approximation, we use the odeplot command in the plots package.

> with(plots):
> odeplot(sol,[x,y(x)],0..2);
Warning, the name changecoords has been redefined

You can find more about the classical numeric approximations by searching Help > Maple Help and entering the word "classical," from which you can get to
dsolve/numeric/classical - numerical solution of ordinary differential equations.

We will also use two other classical Maple methods in Chapter 3, in which [foreuler] is replaced by [heunform] or [rk4]. The alternative methods give considerably better approximations for a given stepsize, but are considerably more complicated and will not be done by hand.

As a side comment, when you use the DEplot command, one might ask what the stepsize is that is used in the Euler method. Maple's default is 1/20-th the x-range. That is frequently too large, and details are missed. Note that if your plots have sharp corners in them (something that cannot happen if the curve has a derivative) then you may want to set your own stepsize by including a ",stepsize = 0.01" right before the final parenthesis in the DEplot command.

Moreover, you should be aware of division by zero when you try to use the DEplot form. The best procedure is to plot just a very small range on either side of your initial value, checking that your syntax is correct and getting a graph; and then gradually increase the x-range of your plot. If something "breaks," and the plot fails to appear with the larger x-range, then you can go back and study the function to see why. Usually, it is because you tried to plot through a vertical asymptote of some part of the differential equation.

Day 4:

In Day 4 we quickly reviewed Maple worksheets showing how to find (usually implicit) solutions for variables separable (Section 2.2) and (always explicit) solutions for first order linear equations (Section 2.3).

Here is an example of a variables separable problem, p. 46 #15.

Begin by writing the equation with dy/dx. Then separate the variables. 
Note that there is no y(x) in the Maple integral - just y.
> restart:
> de:=y(x)^(-2)*diff(y(x),x)=-exp(cos(x))*sin(x);

                       d
                       -- y(x)
                       dx
                 de := ------- = -exp(cos(x)) sin(x)
                            2
                        y(x)

> eq:=int(1/y^2,y)=int(-exp(cos(x))*sin(x),x)+c;

                    eq := - 1/y = exp(cos(x)) + c

In general you will not be able to isolate the dependent variable in a 
solution to a variables separable problem, and you should check using the 
IMPLICIT solution method. In this particular problem, however, one can 
isolate y to get an explicit solution.
> sol:=y(x)=solve(%,y);

                                         1
                    sol := y(x) = ----------------
                                  -exp(cos(x)) - c

If you looked in the back of the book (which is NOT how you check your 
answers, right?!), you would notice that the book has a "+c", rather than 
a "-c". The negative of an unknown constant is still an unknown constant.

As an example of the technique for first order linear equations, we do p. 55, #16.

Always enter the original form of the differential equation. That is the 
form that we will check with, since doing so will catch manipulation 
errors in passing to the standard form equation. 

Recall that standard form for FOL means that the coefficient of the 
derivative is one and that the dependent variable is on the lhs.

> restart:
> original:=(x^2+1)*diff(y(x),x)=x^2+2*x-1-4*x*y(x);

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

> stform:=diff(y(x),x)+4*x/(x^2+1)*y(x)=(x^2+2*x-1)/(x^2+1);

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

The coefficient of the dependent variable always includes the sign. 
Otherwise, the integral will be off by a minus, and the value of mu will 
be the reciprocal of the correct one.

> p:=4*x/(x^2+1);
> int(%,x);
> mu:=exp(%);

                                   4 x
                             p := ------
                                   2
                                  x  + 1


                                   2
                             2 ln(x  + 1)


                                   2     2
                           mu := (x  + 1)

Applying the formula from p. 50, equation 8, we get

> sol:=y(x)=1/mu*(int(mu*rhs(stform),x)+c);

                                     5        4    2
                           -x + 1/5 x  + 1/2 x  + x  + c
             sol := y(x) = -----------------------------
                                       2     2
                                     (x  + 1)

Checking the explicit answer:

> subs(sol,original):
> simplify(%);

         2    6                   4
    -15 x  - x  + 20 x c + 5 - 5 x  - 10 x
  - -------------------------------------- =
                     2     2
                 5 (x  + 1)

               2    6                   4
          -15 x  - x  + 20 x c + 5 - 5 x  - 10 x
        - --------------------------------------
                           2     2
                       5 (x  + 1)

Rather than checking that all six terms in the numerator on the lhs are 
equal, respectively, to all six terms on the rhs, we just subtract and 
look for zero.

> lhs(%)-rhs(%);

                                  0

Day 5:

On Day 5, we worked out p. 46 #27 and p. 55 #25 in detail. When definite integrals are used, the lower limits are the initial values, and the upper limits are the variables in which the answer will be expressed. There is no constant of integration.

In connection with these examples, we reviewed the concept of a dummy variable, which was first introduced in your calculus courses. The variable used in the body of an integral can be replaced with any other symbol without changing the value of the integral.

Note that in the Lab Manual, on p. 26 #5, there is a formula similar to NSS p. 50 (8), which is written in terms of definite integrals and which gives the solution to a first order linear initial value problem in ONE step.

We also went through the Maple worksheet on Bernoulli equations. Here is an additional example to show that the patterns which we observed in class occur on all Bernoulli problems. The Bernoulli method is a 3 by 5 card.

Note that the Bernoulli method is the only technique that we will learn in Chapter 2, besides the variables separable and first order linear ones that you learned in Calculus 2.

 
Bernoulli Equation Example:

Note that the equation is nonlinear, does not separate, but has the 
Bernoulli form in which the lhs looks like a first order linear, while the 
rhs is a function of t times a power of y.

> de:=diff(y(t),t)+y(t)=t*y(t)^3;

                         /d      \                3
                   de := |-- y(t)| + y(t) = t y(t)
                         \dt     /

The Bernoulli substitution comes from the power of y.
> subst:=y(t)=v(t)^(1/(1-3));

                                          1
                       subst := y(t) = -------
                                           1/2
                                       v(t)

At first glance, the new equation in v(t) appears to be nonlinear. BUT - 
the power of v(t) is the same on both sides of the equation; and we can 
multiply to get rid of it.
> subs(subst,de):
> simplify(%);
> %*v(t)^(3/2);

                       /d      \
                      -|-- v(t)| + 2 v(t)
                       \dt     /               t
                  1/2 ------------------- = -------
                                3/2             3/2
                            v(t)            v(t)


                           /d      \
                      -1/2 |-- v(t)| + v(t) = t
                           \dt     /

Put the new first order linear equation in standard form.
> st:=%*(-2);

                         /d      \
                   st := |-- v(t)| - 2 v(t) = -2 t
                         \dt     /

Create the integrating factor mu:
> p:=-2;
> int(%,t);
> mu:=exp(%);

                               p := -2


                                 -2 t


                           mu := exp(-2 t)

Apply the formula on p. 50, (8). Then make it pretty.
> 1/mu*(int(mu*rhs(st),t)+c);
> expand(%);
> vsol:=v(t)=simplify(%);

                     1/2 (1 + 2 t) exp(-2 t) + c
                     ---------------------------
                              exp(-2 t)


                                         2
                         1/2 + t + exp(t)  c


                 vsol := v(t) = 1/2 + t + exp(2 t) c

Go back through the substitution to get the answer to the original 
problem.
> ysol:=subs(vsol,subst);

                                          1
               ysol := y(t) = -------------------------
                                                    1/2
                              (1/2 + t + exp(2 t) c)

Check that the answer is a valid explicit solution to the original 
problem.
> subs(ysol,de):
> simplify(%);

                          4 t
  ---------------------------------------------------- =
                                                   1/2
  (1 + 2 t + 2 exp(2 t) c) (2 + 4 t + 4 exp(2 t) c)

                                4 t
        ----------------------------------------------------
                                                         1/2
        (1 + 2 t + 2 exp(2 t) c) (2 + 4 t + 4 exp(2 t) c)

Day 6:

In Day 6, we began work in Chapter 3 of NSS. In particular, we solved examples p. 98 #4 and p. 99 #14.

With mixing problems, such as #4, the unknown, q(t), is always the mass, never the concentration. The differential equation is

dq/dt = conc-in * rate-in - conc-out * rate-out.

Three out of the four entries in the equation are read from the problem. The other, the concentration out, is computed as a ratio q(t)/volume, where the volume in the tank at time t is the volume initially in the tank, plus t times the difference between the in-flow and out-flow rates.

#14 is an example of a Malthusian population model. The differential equation is:

dp/dt = k p

Once Maple has solved for the solution, we are presented with a problem, since the solution has a variable k that is usually not given in the problem. One needs to create one equation using data from the problem, and then solve that one equation to find k. Thus, problems such as these typically have an initial population, which is the initial value, plus the population at one additional time, which allows us to take the solution to the differential equation and create an auxiliary equation with which to solve for k.

To make sure that we as a class always do the same thing on a problem, the statement of the problem will always indicate which of the two models, Malthusian or logistic, is to be applied.

In the logistic model, solutions contain not just one, but two variables, a and b; so in addition to the initial population, the population at two additional times is given, to allow us to create two equations to solve for the two variables. Maple syntax to solve two equations is:

absol:=solve({eq1,eq2},{a,b});

The Malthusian model also applies to two other situations, the growth of money (leading to the formula for continuously compounded interest) and radioactive decay. (Do NOT forget that we must have a differential equation. No credit for memorized formulas for compounded interest or radioactive decay!).

In both of these situations, the value for k is often determined by a time to double or a time to halve. If p(0) = P is the initial amount, then we get k from the time to double by 2P = P exp(k TD). Similarly, if p(0) = P is the initial amount, then we get k from the time to halve by 1/2 P = P exp(k TH).

Day 7:

Here is another example of a mixing problem.

p. 99 #7.

q1(t) is the number of kilograms in the first tank at time t.

> de1:=diff(q1(t),t)=0*3-3*q1(t)/60;
> inits1:=q1(0)=Q;
> sol1:=dsolve({de1,inits1},q1(t));
                         d            1       
                 de1 := --- q1(t) = - -- q1(t)
                         dt           20      
                      inits1 := q1(0) = Q
                                      /  1   \
                 sol1 := q1(t) = Q exp|- -- t|
                                      \  20  /

q2(t) is the number of kilograms in the tank at time t.

> de2:=diff(q2(t),t)=Q*exp(-t/20)/60*3-q2(t)/60*3;
> inits2:=q2(0)=0;
                 d          1       /  1   \   1       
         de2 := --- q2(t) = -- Q exp|- -- t| - -- q2(t)
                 dt         20      \  20  /   20      
                      inits2 := q2(0) = 0

Here is the solution to the amount of salt in tank 2 as a funtion of time.

> sol2:=dsolve({de2,inits2},q2(t));
                               1     /  1   \    
               sol2 := q2(t) = -- exp|- -- t| Q t
                               20    \  20  /    

The water in tank 2 is the saltiest when the derivative is 0.

> diff(rhs(%),t);
> eq:=simplify(%)=0;
               1     /  1   \       1       /  1   \
            - --- exp|- -- t| Q t + -- Q exp|- -- t|
              400    \  20  /       20      \  20  /
                      1       /  1   \             
             eq := - --- Q exp|- -- t| (t - 20) = 0
                     400      \  20  /             
> tsol:=solve(eq,{t});
                        tsol := {t = 20}

The maximum salt in tank 2 is attained at time t = 20 minutes, since the 
derivative changes from positive to negative. (First derivative test.)

The ratio of the maximum salt in tank2 to the original amount of salt in 
tank 1 is

> subs(tsol,rhs(sol2))/Q;
                            exp(-1)

In Day 7 we looked at examples of Newton's Law of Cooling and Newton's Law of Motion. Here are some examples showing the pattern.

p. 99 #15

Note that at first the solution is rather useless, since there are two 
variables that we don't know.
> restart:
> de:=diff(p(t),t)=a*p(t)-b*p(t)^2;
> inits:=p(0)=300;
> sol:=dsolve({de,inits},p(t));

                         d                        2
                   de := -- p(t) = a p(t) - b p(t)
                         dt


                         inits := p(0) = 300


                                       300 a
         sol := p(t) = -------------------------------------
                       300 b + exp(-a t) a - 300 exp(-a t) b

Create two equations with the two extra population levels to solve for a 
and b.
> eq1:=1200=subs(t=5,rhs(sol));
> eq2:=1500=subs(t=10,rhs(sol));
> absol:=solve({eq1,eq2},{a,b});

                                       300 a
         eq1 := 1200 = -------------------------------------
                       300 b + exp(-5 a) a - 300 exp(-5 a) b


                                       300 a
        eq2 := 1500 = ---------------------------------------
                      300 b + exp(-10 a) a - 300 exp(-10 a) b


                                            11
             absol := {a = 1/5 ln(15), b = ----- ln(15)}
                                           84000

With the complete solution, find the population in 2010.
> sol:=subs(absol,sol);
> subs(t=40,%);
> evalf(%);

                                     60 ln(15)
      sol := p(t) = -------------------------------------------
                    11
                    --- ln(15) + 9/56 exp(-1/5 ln(15) t) ln(15)
                    280


                                  60 ln(15)
           p(40) = ---------------------------------------
                   11
                   --- ln(15) + 9/56 exp(-8 ln(15)) ln(15)
                   280


                         p(40) = 1527.272725

p. 100 #25

Note that the solution is initially useless because it contains an unknown 
variable.
> de:=diff(p(t),t)=k*p(t);
> inits:=p(0)=P;
> sol:=dsolve({de,inits},p(t));

                              d
                        de := -- p(t) = k p(t)
                              dt


                          inits := p(0) = P


                       sol := p(t) = P exp(k t)

Create one equation to solve for the one unknown. Then make a useful 
solution.
> eq:=1/2*P=subs(t=5600,rhs(sol));
> ksol:=solve(%,{k});
> sol:=subs(ksol,sol);

                      eq := P/2 = P exp(5600 k)


                     ksol := {k = -1/5600 ln(2)}


                 sol := p(t) = P exp(-1/5600 ln(2) t)

Now solve for the age of the sample.
> eq1:=2/100*P=rhs(sol);
> tsol:=solve(%,{t});
> evalf(%);

                         P
                 eq1 := ---- = P exp(-1/5600 ln(2) t)
                         50


                                   5600 ln(50)
                      tsol := {t = -----------}
                                      ln(2)


                          {t = 31605.59466}

p. 107 #5

Note the initially useless solution.
> de:=diff(T(t),t)=k*(16-T(t));
> inits:=T(0)=34.5;
> sol:=dsolve({de,inits},T(t));

                          d
                    de := -- T(t) = k (16 - T(t))
                          dt


                         inits := T(0) = 34.5


                  sol := T(t) = 16 + 37/2 exp(-k t)

Create an equation to solve for k, then get a useful solution.
> eq:=33.7=subs(t=1,rhs(sol));
> ksol:=solve(%,{k});
> sol:=subs(ksol,sol);

                    eq := 33.7 = 16 + 37/2 exp(-k)


                     ksol := {k = 0.04420609250}


            sol := T(t) = 16 + 37/2 exp(-0.04420609250 t)

Solve for t. The time of death was 2.87 hours before noon.
> eq1:=37=rhs(sol);
> tsol:=solve(%,{t});

             eq1 := 37 = 16 + 37/2 exp(-0.04420609250 t)


                      tsol := {t = -2.867290422}

p. 108 #11.

Note that there are two Newton Cooling terms, since there are two heat 
flows in this problem. One heat flow is from the air outside the van to 
the air inside the van, and a second from the air inside the van to the 
A/C unit itself.

Note that we do not need to solve for K and K[u], since we have time 
constant data, rather than the temperatures at two additional times. The 
time constant is how long one has to wait for the temperature 
difference to decay to 1/e of its initial value, i.e., the
"e-th life."

With radioactive decay, we solve for  k  from the half life: -ln(2) = - 
k*t, so that k = ln(2)/t. Had we used the "e-th life" there, the result 
would have been  - ln(e) = - k*t , or  - 1 = - k*t, so that  k = 1/t, 
which is what we do with time constant problems. (By substituting 
p(t) = A - T(t), Newton's Law of Cooling  dT/dt = k (A - T) becomes dp/dt 
= - k p, which is radioactive decay!)

Note p. 106, K[1] = K + K[u], which allows us to solve for K[u]:  
3 = 1/2 + K[u], so K[u] = 5/2.

> de:=diff(T(t),t)=1/2*(35-T(t))+5/2*(16-T(t));
> inits:=T(0)=55;
> sol:=dsolve({de,inits},T(t));

                          d
                    de := -- T(t) = 115/2 - 3 T(t)
                          dt


                          inits := T(0) = 55


                sol := T(t) = 115/6 + 215/6 exp(-3 t)

Now solve for the time at which the van reaches 27 degrees. Since the time 
constant data was in hours, the 0.5 represents hours, too.
> eq:=27=rhs(sol);
> tsol:=solve(%,{t});
> evalf(%);

                  eq := 27 = 115/6 + 215/6 exp(-3 t)


                                           47
                      tsol := {t = -1/3 ln(---)}
                                           215


                          {t = 0.5068301419}

p. 115 #3. 

The acceleration is negative since we fall down, not up, in the normal 
coordinate system. The velocity is negative, so to make the air resistance 
positive, to work against gravity, we need a minus before 50 v(t).
> de:=500*diff(v(t),t)=500*(-9.81)-50*v(t);
> inits:=v(0)=0;
> sol:=dsolve({de,inits},v(t));

                         /d      \
               de := 500 |-- v(t)| = -4905.00 - 50 v(t)
                         \dt     /


                          inits := v(0) = 0


                                981   981        t
                sol := v(t) = - --- + --- exp(- ----)
                                10    10         10

The distance traveled is the integral of the velocity. Use the definite 
integral, and solve for the landing time T.
> eq:=-1000=int(rhs(sol),t=0..T);
> tsol:=solve(%,{T});
> evalf(%);

                            981 T              T
            eq := -1000 = - ----- - 981 exp(- ----) + 981
                             10                10


                                -1981     19810
  tsol := {T = 10 LambertW(-exp(-----)) + -----},
                                 981       981

                                  -1981     19810
        {T = 10 LambertW(-1, -exp(-----)) + -----}
                                   981       981


                {T = 18.64374931}, {T = -11.55158478}

Throw out the negative answer.

You should now be able to do all the problems in the HW list through those on p. 118.

You can ignore Section 3.5. We won't have time for it, since I want to get to Chapters 4 and 6, where we solve differential equations of orders higher than one.

In addition, you should be able to do the HW problems in Section 3.6 and 3.7, on p. 132 and p. 142. (Ignore all problem instructions about tolerances.) The modified Euler and Runge-Kutta approximations will ONLY be done in Maple in this class (no hand work!); and that makes the problems easy, since there is only one word that is different in the dsolve/numeric command and no difference at all in the odeplot command.

In Calculus 1, when you used Riemann sums to approximate the numeric value of an integral, you learned several approximation methods, like the left hand rule, the trapezoidal rule, and Simpson's rule. Simpson's rule, based on approximating sections of the curve by parabolas, was the most complicated to understand; but it gave the best approximation for a given stepsize partition.

A similar situation arises in differential equations, in that there are numerous approximation methods. The forward Euler method is easy to understand, but does not give a good approximation unless one uses a very small stepsize. The Runge-Kutta method is hard to understand, but gives a better approximation.

What you need to carry away is that the forward Euler method is order 1, meaning that if you divide the stepsize by a factor of 10 (which of course takes ten times as many steps and ten times as long to run), the error in the approximation will be divided by about 10.

In contrast, with the modified Euler, dividing the stepsize by a factor of 10 will divide the error of the approximation not by 10, but by 10^2. Hence, the modified Euler method is said to be "of order 2."

Similarly, the 4-th order Runge-Kutta method has the property that if we divide the stepsize by 10, the error in the approxmation is divided not by 10, but by 10^4.

In SAT Texian, rk4:differential equations::Simpson's rule:integrals.

Here are some examples:

p. 133 #9

The only difference in the Maple syntax between forward Euler, modified 
Euler, and Runge-Kutta is the one word.

A useful trick: Do the lookup yourself and notice that what I did was to 
do a global text seach on classical, to get 
"dsolve/numeric/classical - numerical solution of ordinary differential 
equations"
and then I copied one of their examples into my worksheet, replacing their 
data by mine. That helps me get the syntax right.

Maple is a programming language, like C or MatLab, and you can write for 
and do-while loops in it. If you want to use a bit of fancy Maple, we can 
make the array list with one of the loops. (Or, you can simply type in the 
values, if you prefer. Either way works.)

> 0.2*p$p=0..10;
> mylist:=[%];
> 

         0., 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0


   mylist := [0., 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0]

An array answer.
> de:=diff(y(x),x)=x+3*cos(x*y(x));
> inits:=y(0)=0;

                        d
                  de := -- y(x) = x + 3 cos(x y(x))
                        dx


                          inits := y(0) = 0

Here is the copied syntax from Maple Help, but with my data replacing 
theirs:
> ans := dsolve({de, inits}, numeric,
>                method=classical[heunform],
>                output=array(mylist),
>                stepsize=0.02);

                       [          [x, y(x)]          ]
                       [                             ]
                       [[0.              0.         ]]
                       [[                           ]]
                       [[0.2    0.619074340748844932]]
                       [[                           ]]
                       [[0.4    1.25021747100689407 ]]
                       [[                           ]]
                       [[0.6    1.77373367436592200 ]]
                       [[                           ]]
                       [[0.8    2.03954325981254403 ]]
                ans := [[                           ]]
                       [[1.0    2.05163395949820471 ]]
                       [[                           ]]
                       [[1.2    1.92339233977137014 ]]
                       [[                           ]]
                       [[1.4    1.74747202844878880 ]]
                       [[                           ]]
                       [[1.6    1.57233033146240109 ]]
                       [[                           ]]
                       [[1.8    1.41901669728271273 ]]
                       [[                           ]]
                       [[2.0    1.29418168290067848 ]]

Then we plot. (Of course, the plot will only show in the Maple worksheet, 
not in the text version in News of the Day.) Note that the
odeplot syntax for modified Euler and Runge-Kutta does not change at all 
from that which we learned with the forward Euler method.
> with(plots):
> odeplot(ans,[x,y(x)],0..2);
Warning, the name changecoords has been redefined

p. 142 #9

We're going to deviate a bit from the book, so that you can compare dsolve 
and dsolve/numeric, in three of the classsical dsolve/numeric forms. 
Again, note that the syntax for the three methods differs only in the 
single word.
> de:=diff(y(x),x)=x+1-y(x);
> inits:=y(0)=1;
> sol0:=dsolve({de,inits},y(x));
> 
sol1:=dsolve({de,inits},y(x),numeric,method=classical[foreuler],stepsize=0.25);
> 
sol2:=dsolve({de,inits},y(x),numeric,method=classical[heunform],stepsize=0.25);
> 
sol3:=dsolve({de,inits},y(x),numeric,method=classical[rk4],stepsize=0.25);

                           d
                     de := -- y(x) = x + 1 - y(x)
                           dx


                          inits := y(0) = 1


                      sol0 := y(x) = x + exp(-x)


               sol1 := proc(x_classical)  ...  end proc


               sol2 := proc(x_classical)  ...  end proc


               sol3 := proc(x_classical)  ...  end proc

A floating point decimal approximation of the formula solution:
> subs(x=1,rhs(sol0));
> evalf(%);

                             1 + exp(-1)


                             1.367879441

Now compute the forward euler, modified Euler, and finally the fourth 
order Runge-Kutta method, all evaluated at x = 1. Notice which one is 
closest when all use the same stepsize.
> sol1(1);
> sol2(1);
> sol3(1);

                 [x = 1., y(x) = 1.31640625000000000]


                 [x = 1., y(x) = 1.37252902984619140]


                 [x = 1., y(x) = 1.36789419940674861]

Day 8:

In Day 8 we began Chapters 4 (second order linear) and 6 (linear higher than second order ). Since differential equations in both chapters are handled exactly the same way, both will be studied at the same time and homework from both chapters can be done at the same time.

There are precisely two kinds of equations that we will learn to solve, constant coefficient (the coefficients are constants - have no x's) and Cauchy-Euler (the coefficients are a constant times a power of x, where the power exactly matches the order of the derivative of y for which it is a coefficient.

Here are examples of each of the two types, taken from p. 15 in the HW.

Both are homogeneous, meaning that the right hand side is zero. We will begin with homogeneous equations, since the solutions to linear equations that are not homogeneous are built from the homogeneous solutions. Nonhomogeneous solutions will be covered next week.

Type:  Linear homogeneous (rhs = 0) constant coefficient equation.
> restart:
> de:=diff(y(x),x$3)+3*diff(y(x),x$2)+2*diff(y(x),x)=0;

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

The coefficients are 1, 3, and 2, all of which are constants. For constant 
coefficient equations, the form of an explicit solution is always y(x) = 
exp(m*x).
> sol:=y(x)=exp(m*x);
> subs(sol,de);
> simplify(%);

                        sol := y(x) = exp(m x)


       / 3          \     / 2          \
       |d           |     |d           |     /d          \
       |--- exp(m x)| + 3 |--- exp(m x)| + 2 |-- exp(m x)| = 0
       |  3         |     |  2         |     \dx         /
       \dx          /     \dx          /


                                 2
                    m exp(m x) (m  + 3 m + 2) = 0

The polynomial has the same degree as the order of the differential 
equation. The polynomial equation is called the characteristic equation. 
We let Maple solve the cubic equation.
> char_eq:=%/exp(m*x);
> solve(%,{m});

                                  2
                   char_eq := m (m  + 3 m + 2) = 0


                     {m = 0}, {m = -1}, {m = -2}

The general solution (which includes all three constants of integration) 
is
> sol:=y(x)=c1*exp(0*x)+c2*exp(-1*x)+c3*exp(-2*x);

             sol := y(x) = c1 + c2 exp(-x) + c3 exp(-2 x)

Type: Linear homogeneous(rhs = 0) Cauchy-Euler equation.
> de:=3*x^2*diff(y(x),x$2)+11*x*diff(y(x),x)-3*y(x)=0;

                    / 2      \
                  2 |d       |        /d      \
         de := 3 x  |--- y(x)| + 11 x |-- y(x)| - 3 y(x) = 0
                    |  2     |        \dx     /
                    \dx      /

The coefficient of the second derivative is 3x^2, which has a second power 
of x; the coefficient of the coefficient of the first derivative
is  11x, which has a first power of x; and the coefficient of the "zeroth 
derivative," y(x) has no powers of x. For Cauchy-Euler equations, the form 
of an explicit solution is always y(x) = x^m.
> sol:=y(x)=x^m;
> subs(sol,de);
> simplify(%);

                                          m
                           sol := y(x) = x


                    / 2    \
                  2 |d    m|        /d   m\      m
               3 x  |--- x | + 11 x |-- x | - 3 x  = 0
                    |  2   |        \dx   /
                    \dx    /


                        m  2      m        m
                     3 x  m  + 8 x  m - 3 x  = 0

The polynomial has the same degree as the order of the differential 
equation. The polynomial equation is called the indicial equation. We let 
Maple solve the quadratic equation.
> %/x^m;
> ind_eq:=expand(%);
> solve(%,{m});

                        m  2      m        m
                     3 x  m  + 8 x  m - 3 x
                     ----------------------- = 0
                                m
                               x


                                  2
                     ind_eq := 3 m  + 8 m - 3 = 0


                         {m = 1/3}, {m = -3}

The general solution (which includes both constants of integration) is
> sol:=y(x)=c1*x^(1/3)+c2*x^(-3);

                                      (1/3)    c2
                    sol := y(x) = c1 x      + ----
                                                3
                                               x

If there are initial values, so that we need to solve for the constants of the equation to get the unique solution, then we do that as in p. 15 #22.

p. 15 #22. IVP: Second order linear homogeneous constant coefficient 
differential equation, with initial values y(0) = 2,  y'(0) = 1.
> de:=diff(y(x),x$2)+diff(y(x),x)-2*y(x)=0;

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

> sol:=y(x)=exp(m*x);
> subs(sol,de);
> simplify(%);
> 

                        sol := y(x) = exp(m x)


           / 2          \
           |d           |   /d          \
           |--- exp(m x)| + |-- exp(m x)| - 2 exp(m x) = 0
           |  2         |   \dx         /
           \dx          /


                                 2
                      exp(m x) (m  + m - 2) = 0

> char_eq:=%/exp(m*x);
> solve(%,{m});

                                  2
                      char_eq := m  + m - 2 = 0


                          {m = 1}, {m = -2}

The general solution is
> sol:=y(x)=c1*exp(1*x)+c2*exp(-2*x);

                sol := y(x) = c1 exp(x) + c2 exp(-2 x)

There are two constants of integration since we had to integrate twice to 
go from y'' to y. Since there are two(2) constants of integration to solve 
for, we will need two(2) initial conditions. Set up the equations and 
solve simultaneously.
> 2=subs(x=0,rhs(sol));
> eq1:=simplify(%);
> diff(sol,x);
> 1=subs(x=0,rhs(%));
> eq2:=simplify(%);

                      2 = c1 exp(0) + c2 exp(0)


                          eq1 := 2 = c1 + c2


                 d
                 -- y(x) = c1 exp(x) - 2 c2 exp(-2 x)
                 dx


                     1 = c1 exp(0) - 2 c2 exp(0)


                         eq2 := 1 = c1 - 2 c2

Solve simultaneously.
> csol:=solve({eq1,eq2},{c1,c2});

                     csol := {c2 = 1/3, c1 = 5/3}

Find the unique solution to the IVP.
> subs(csol,sol);

                  y(x) = 5/3 exp(x) + 1/3 exp(-2 x)

The reason that the solutions for these linear homogeneous differential equations are so easy is that the roots of the characteristic (respectively, indicial) equations are all different (no repeated roots) and all real (not complex numbers). In the next couple of weeks we will show how the technique indicated here extends to the case of repeated eigenvalues (roots of the polynomial equations) and complex eigenvalues.

The problem is, that in order to do the extension we need to know a lot more of how we came by these solutions. And the basis of that knowledge is theory.

So, on Day 8, we begin by asking some easy questions.

Question #1

When is an equation "homogeneous"?

Answer: Decide what the dependent variable is. Move all terms that involve the dependent variable to the lhs, move all terms that involve *only* the independent variable to the rhs.

x^2 y'' + sin(x) = cos(x) y - 3

is nonhomogeneous since after the rearrangement we get

x^2 y'' - cos(x) y = - sin(x) - 3,

and -sin(x) - 3 is not zero.

Question #2

When is an equation "linear"?

In the example that we worked for Question #1, the first step was to separate the equation into terms having y and terms not having y. The rhs is the terms that do not have y. They determine if the equation is homogeneous or not - but they have NO effect on whether or not the equation is linear. ONLY the lhs determines if the equation is linear.

In the example above, L[y] = x^2 y'' - cos(x) y. (Note that the rhs is *not* included in the operator L.) L[y] is a differential operator, since it describes what derivative operations we do on the dependent variable, y.

The operation L[y] will be a linear operation if two things are 
true:

1) L[y1 + y2] = L[y1] + L[y2]
2) L[c*y] = c*L[y].

If either one of those statements is false, then the operator is not linear, and therefore the differential equation is not linear.

You should be doing these by hand, not by Maple.

Example 1: yy' = x is not homogeneous and not linear.

The first step is to separate the have y's from the have nots. (That is already done here.)

Since the have-no-y side is different than 0, the equation is not homogeneous.

The operator to check is L[y] = y*y'. Note that nothing that does not have a y is included in the L.

To be linear, we would need L[y1 + y2] = L[y1] + L[y2]. But
L[y1 + y2] =
(y1 + y2)*(y1 + y2)'=
y1y1' + y1y2' + y2y1' + y2y2',
whereas L[y1] = y1y1' and L[y2]=y2y2'.

Since the operator fails the first test, we do not need to check the second case; and we can say that L is not linear, hence the differential equation is not linear.

Example 2: y' = xy' + x is not homogeneous, but it is linear.

The first step is to separate the have y's from the have nots.

y' - xy = x.

The equation is not homogeneous, since the have-not side is not zero.

The equation is linear. The operator is L[y] = y' - xy. The x is not included in the L.

Note that

L[cy] =
(cy)' - x(cy) =
cy' - cxy =
c(y' - xy) =
c*L[y]

To understand why is is possible to add answers and get new answers in the case of linear, homogeneous differential equations (i.e., to go from two answers to infinitely many by taking linear combinations with the parameters c), we need to understand exactly what linearity does for us. The old definition of "you can do this, but you cannot do that" is not enough. We need the two properties of a linear operator to understand. And in that vein, you should be prepared to decide if a differential equation is linear and/or homogeneous, based on the above criteria.

Question #3

When is a linear IVP guaranteed a unique solution?

There are three Existence Uniqueness Theorems in the book, so far:

The EU Theorem on p. 12 applies to both nonlinear and linear first order IVP's. The EU Theorem on p. 54 applies only to linear first order IVP's. And the EU Theorems on p. 165 and p. 318 apply to linear first and higher order IVP's, so they extend the theorem on p. 54.

In p. 15 #23, we had a formula solution, so we could actually solve for the values of the c's, and show that there was a unique solution. What the theorems do for us is to guarantee that there is only one solution to an IVP, even when we do not have a formula solution, meaning that we cannot actually solve for values of c.

It is easiest to study the three theorems as a triad, since we can then see what is similar and what is different. Notice the patterns!

p. 12 EU Theorem:

Standard form, y' = f(x,y).

Functions that must be continuous at the start: f and diff(f,y). (You need to show me the f and diff(f,y)!)

If we start at a good point, there is one and only one solution -- for a ways.

Remark: Sometimes the HW shows examples to help explain the EU theorems.

p. 15 #18 shows that the solution can blow up, but for NONLINEAR 
equations, exactly where depends on the initial value, not on the de.
(Observe that for linear equations, blowup ARE visible from the de!) 
p. 15 #29a shows that if we do not have both f and diff(f,y) continuous at 
the start, there can be two solutions to an IVP.
p. 15 #31b shows that if we do not have both f and diff(f,y) continuous at
the start, there can be no solutions to an IVP.

p. 54 EU Theorem:

Standard form, y' + p(x)y = g(x). Lead coefficient (on y') is 1, and all y's are on lhs.)

Functions that must be continuous at the start: p(x) and g(x).

If we start at a good point, there is one and only one solution that matches the initial data.

The solution is guaranteed not to blow up until we get to the first discontinuity of p(x) or g(x), (starting from the initial value). After it blows up, it does not start up again on the other side of the discontinuity.

p. 165 and p 318 EU Theorem:

Standard form: lead coefficient is 1 on the n-th derivative, all y's and y's, y'', etc., on the lhs.

Functions that must be continuous are all the coefficient functions and the function on the rhs.

Initial data means that the values of the function and the derivative are all at the SAME x value.

The theorem does NOT apply if you have two values of y at different x's, rather than a y and a y' at the same x. (See p. 167 #30 HW which helps you to understand the difference between initial value data and boundary value data. Boundary value data *can* allow multiple solutions; initial value data cannot.)

There is one and only one solution that matches the initial data, and it will not blow up until we get to the first discontinuity of one of the coefficient functions or a discontinuity of the rhs.

Example:
p. 324 #3. y''' - y'' + sqrt(x-1) y = tan(x).

The equation is in standard form, since the lead coefficient on y''' is 1 and y and all derivatives are on the lhs.

The starting point (5,1) is good since the coefficient functions 1, -1, sqrt(x-1), and the rhs function tan(x) are all continuous at x = 5.

Hence, we are guaranteed that there is one and only one solution that matches the initial data. (We do not have a formula solution, since the equation is neither constant coefficient nor Cauchy-Euler, so we cannot solve for values of c).

The starting point is x = 5. The solution is valid from 3 Pi/2 = 4.71 to 5 Pi/2 = 7.85. (It cannot just jump over the discontinuities!). We go up and down from the starting point to the first discontinuity of any of the functions.

Had we started at x = 1.5, which is also a point of continuity for each of the functions, then the solution would only have been valid on the interval (1,Pi/2) = (1,1.57), since the square root would have had problems on the left of 1.5, before the tangent function had problems at - Pi/2. We go up and down from the starting point to the first discontinuity of any of the functions.

What HW problems should I attempt in Chapters 4 and 6?

p. 167, all in HW list; p. 189, #19, 21, 26, 33, 35, 39.

In general, you should be able to give the general solution to any linear, homogeneous, constant coefficient or Cauchy-Euler as long as the characteristic/indicial equation has distinct real roots. If the roots repeat, you do not need to do the problem since we have not yet covered that material. Similarly, you do not yet know what to do with complex roots.)

Note that some of the problems deal with the EU Theorems or linearity via linear differential operators. You *should* be able to do problems that ask you about those two topics.

As an example to help you, here is p. 169 #49a.

de:=3*diff(y(x),x$3)+18*diff(y(x),x$2)+13*diff(y(x),x)-19*y(x)=0;
> sol:=y(x)=exp(m*x);
> subs(sol,de);
> simplify(%);
           / 3      \      / 2      \
           |d       |      |d       |      /d      \
   de := 3 |--- y(x)| + 18 |--- y(x)| + 13 |-- y(x)| - 19 y(x) = 0
           |  3     |      |  2     |      \dx     /
           \dx      /      \dx      /


                        sol := y(x) = exp(m x)


    / 3          \      / 2          \
    |d           |      |d           |      /d          \
  3 |--- exp(m x)| + 18 |--- exp(m x)| + 13 |-- exp(m x)|
    |  3         |      |  2         |      \dx         /
    \dx          /      \dx          /

         - 19 exp(m x) = 0


                            3       2
               exp(m x) (3 m  + 18 m  + 13 m - 19) = 0

> char_eq:=%/exp(m*x);
> fsolve(%,{m});

                             3       2
               char_eq := 3 m  + 18 m  + 13 m - 19 = 0


      {m = -4.831922263}, {m = -1.869273823}, {m = 0.7011960861}

> 
sol:=y(x)=c1*exp(-4.831922263*x)+c2*exp(-1.869273823*x)+c3*exp(.7011960861*x);

  sol := y(x) = c1 exp(-4.831922263 x) + c2 exp(-1.869273823 x)

         + c3 exp(0.7011960861 x)

Here's another.

Cauchy-Euler equation, with initial values:

> de:=x^3*diff(y(x),x$3)-3*x^2*diff(y(x),x$2)+6*x*diff(y(x),x)-6*y(x)=0;
> inits:=y(1)=4,D(y)(1)=9,(D@@2)(y)(1)=14;
         /  3      \        /  2      \                              
       3 | d       |      2 | d       |       / d      \             
de := x  |---- y(x)| - 3 x  |---- y(x)| + 6 x |--- y(x)| - 6 y(x) = 0
         |   3     |        |   2     |       \ dx     /             
         \ dx      /        \ dx      /                              
      inits := y(1) = 4, D(y)(1) = 9, @@(D, 2)(y)(1) = 14

Solution form for Cauchy-Euler equations:

> sol:=y(x)=x^r;
                                       r
                        sol := y(x) = x 

Compute indicial equation, solve for values of r:

> subs(sol,de):
  simplify(%);
  %/x^r:
  expand(%);
               r  3      r  2       r        r    
              x  r  - 6 x  r  + 11 x  r - 6 x  = 0
                     3      2               
                    r  - 6 r  + 11 r - 6 = 0

Eigenvalues:

> solve(%,r);
                            1, 2, 3

General solution to differential equation

> sol1:=y(x)=c1*x+c2*x^2+c3*x^3;
                                        2       3
              sol1 := y(x) = c1 x + c2 x  + c3 x 

Initial value problem solution:

> eq1:=4=subs(x=1,rhs(sol1));
  diff(sol1,x);
  eq2:=9=subs(x=1,rhs(%));
  diff(sol1,x$2);
  eq3:=14=subs(x=1,rhs(%));
                    eq1 := 4 = c1 + c2 + c3
                 d                             2
                --- y(x) = c1 + 2 c2 x + 3 c3 x 
                 dx                             
                  eq2 := 9 = c1 + 2 c2 + 3 c3
                     2                      
                    d                       
                   ---- y(x) = 2 c2 + 6 c3 x
                      2                     
                    dx                      
                    eq3 := 14 = 2 c2 + 6 c3

Initial value solution:

> csol:=solve({eq1,eq2,eq3},{c1,c2,c3});
  subs(csol,sol1);
                csol := {c3 = 2, c1 = 1, c2 = 1}
                                  2      3
                      y(x) = x + x  + 2 x 

Forbidden fruit:

> dsolve({de,inits},y(x));
                                  2      3
                      y(x) = x + x  + 2 x

Day 9:

Exam 1 was celebrated.

Day 10:

Solutions to linear homogeneous differential equations as vectors

We also began to discuss the concept of linear independance of solutions to linear, homogeneous differential equations.

In Day 8 we showed how to use a guess for the form of solutions to constant coefficient and Cauchy-Euler equations to get a finite number of solutions. However, since a differential equation has infinitely many solutions, having two or three solutions is not nearly the number of solutions that we need.

To get infinitely many solutions for linear, homogeneous equations, we can take linear combinations of the finitely many solutions. For example, with y" - 3 y' + 2y = 0, we get two solutions exp(x) and exp(2x) by substituting the explicit solution y(x) = exp(mx) into the differential equation, solving the resulting characteristic equation m^2 - 3m + 3 = 0, to get m = 1 and m = 2. By taking a linear combination, y(x) = c1 exp(x) + c2 exp(2x), we then get infinitely many solutions, since any choice of the parameters c1 and c2 yields a valid solution.

To understand why we are allowed to add solutions and get new solutions - as well as to understand why something different needs to be done when we have repeated eigenvalues in the case y" - 2y' + y = 0, where we get m = 1 and m = 1 as repeated eigenvalues, we need to observe that solutions to linear, homogeneous differential equations behave like (are) vectors.

In the plane, we can start with the two vectors i and j, which point in different directions and are independent of each other. Each linear combination v = c1 i + c2 j is again a vector; and since i and j are pointing in different directions (one is not a scalar multiple of the other), we get all possible vectors in the plane by taking such linear combinations.

In the case of a linear, homogeneous differential equation, such as y" - 3 y' + 2y = 0, we can start with two solutions exp(x) and exp(2x), which are independent of each other in the sense that exp(2x) is not a scalar multiple of exp(x). Each linear combination y(x) = c1 exp(x) + c2 exp(2x) is again a solution; and since exp(x) and exp(2x) are "pointing in different directions," (one is not a scalar multiple of the other), we get all possible solutions by taking all such linear combinations.

To see that solutions of all linear homogeneous differential equations (NOT just constant coefficient and Cauchy-Euler) can be added to again get solutions or multiplied by a number to again get a solution, we use the idea of a linear operator L[y] that we learned in Day 8.

To say that an operator L is linear means that L(c1 y1 + c2 y2) = c1 L[y1] + c2 L[y2].

If we have a homogeneous linear equation, we can write it as L[y] = 0, for a linear differential operator L.

If y1 and y2 are solutions to L[y] = 0, then L[y1] = 0 and L[y2] = 0, since when we subs and simplify, the result has to match the rhs.

However, it then follows that L[c1 y1 + c2 y2] = c1 L[y1] + c2 L[y2] = c1 0 + c2 0 = 0, so the linear combination c1 y1 + c2 y2 is again a solution to L[y] = 0 for all choices of the parameters c1 and c2.

(In a nutshell, this explains why we can solve linear equations but seldom solve nonlinear equations. Even if we were lucky enough to get two solutions to a nonlinear equation, we cannot take linear combinations of the two solutions that we have in order to get even more solutions.)

In other words, when we have a linear homogeneous differential equation,we can add solutions to get a new solution or multiply a solution by a constant to get a new solution.

But that behavior is just like what we observe with vectors in the plane!

Linear Independence of Vectors and Solutions

If we let our imagination run, we note further analogies.

In the plane, if we start with i and j, we can get all vectors by taking linear combinations c1 i + c2 j. However, i and j are not the only pair of vectors that we can use to get all the vectors in the plane. We can use v1 = 1/2 i + 1/2 j and v2 = 1/2 i - 1/2 j. All vectors in the plane can be represented as c1 v1 + c2 v2 for the appropriate choice of c1 and c2.

The same thing happens with homogeneous linear differential equations. All solutions to the equation y" - y = 0 can be written as c1 exp(x) + c2 exp(-x). But all solutions can also be written as a linear combination of the hyperbolic trig functions: c1 cosh(x) + c2 sinh(x).

There are infinitely many pairs of solutions that form a fundamental solution set, meaning that ALL solutions can be gotten by taking linear combinations of those two solutions.

But not every pair of vectors in the plane forms a fundamental set. Nor is every pair of solutions a fundamental solution set.

The pair formed by i and 2i do not form a fundamental set, since we will never get the vector j by taking a linear combination of i and 2i.

Similarly, for the homogeneous linear equation y" - y = 0, exp(x) and 2 exp(x) are both solutions, but one cannot write the solution exp(-x) as a linear combination of exp(x) and 2 exp(x).

Now if we are dealing with the plane, which is two dimensional, we can give a complete description of whether or not a pair of vectors forms a fundamental set in terms of scalar multiples. A pair {v1, v2} forms a fundamental set, meaning that all vectors can be written as a linear combination of v1 and v2 if and only if v2 is not a multiple of v1.

However, if we are in three dimensions, we need a slightly more general idea than just scalar multiple. Three vectors {v1, v2, v3} form a fundamental set if and only if there are no scalars c1, c2, and c3 for which the linear combination c1 v1 + c2 v2 + c3 v3 is the zero vector - except for the obvious scalars c1 = 0, c2 = 0, and c3 = 0.

Three vectors in 3 dimensions that have this property that the only linear combination which is equal to the zero vector is the obvious linear combination - are said to be linearly independent.

For example, the triple of vectors {i, j, k} is linearly independent since the only linear combination c1 i + c2 j + c3 k = 0 is the obvious linear combination with coefficients c1 = 0, c2 = 0, and c3 = 0.

When we have a third or higher order linear differential equation, our analogy should also change. The equation y''' - 6 y" + 11 y' - 6y = 0 has characteristic equation m^3 - 6m^2 + 11m - 6 = 0, and three eigenvalues m = 1, m = 2, and m = 3. Hence we have three solutions, {exp(x), exp(2x), and exp(3x)}. To see if all solutions to the equation can be made from linear combinations of these three solutions, we need to know if the set of solutions is a fundamental solution set. So the question is, is there any set of coefficients {c1, c2, and c3} except the obvious set c1 = 0, c2 = 0, and c3 =0, for which c1 exp(x) + c3 exp(2x) + c3 exp(3x) is identically zero?

(Being identically zero means that the expression is zero not just for one x or two x's, but for EVERY single x. The graph does not just cross the x-axis at a couple of points; the graph IS the x-axis.)

The Wronskian Determinant

There is a nice criterion to determine if three solutions are LINEARLY INDEPENDENT - meaning that the ONLY linear combination coefficients that will produce a graph that lies wholly on the x-axis is the obvious set where all coefficients are zero - and that criterion is called the Wronskian determinant. (See Wronski, Jozef Maria.)

To form the Wronskian matrix, we place the solutions on the top row, their derivatives on the second row, their third derivatives on the third row, etc., until we have a square matrix. We then compute the determinant of that square matrix, simplify - and if the Wronskian determinant is different than zero, we will know that the solutions are linearly independant, and that the solutions form a fundamental set, meaning that ALL solutions can be made from linear combinations of those solutions.

Remark: if you reverse the order of two solutions in the Wronskian matrix, the sign of the determinant will change. However, changing the sign does not affect whether or not the determinant is zero.

The Maple code for how to compute the Wronskian determinant is shown in the solutions to the third order linear homogeneous equation on Exam 1.

For example, if we wished to see if the solutions {exp(x), exp(2x), exp(3x)} form a fundamental solution set - meaning that we get all solutions by taking linear combinations of those three solutions - then we need to compute the Wronskian determinant and show that it is not zero.

> with(VectorCalculus):
> with(LinearAlgebra):Warning, the assigned names <,> and <|> now have 
a global binding

Warning, these protected names have been redefined and unprotected: *, 
+, ., D, Vector, diff, int, limit, series

Warning, the names &x, CrossProduct and DotProduct have been rebound

> FS:=[exp(x),exp(2*x),exp(3*x)];

                  FS := [exp(x), exp(2 x), exp(3 x)]

> W:=Wronskian(FS,x);

                   [exp(x)     exp(2 x)      exp(3 x) ]
                   [                                  ]
              W := [exp(x)    2 exp(2 x)    3 exp(3 x)]
                   [                                  ]
                   [exp(x)    4 exp(2 x)    9 exp(3 x)]

The determinant is not zero.
> Determinant(%);
> simplify(%);

                      2 exp(x) exp(2 x) exp(3 x)


                              2 exp(6 x)

Since the value 2exp(6 x) is not zero, the three solutions are linearly independant and form a fundamental solution set, so that ALL solutions to y''' - 6 y" + 11 y' - 6y = 0 can be written as y(x) = c1 exp(x) + c2 exp(2x) + c3 exp(3x).

On all future quizzes and exams, when you give the general solution to a constant coefficient or Cauchy-Euler equation, you will need to check that the Wronskian determinant of the solutions is not zero.

To see that the Wronskian determinant can be zero, meaning that we do not get all solutions, consider the constant coefficient equation y" - 2y' + y = 0. The characteristic equation is m^2 - 2m + 1 = (m-1)^2 = 0, so the eigenvalues are repeated: m = 1 and m = 1.

The resulting solutions {exp(x), exp(x)} are LINEARLY DEPENDENT, since the second solution exp(x) is a scalar multiple (by the scalar 1) of the first solution. We can see the same thing by noting that there is a linear combination of exp(x) and exp(x) which does not use only 0 coefficents, but which is identically zero: (1) exp(x) + (-1) exp(x). Finally, we see that the Wronskian determinant is zero, since

 
> with(VectorCalculus):
> with(LinearAlgebra):Warning, the names &x, CrossProduct and DotProduct 
have been rebound

Warning, the names &x, CrossProduct and DotProduct have been rebound

> FS:=[exp(x),exp(x)];

                        FS := [exp(x), exp(x)]

> W:=Wronskian(FS,x);

                            [exp(x)    exp(x)]
                       W := [                ]
                            [exp(x)    exp(x)]

The determinant is zero.
> Determinant(%);

                                  0

Since the solutions are linearly dependent, they do not form a fundamental solution set, and we cannot get all solutions by taking linear combinations c1 exp(x) + c2 exp(x).

Day 11:

On Day 11 we answered the following four questions:

Cramer's Rule

Question #1: What is Cramer's rule for solving systems of simultaneous 8th grade linear equations?

Cramer, Gabriel (1704-1752)

Cramer's rule allows one to solve systems of simultaneous linear algebraic equations by mean of determinants, provided the determinant of the coefficient matrix is different than zero. For example,

has solutions x = 2 and y = 3. We can find these values as ratios of determinants:

| 23  5 |
| 15  1 |
_________
|  4  5 |
|  6  1 |

and

| 4  23 |
| 6  15 |
_________
| 4   5 |
| 6   1 |

We will again have reason to recall Cramer's rule when we solve nonhomogeneous linear differential equations using variation of parameters.

We can see that the functions x and x^2 are linearly independent since the Wronskian determinant

| x  x^2 |
| 1  2 x |

is different than zero, so we can apply Cramer's rule. For if c1 x + c2 x^2 is identically zero, then its derivative c1 1 + c2 2x is also identically zero. Hence, the value of c1 is

| 0  x^2 |
| 0  2 x |
__________
| x  x^2 |
| 1  2 x |

which is 0, while the value for c2 is

| x    0 |
| 1    0 |
__________
| x  x^2 |
| 1  2 x |

which is also 0. Hence, the only coefficients c1 and c2 that will make the linear combination c1 x + c2 x^2 be identically zero are the obvious coefficients c1 = 0 and c2 = 0.

Question #3:. How does the EU Theorem guarantee that if the Wronskian determinant is nonzero, then we have all possible solutions?

The EU Theorem says that at a good x0, there will be exactly one solution through the point (x0,y0) with slope y'(x0). If someone gives us any initial data numbers, say y = 24, y' = 47, then there will be exactly one solution that matches those numbers at the given x0.

But *every* solution has to have *some* y and y' at x0. And by Cramer's Rule if the Wronskian determinant is different than zero at x0, then we can find a linear combination that matches those two numbers, then by uniqueness, every solution must be a linear combination (the coefficients that give that solution are the values of c1 and c2 that Cramer's Rule gives us.)

What does one do if the polynomial gives repeated roots? For the constant coefficient case, if one gets eigenvalues m = 1, m = 1, m = 1, m = 2, m = 2, m = 2, m = 2, say for a seventh order equation, then a Fundamental Solution set is

whereas, if we got the same values of m for a Cauchy-Euler equation, we would get as a Fundamental Solution set

FS:=[x,x*ln(x),x*ln(x)^2,x^2,x^2*ln(x),x^2*ln(x)^2,x^2*ln(x)^3];

A formal proof that the procedure indicated works is given in NSS, p. 329 and p. 330. You can recall what to multiply by for repeated eigenvalues by computing the partial with respect to m of exp(mx) for constant coefficent equations and the partial with respect to m of x^m for Cauchy-Euler equations.

Day 12:

On Day 12 we derived the formula for solutions to constant coefficient and Cauchy-Euler equations with complex eigenvalues.

Complex Conjugate Pairs

Note that if the coefficents of the original differential equation are real, which for us they always will be, then the coefficients of the characteristic/indicial equation will also have real coefficients.

One can easily check, by computing, that the conjugate of a product of complex numbers is the product of their conjugates. Moreover, one can also confirm that the conjugate of a sum (difference) of complex numbers is the sum (difference) of their conjugates.

It then follows that the conjugate of the characteristic (indicial) equation is the characteristic (indicial) equation with m replaced by the conjugate of m.

However, since the equation has real coefficients, taking the conjugate of the characteristic equation doesn't change anything. Therefore, the conjugates of the eigenvalues must already have been eigenvalues. Therefore, eigenvalues always come as complex conjugate pairs.

We use the fact that eigenvalues always come in complex conjugate pairs to see that solutions to the differential equation always come in complex conjugate pairs.

The key is

exp(i theta) = cos(theta) + i sin(theta),

which comes from the Maclaurin series for exp(z), with z = i theta - the even terms going into the even function cosine, and the odd terms going into the odd function sine.

The fact that the solutions for m come in complex conjugate pairs, with equal plus and minus values of i, explains why we also get complex conjugate pairs of solutions to the differential equation:

exp(m x) = exp( (a +/- bi)x ) = exp(a x) exp(i (+/-bx ))

where theta = +/- bx, so that by Euler's formula we get complex solutions y = exp(a x) (cos( bx ) +/- i sin( bx )).

Real and Imaginary Parts

For any complex conjugate pair z = f + gi and it's conjugate, z = f - gi, we can get the real part f (with no i!) by taking a linear combination 1/2 z + 1/2 z conjugate.

We can also get the imaginary part g (again, with no i!) by the linear combination 1/(2i) z - 1/(2i) z conjugate.

Since these are linear combinations of answers to linear homogeneous differential equations, the real part y = exp(a x) cos( bx ) and the imaginary part y = exp(a x) sin( bx ) are again answers.

Of course, the same thing happens with Cauchy-Euler equations, only now wehave x^( (a +/- bi) ) = exp ( (a +/- bi) ln(x) ) = x^a (cos(b ln(x)) +/- sin(b ln(x)), and again, the real and imaginary parts are linearly independent and form a fundamental solution set.

Caveat: it is not true that the cos comes from one of the m's and the sin from its conjugate m in the pair. You did not come just from your dad if you are male, nor just from your mom if you are female. Both solutions come from the pair.

It does not matter which sign you use on the sin term, since the sign is absorbed in the linear combination when you write the final answer.

Abel's Formula

For some of the problems, you need to know Abel's formula. Abel observed that if we have a second order linear homogeneous differential equation, and two solutions y1(x) and y2(x), then the Wronskian determinant is (obviously) also a function of W(x). He wondered if W(x) was, itself, the solution to a differential equation.

He was able to show (and really good 308 students can, too!) that if the differential equation is y" + p(x)y' + q(x)y = 0, then W(x) is the solution to W' + p(x)W = 0, where the p(x) is the same function!

Since we know how to solve the first order linear equation W' + p(x)W = 0 to get W(x) = C exp( - int(p(t),t=x0..x) ), we observe a striking result. If the constant C is zero, then the Wronskian determinant is zero for every x. And if the constant C is not zero, then since an exponential can never be zero, it follow that W(x) cannot be zero for any x: it is either zero for no x's or all x's!

You may find it useful to look at the first of the Exam 2 review sheets which are available in the link below. In poker1.mws, you can go down to "Putting it all together: repeated complex eigenvalues."

Day 13:

On Day 13, we began non-homogeneous linear differential equations, which can all be written as L[y] = rhs, with L a linear differential operator and rhs a non-zero right hand side.

For non-homogeneous linear differential equations, if we have any (one particular) non-homogenous solution, then we can get all solutions by adding all homogeneous solutions to the one non-homogeneous solution.

To see that adding homogeneous solutions to a non-homogeneous solution again gives a non-homogeneous solution, we note that if yp is a non-homogeneous solution (meaning L[yp] = rhs ) and yh is a homogeneous solution (meaning L[yh] = 0 ), then L[yp + yh] = L[yp] + L[yh] = rhs + 0 = rhs. Hence L[yp + yh] is again a non-homogeneous solution, since it makes the correct rhs.

Saying that yp + hy is again a nonhomogeneous solution does not, a priori, say that all nonhomogeneous solutions can be made from yp. However, that is also true.

To see that all non-homogeneous solutions can be made from any one particular non-homogeneous solution, note that if y is any non-homogenous solution and yp is the one particular non-homogeneous equation that we referred to earlier, then L[y] = rhs and L[yp] = rhs, so L[y - yp] = L[y] - L[yp] = rhs - rhs = 0. Hence, y - yp is a nonhomogeneous solution: y - yp = yh. But if y - yp = yh, then y = yp + yh, so every non-homogenous solution y can be written in terms of any fixed particular solution.

Note that the Wronskian determinant is necessary to make sure that we have all homogeneous solutions, but the Wronskian is not needed nor used with non-homogeneous solutions.

Since every non-homogeneous solution can be gotten from any one particular non-homogeneous solution by superposition (add one solution to another solution), the key element in the solution of a linear non-homogeneous differential equation is to find some particular non-homogeneous solution.

There are two ways to find a particular solution: one method only works some of the time (undetermined coefficients) but gives nice answers to many engineering problems, while the other method (variation of parameters) always works but gives answers that are not always so nice. We will learn and be responsible for both.

Undetermined Coefficients

The first of the two methods is called undetermined coefficients. It works only for constant coefficient differential equations, and then only if the rhs is in one of five possible forms (rhs = polynomial, rhs = polynomial times exponential, rhs = polynomial times sine or cosine, rhs = polynomial times exp times sine or cosine, or sums of these forms). Note that division is not allowed, and that the method will not work with a rhs that has a discontinuity.

There is some room for cleverness: a rhs of sin(x) cos(x) is not one of the five types, nor is 2^x, nor is sin(t)/exp(x). However, when they are rewritten as 1/2 sin(2 x), exp(a x) with a = ln(2), and as exp(-x) sin(x), respectively, then the forms match and the method of undetermined coefficients may be used. However no slight of hand nor trig identity will change a discontinuous rhs into one of the five forms, so a rhs of tan(x) can never be used for undetermined coefficients. We will need variation of parameters for that sort of rhs.

The first thing that one does to solve a non-homogeneous differential is to solve the homogeneous version. Always, this step first. Even if the problem asks you only to look at the form of the particular solution, you must solve the homogeneous first. No exceptions.

After solving the homogeneous equation, we then ask three questions to determine the form of a particular solution.

Q1: What degree polynomial?

Q2: Twins?

Q3: Any deaths?

The Method of Undetermined Coefficients

1) Solve the homogeneous equation.

What degree polynomial?

2) Determine the initial form of the particular solution by WYSIWYG: ask what degree polynomial. If the rhs has a polynomial of degree n, then the initial guess for the form of yp will have a polynomial of degree n with as many variables as we can fit in. (The only thing that is important in the form is the degree, so a rhs of -t^2 or of t^2 +1 or of 2*t^2 all get the same initial form: yp will have a*t^2+b*t+c.)

If you do not see a polynomial, ask if there is a hidden polynomial of degree zero, called 1.

If the rhs has exp(3 t) then the initial form has exp(3 t) multiplying the polynomial.

Twins?

4) If the rhs has either just sin(4 t) or just cos(4 t) or both, then the initial guess will have two terms, a polynomial times sin(4 t) plus a different polynomial times cos(4 t). Both twins are invited to the party, and they each get their own piece of polynomial birthday cake. (Of course, if the rhs has either exp(3 t) sin(4 t) or exp(3 t) cos(4 t) or both, then the initial guess contains two terms, a polynomial times exp(3 t) sin(4 t) plus a different polynomial times exp(3 t) cos(4 t).)

Deaths?

5) If any of the homogeneous solutions are exactly contained in the initial guess for yp, then the whole yp solution is multiplied by the smallest power of t so that no unknown coefficients are lost.

For example, with the differential equation y" - 2y' + y = t exp(t), which can be written as L[y] = t exp(t), with L[y] = y" - 2y' + y, then a fundamental set of homogeneous solutions is { exp(t) , t exp(t) }.

Since the rhs is a polynomial of degree one times exp(t), the initial guess is also a polynomial of degree one time exp(t): yp = (at + b)exp(t) (there are no twins).

By using the distributive law, we see that yp = a t*exp(t) + b exp(t). Note that the functions multiplying a and b, t*exp(t) and exp(t), are homogeneous solutions. (Seeing that they are homogeneous solutions is only possible *after* we see the homogeneous solutions. That's why we solve for them up front.) Hence, a modification of the initial form is necessary. If we subs and simplify, ie., apply the operator L to our initial guess, then we get

L[yp] = a L[t*exp(t)] + b L[exp(t)] = a*0 + b*0 = t*exp(t),

which cannot be solved for a and b because a and b are not available in the resulting equation 0 = rhs. The functions exp(t) and t exp(t) were annhilated by L, meaning L[exp(t)] = 0 and L[t exp(t)] = 0. That's OK, but the a and b were also killed when the functions that they multiplied were killed.

If any undetermined coefficient letters disappear, we need to multiply by one or more powers of t. Note that multiplying by a power of t changes the functions that multiply the undetermined coefficients, and there may no longer be an overlap with the homogeneous solutions.

If we only multiply the initial yp by one copy of t, then there will be an a in the equation to solve, but not a b: ie., if the modified guess is yp = (at + b)exp(t)*t = a t^2*exp(t) + b t*exp(t), then L[yp] = a L[t^2*exp(t)] + b L[t*exp(t)] = a L[t^2*exp(t)] + b*0, since t*exp(t) is a homogeneous solution and L[homogeneous solution] = 0.

Therefore, we must multiply by still one more copy of t, to get a modified form for a particular solution yp = (at + b)exp(t)*t^2. Now, neither a nor b is lost when we subs and simplify, so both variables will be there to solve for.

Yes, Virginia, one can get in trouble by multiplying by too high a power of t. We need the smallest power of t that assures no letters are lost.

A Multiple Term RHS

Note that modifications are done to the initial guess for a solution, not just a part of that solution; but only those solutions that need it. For example, if the differential equation has eigenvalues +/- i, so that the homogeneous solutions are sin(t) and cos(t), then a rhs = t + t*sin(t) would be broken into two parts, a rhs of t and a rhs of t*sin(t). The initial guesses for the two parts are made separately.

The initial guess for the rhs of t is yp = at + b, and that initial guess will need no modification since no homogeneous solutions are involved (no sin or cos in at + b) and no letters are lost when we subs and simplify.

On the other hand, the rhs of t*sin(t) has initial guess of (ct + d)sin(t) + (et + f)cos(t), because the polynomial is of degree one, and there are twins, both of which must be included, each with its own polynomial. Moreover, the t*sin(t) initial guess must be modified, since two letters will be lost in passing through L:

L[c t*sin(t) + d sin(t) + e t*cos(t) + fcos(t)] =
c L[t*sin(t) + 0 + e L[t*cos(t) + 0.

Hence d and f are not there to solve for.

In particular, the correct final form form for yp is then yp = at + b + ((ct + d)sin(t) + (et + f)cos(t))*t. The first term needs no adjustment, while the second term needs only one copy of t, since the functions in ((ct + d)sin(t) + (et + f)cos(t))*t are t^2*sin(t), t*sin(t), t^2*cos(t), and t*cos(t), none of which are homogeneous solutions.

Once the form is adjusted by the proper power of t, then we subs and simplify the particular solution to solve for the undetermined coefficients.

Here is an example with all the Maple code.

A nonhomogeneous constant coefficient linear equation:

(Recall that the method only applies if the equation is constant 
coefficient and the right hand side matches the allowable kinds of rhs.)

> de:=diff(y(t),t$2)+y(t)=t+t*sin(t);

                      / 2      \
                      |d       |
                de := |--- y(t)| + y(t) = t + t sin(t)
                      |  2     |
                      \dt      /

Solve the homogeneous solution.
> m^2+1=0;
> solve(%,{m});

                               2
                              m  + 1 = 0


                          {m = I}, {m = -I}

The check your homogeneous solution
> solh:=y(t)=c1*sin(t)+c2*cos(t);
> subs(solh,lhs(de)):
> simplify(%);

                 solh := y(t) = c1 sin(t) + c2 cos(t)


                                  0

Check fundamental.
> with(VectorCalculus):
> with(LinearAlgebra):
> FS:=[sin(t),cos(t)]:
> W:=Wronskian(FS,t);
> Determinant(%):
> simplify(%);
Warning, the names &x, CrossProduct and DotProduct have been rebound

Warning, the names &x, CrossProduct and DotProduct have been rebound


                            [sin(t)    cos(t) ]
                       W := [                 ]
                            [cos(t)    -sin(t)]


                                  -1

After getting all the homogeneous solutions, we are finally ready to look 
for a (NOT "the") particular solution.

The two rhs terms are evaluated separately. 

The first term is a first degree polynomial, no twins, so our guess is 
at + b. No letters are annihilated, so no extra t's are needed to prevent 
loss of coefficients.

The second term has a polynomial of degree one and twins. The initial 
guess is  (ct + d) sin(t) + (et + f) cos(t), which we can multiply out as 
c t sin(t) + d sin(t) + e t cos(t) + f cos(t). The coefficents of d and f 
are the functions sin(t) and cos(t), which are homogeneous solutions. 
Hence, when we subs and simplify (apply the linear opeator L), then we 
will have no d and f to solve for. Therefore, we
need to multiply the entire expression by t. (Note that after we multiply 
by t, the coefficient functions of d and f are no longer sin(t)
and cos(t), they are t sin(t) and t cos(t) ... and those two functions are 
*not* homogensous solutions.)

Therefore the final guess is (WARNING: use lower case letters. D has a 
special meaning in Maple!)


> solp:=y(t)=a*t+b+t*((c*t+d)*sin(t)+(e*t+f)*cos(t));

   solp := y(t) = a t + b + t ((c t + d) sin(t) + (e t + f) cos(t))

Subs and simplify the form of the particular solution.
> subs(solp,de);
> simplify(%);

  / 2                                                     \
  |d                                                      |
  |--- (a t + b + t ((c t + d) sin(t) + (e t + f) cos(t)))| + a t + b
  |  2                                                    |
  \dt                                                     /

         + t ((c t + d) sin(t) + (e t + f) cos(t)) = t + t sin(t)


  2 c sin(t) + 4 cos(t) c t + 2 cos(t) d + 2 e cos(t) - 4 sin(t) e t

         - 2 sin(t) f + a t + b = t (1 + sin(t))

The six equations to solve for a, b, c, d, e and f are

Match coefficient of sin(t):           2c = 0   
Match coefficient of t sin(t):        -4e = 1
Match coefficient of cos(t):      2d + 2e = 0
Match coefficient of t cos(t):         4c = 0 
Match coefficients of t:                a = 1
Match constants:                        b = 0

You should be able to set up and solve those equations, but Maple will 
even solve for us if we make this into an identity.

> identity(%,t);
> abcdefsol:=solve(%,{a,b,c,d,e,f});

  identity(2 c sin(t) + 4 cos(t) c t + 2 cos(t) d + 2 e cos(t)

         - 4 sin(t) e t - 2 sin(t) f + a t + b = t (1 + sin(t)), t)


     abcdefsol := {a = 1, b = 0, c = 0, f = 0, e = -1/4, d = 1/4}

Now that we have determined the undetermined coefficients, a, b, c, d, e, 
and f, we get a particular homogeneous solution.
> subs(abcdefsol,solp);

               y(t) = t + t (1/4 sin(t) - 1/4 t cos(t))

To get the general solution, we add all the homogeneous solutions. (The 
last object is an equation; and we cannot add a number to an equation, so 
we hum a few bars and fake an equation to add.)

The general solution is
> gsol:=%+(0=rhs(hsol));

  gsol := y(t) =

        t + t (1/4 sin(t) - 1/4 t cos(t)) + c1 sin(t) + c2 cos(t)

WARNING: A common error with nonhomogeneous equation initial values is to try to use the initial values too soon. You only work with the initial values *after* you have the nonhomogeneous general solution. (Students sometimes try to apply the initial values after the first step, when all they have is the homogeneous general solution. That is wrong. The initial values apply to the nonhomogeneous solution, not the homogeneous solution; and if you apply them too soo, then the part due to a particular solution is not included in the calculations.)

Day 14:

On Day 14 we reviewed the Method of Undetermined Coefficients for finding a first (particular) nonhomogeneous solution.

The method applies to constant coefficient equations for which the right hand side can be written in the form p(x) exp(ax) sin(bx) or p(x) exp(ax) cos(bx). There will always be a polynomial, but the exponential and trig functions need not appear.

The first step in finding the form of the first nonhomogeneous polynomial is to find the DEGREE of the existing polynomial, since the existing polynomial will be replaced by a polynomial having the same degree, but with as many letters as we can fit into a polynomial with that degree.

The second step is to make sure that if the original rhs contains either sine or cosine, then both sine and cosine are in the form of the particular solution, each with a different polynomial of the correct degree.

The third step is to check that no letters will be lost, since otherwise, they will not be in the equations that we use to solve for the undetermined coefficients. We will be unable to solve for a coefficient in the form if the function that multiplies that coefficient is annihilated by the linear operator L[y] from the nonhomogeneous equation, L[y] = 0. In other words, we will be unable to solve for a coefficient in the form if the function that multiplies that coefficient is a homogeneous solution.

After one has written an initial form for a particular solution, using the degree of the polynomial and twins, one will need to multiply the entire form by the smallest power of x so that none of the coefficients are multiplied by homogeneous solutions.

If there is a repetition of the same function in the form, we omit the second occurrence.

Here is an example:

Solve the nonhomogeneous constant coefficient differential equation L[y] = 
sin(x) - x cos(x), where L[y] = y" + y.

> de:=diff(y(x),x$2)+y(x)= sin(x)-x*cos(x);

                   / 2      \
                   |d       |
             de := |--- y(x)| + y(x) = sin(x) - x cos(x)
                   |  2     |
                   \dx      /

> m^2+1=0;
> solve(%,{m});

                               2
                              m  + 1 = 0


                          {m = I}, {m = -I}

A fundamental solution set is FS = { sin(x) , cos(x) }.

Check that sin(x) and cos(x)  are solutions
> hsol:=y(x)=c1*sin(x)+c2*cos(x);
> lhs(de)=0;
> subs(hsol,%):
> simplify(%);

                 hsol := y(x) = c1 sin(x) + c2 cos(x)


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


                                0 = 0

Check that sin(x) and cos(x) are linearly independent
> with(VectorCalculus):
> with(LinearAlgebra):
> FS:=[sin(x),cos(x)];
> W:=Wronskian(FS,x);
> Determinant(%):
> simplify(%);
Warning, the names &x, CrossProduct and DotProduct have been rebound

Warning, the names &x, CrossProduct and DotProduct have been rebound


                        FS := [sin(x), cos(x)]


                            [sin(x)    cos(x) ]
                       W := [                 ]
                            [cos(x)    -sin(x)]


                                  -1

Now that we have completed the general solution to the associated 
homogeneous equation (and only then), we then look at the right hand side 
of the nonhomogeneous differential equation.

The polynomial that multiplies the sine is 1, which is a polynomial of 
degree 0. Hence, the polynomial 1 that multiplies sine will be replaced by 
the polynomial a.  Since both sin(x) and cos(x) must be included, we get 
FROM JUST THE SINE TERM the form 

y(x) = a sin(x) + b cos(x).

Since letters a and b will be lost when we subs and simplify (because the 
functions that multiply them are homogeneous solutions), we must
modify the form that comes from the sin(x) to get  

y(x) = x ( a sin(x) + b cos(x) ) = a x sin(x) + b x cos(x)

The polynomial that multiplies the cosine is x, which is a polynomial of 
degree 1. Hence, the polynomial  x that multiplies the cosine will be
replaced by the polynomial  cx + d, also of degree one but containing as
many letters as we can put in a degree one polynomial.
 
Since both sin(x) and cos(x) must be included, we get FROM JUST THE COSINE TERM
the form 

y(x) = (cx + d) sin(x) + (ex + f) cos(x) =
       c x sin(x) + d sin(x) + e x sin(x) + f cos(x).

Since the letters d and f would be lost when we subs and simplify 
(because the functions that multiply them are homogeneous solutions), we 
must modify the form that comes from the cosine to get 

y(x) = x ( (cx + d) sin(x) + (ex + f) cos(x) ) =
       c x^2 sin(x) + d x sin(x) + e x^2 cos(x) + f x cos(x).

Finally, we need to look at interactions between the forms coming from 
sin(x) and x cos(x), respectively.

Since a and d both multiply x sin(x) in the forms, we omit the duplicate, 
effectively choosing d to be zero. Ditto, since both b and f multiply x cos(x),
we omit the duplicate f. 

Thus, the final form (with duplicates d and f omitted) is

> solp:=y(x)=a*x*sin(x)+b*x*cos(x)+c*x^2*sin(x)+e*x^2*cos(x);

                                              2             2
  solp := y(x) = a x sin(x) + b x cos(x) + c x  sin(x) + e x  cos(x)

Apply the linear operator L using subs and simplify.

> subs(solp,de):
> simplify(%);

  2 a cos(x) - 2 b sin(x) + 2 c sin(x) + 4 c x cos(x) + 2 e cos(x)

         - 4 e x sin(x) = sin(x) - x cos(x)

Let Maple do the pattern matching with the identity command, then 
determine the undetermined coefficients.

> identity(%,x);
> abcesol:=solve(%,{a,b,c,e});

  identity(2 a cos(x) - 2 b sin(x) + 2 c sin(x) + 4 c x cos(x)

         + 2 e cos(x) - 4 e x sin(x) = sin(x) - x cos(x), x)


            abcesol := {e = 0, a = 0, c = -1/4, b = -3/4}

Find a particular nonhomogeneous solution

> subs(abcesol,solp);

                                             2
                 y(x) = -3/4 x cos(x) - 1/4 x  sin(x)

To get the general nonhomogeneous solution, we would then add all 
homogeneous solutions.

> gnhsol:=%+(0=rhs(hsol));

  gnhsol :=

                                    2
        y(x) = -3/4 x cos(x) - 1/4 x  sin(x) + c1 sin(x) + c2 cos(x)

Alternatively, if in solving for a particular solution you do not want to 
consider the interaction between the forms for sin(x) and for - x cos(x), 
you can solve separately for a particular solution to match just sin(x) and a
particular solution to match just - x cos(x), determining the coefficients 
separately in different calculations.

Using the form y(x) = x ( a sin(x) + b cos(x) ) on the L[y] = sin(x) 
problem, we get a particular solution y(x) = - 1/2 x cos(x). 

Re-solving, using the form   y(x) = (cx + d) sin(x) + (ex +f) cos(x) on 
the L[y] = - x cos(x) problem, we get y(x) =  - 1/4 x (cos(x) + x sin(x)). 

Superimposing the two particular solutions, we get a particular solution
for the sum of the separate right hand sides.

Since L[ -1/2 x cos(x) ] = sin(x) and 

L[ - 1/4 x (cos(x) + x sin(x)) ] = - x cos(x), 

then by linearity (superposition)
 
L[ -1/2 x cos(x) - 1/4 x (cos(x) + x sin(x)) ] = 
L[ -1/2 x sin(x)] + L[ - 1/4 x (cos(x) + x sin(x)) ] = 
sin(x) - x cos(x). 

Day 15:

On Day 15 we learned how to use the method of variation of parameters to solve for a particular solution.

As before, we will begin by solving for the homogeneous solutions. However, we no longer make a guess for the form of a particular solution. The particular solution will always have the form y(x) = v1(x) y1(x) + v2(x) y2(x), where y1(x) and y2(x) form a fundamental set of homogeneous solutions. Your assignment will be to determine what v1(x) and v2(x) can be.

Notice that the form of a particular nonhomogeneous solution is the same as the form of a homogeneous solution, except that the parameters c1 and c2 are no longer constant: they are allowed to be functions of x and hence vary. That's why this method is called variation of parameters.

We will find v1'(x) and v2'(x) from two simultaneous linear equations:

v1'(x) y1(x) + v2'(x) y2(x) = 0
v1'(x) y1'(x) + v2'(x) y2'(x) = rhs in standard form

Note that the coefficients of the derivatives of v form the Wronskian matrix,

 
[ y1(x)  y2(x) ]
[ y1'(x) y2'(x)]

and the right hand side of the system is all zero's down to the rhs of the (standard form) differential equation.)

(If we had a third order equation, there would be two zeros and the rhs of the standard form differential equation.)

To make sure that you understand the requirement that we use the rhs from the standard form equation, compare p. 195, (9) where the rhs is

[ 0  ]
[ g/a]

although the original equation is ay" + by' + cy = g(x).

Once we have solved the system of equations for v1'(x) and v2'(x), we will integrate to find v1(x) and v2(x), taking the dot product of [v1(x),v2(x)] and FS = [y1(x),y2(x)] to get a particular solution y(x) = v1(x) y1(x) + v2(x) y2(x).

We did one example of variation of parameters using Cramer's Rule, learned on Day 10, to solve for v1'(x) and v2'(x). However, it is better to use the matrix inverse in Maple to solve such systems with less typing.

You should download the Maple worksheet entitled variation of parameters and study it.

To give you two other examples to study, here is the matrix method again.

There are no restrictions on the kind of linear differential equation nor 
on the kind of right hand side when we use variation of parameters. Here 
is an example:

Use the method of variation of parameters to find a particular 
nonhomogeneous solution, then solve  the initial value problem given by 
y(1) = 1, y'(1) = 2, y"(1) = 3.

The equation is of type Cauchy-Euler.
> de:=x^3*diff(y(x),x$3)+3*x^2*diff(y(x),x$2)+2*x*diff(y(x),x)=x;

               / 3      \        / 2      \
             3 |d       |      2 |d       |       /d      \
      de := x  |--- y(x)| + 3 x  |--- y(x)| + 2 x |-- y(x)| = x
               |  3     |        |  2     |       \dx     /
               \dx      /        \dx      /

We find the indicial equation.
> sol:=y(x)=x^m;
> subs(sol,lhs(de)=0):
> simplify(%):
> %/x^m:
> expand(%);

                                          m
                           sol := y(x) = x


                               3
                              m  + m = 0

> solve(%,{m});

                      {m = 0}, {m = I}, {m = -I}

The general homogeneous solution.
> solh:=y(x)=c1+c2*sin(ln(x))+c3*cos(ln(x));

          solh := y(x) = c1 + c2 sin(ln(x)) + c3 cos(ln(x))

Verify solution.
> subs(solh,lhs(de)=0):
> simplify(%);

                                0 = 0

Verify linearly independent.
> with(VectorCalculus):
> with(LinearAlgebra):
> FS:=[1,sin(ln(x)),cos(ln(x))];
Warning, the names &x, CrossProduct and DotProduct have been rebound

Warning, the names &x, CrossProduct and DotProduct have been rebound


                  FS := [1, sin(ln(x)), cos(ln(x))]

> W:=Wronskian(FS,x);
> Determinant(%):
> simplify(%);

  W :=

        [1 , sin(ln(x)) , cos(ln(x))]

        [    cos(ln(x))     sin(ln(x))]
        [0 , ---------- , - ----------]
        [        x              x     ]

        [      sin(ln(x))   cos(ln(x))     cos(ln(x))   sin(ln(x))]
        [0 , - ---------- - ---------- , - ---------- + ----------]
        [           2            2              2            2    ]
        [          x            x              x            x     ]


                                   1
                                - ----
                                    3
                                   x

Variation of parameters. 

The system of equation to solve for v1', v2', and v3' are 
W Vp = R, where W has already been computed, and R is the
matrix below. (Note that we need the rhs of the standard form differential 
equation, so we divide x by the x^3 coefficient of y'''.)
> R:=Matrix(<0,0,x/x^3>);

                                  [ 0  ]
                                  [    ]
                                  [ 0  ]
                             R := [    ]
                                  [ 1  ]
                                  [----]
                                  [  2 ]
                                  [ x  ]

If  W Vp = R, then  Vp = 1/W R. Hence v1' is the top entry, v2' is the 
middle entry, and v3' is the bottom entry in the matrix.
> 1/W.R;
> simplify(%);

                    [             1             ]
                    [                           ]
                    [         sin(ln(x))        ]
                    [- -------------------------]
                    [            2             2]
                    [  cos(ln(x))  + sin(ln(x)) ]
                    [                           ]
                    [         cos(ln(x))        ]
                    [- -------------------------]
                    [            2             2]
                    [  cos(ln(x))  + sin(ln(x)) ]


                            [     1     ]
                            [           ]
                            [-sin(ln(x))]
                            [           ]
                            [-cos(ln(x))]

We want v1, v2, and v3, so we integrate everything in the matrix with one 
Maple command. (map does whatever the operation is to the entire box. What 
we are doing to everything is integrating dx.) v1 is the top entry, v2 is 
the middle, and v3 is the bottom entry.
> V:=map(int,%,x);

                  [                 x                  ]
                  [                                    ]
             V := [1/2 cos(ln(x)) x - 1/2 sin(ln(x)) x ]
                  [                                    ]
                  [-1/2 cos(ln(x)) x - 1/2 sin(ln(x)) x]

A particular solution is  y = v1 1 + v2 sin(ln(x) + v3 cos(ln(x)), which 
we can compute as a dot product of [v1,v2,v3] and 
[1,sin(ln(x)),cos(ln(x))]. 

Hence, we convert the vertical matrix V into a horizontal vector, do the same 
for the list FS, and take the dot product using 
(what else?!) a dot.
> convert(%,Vector).convert(FS,Vector);
> solp:=y(x)=simplify(%);

  x + (1/2 cos(ln(x)) x - 1/2 sin(ln(x)) x) sin(ln(x))

         + (-1/2 cos(ln(x)) x - 1/2 sin(ln(x)) x) cos(ln(x))


                          solp := y(x) = x/2

To the particular nonhomogeneous solution, we now add all the homogeneous 
solutions, to get the general nonhomogeneous solution.

> nhgsol:=%+(0=c1+c2*sin(ln(x))+c3*cos(ln(x)));

      nhgsol := y(x) = x/2 + c1 + c2 sin(ln(x)) + c3 cos(ln(x))

Now we solve for c1, c2, and c3 from the initial values.
 
y(1) = 1.
> 1=subs(x=1,rhs(nhgsol)):
> eq1:=simplify(%);

                       eq1 := 1 = 1/2 + c1 + c3

y'(1) = 2.
> diff(nhgsol,x):
> 2=subs(x=1,rhs(%)):
> eq2:=simplify(%);

                         eq2 := 2 = 1/2 + c2

y"(1) = 3
> diff(nhgsol,x$2):
> 3=subs(x=1,rhs(%)):
> eq3:=simplify(%);

                         eq3 := 3 = -c2 - c3

Solve for c1, c2, and c3.
> csol:=solve({eq1,eq2,eq3},{c1,c2,c3});

                csol := {c1 = 5, c3 = -9/2, c2 = 3/2}

Hence, the unique solution to the nonhomogeneous Cauchy-Euler initial 
value problem is

> ivpsol:=subs(csol,nhgsol);

      ivpsol := y(x) = x/2 + 5 + 3/2 sin(ln(x)) - 9/2 cos(ln(x))

To check to see if this is a solution to the nonhomogeneous differential 
equation, we subs and simplify. You can use subs and diff to check that 
the initial values match, but we leave that to the reader.

> subs(ivpsol,de):
> simplify(%);

                                x = x

p. 198 #25 shows that we can apply even apply variation of parameters to a 
nonhomogeneous Bessel equation, provided someone gives us the homogeneous 
solutions.

Note that the equation is not constant coefficient, since x^2 is not a 
constant. Note that the equation is not Cauchy-Euler, since the 
coefficient of y(x) could only be a constant if this were Cauchy-Euler, 
not (x^2 - 1/4).

Since it is neither constant coefficient nor Cauchy-Euler, you are not 
going to be able to find a fundamental set of homogeneous solutions.
But then again, no one asked you to. They gave you the 
solutions. Your job is 1) check that the solutions are valid and 2) find a 
particular solution, and 3) add all homogeneous solutions.

> de:=x^2*diff(y(x),x$2)+x*diff(y(x),x)+(x^2-1/4)*y(x)=x^(5/2);

              / 2      \
            2 |d       |     /d      \     2                (5/2)
     de := x  |--- y(x)| + x |-- y(x)| + (x  - 1/4) y(x) = x
              |  2     |     \dx     /
              \dx      /

Check that the putative solutions given in the book really work.

> sol:=y(x)=c1*x^(-1/2)*cos(x)+c2*x^(-1/2)*sin(x);

                               c1 cos(x)   c2 sin(x)
                 sol := y(x) = --------- + ---------
                                  1/2         1/2
                                 x           x

Verify solutions.

> subs(sol,lhs(de)=0):
> simplify(%);

                                0 = 0

Verify linear independence.

> with(VectorCalculus):
> with(LinearAlgebra):
> FS:=[x^(-1/2)*cos(x),x^(-1/2)*sin(x)];
> W:=Wronskian(%,x);
> Determinant(%):
> simplify(%);
Warning, the names &x, CrossProduct and DotProduct have been rebound

Warning, the names &x, CrossProduct and DotProduct have been rebound


                               cos(x)  sin(x)
                        FS := [------, ------]
                                 1/2     1/2
                                x       x


              [       cos(x)                  sin(x)       ]
              [       ------                  ------       ]
              [         1/2                     1/2        ]
              [        x                       x           ]
         W := [                                            ]
              [     cos(x)   sin(x)         sin(x)   cos(x)]
              [-1/2 ------ - ------    -1/2 ------ + ------]
              [       3/2      1/2            3/2      1/2 ]
              [      x        x              x        x    ]


                                 1/x

Variation of parameters. 

The equations to solve for v1' and v2' are W Vp = R, where W has already been 
computed, and R is the matrix below.
 
(Note that we need the rhs of the standard form differential 
equation, so we divide x^(5/2) by the x^2 coefficient of y".)

> R:=Matrix(<0,x^(5/2)/x^2>);

                                  [ 0  ]
                             R := [    ]
                                  [ 1/2]
                                  [x   ]

If  W Vp = R, then  Vp = 1/W R. Hence v1' is the top entry and v2' is the 
bottom entry in the matrix.

> 1/W.R;
> simplify(%);

                        [      x sin(x)     ]
                        [- -----------------]
                        [        2         2]
                        [  cos(x)  + sin(x) ]
                        [                   ]
                        [     x cos(x)      ]
                        [ ----------------- ]
                        [       2         2 ]
                        [ cos(x)  + sin(x)  ]


                             [-sin(x) x]
                             [         ]
                             [cos(x) x ]

We want v1 and v2,  so we integrate everything in the matrix with one 
Maple command. (map does whatever the operation is to the entire box. What 
we are doing to everything is integrating dx.) 

v1 is the top entry and v2 is the bottom entry.

> V:=map(int,%,x);

                           [-sin(x) + cos(x) x]
                      V := [                  ]
                           [cos(x) + sin(x) x ]

A particular solution is  y = v1 x^(-1/2)*cos(x) + v2 *x^(-1/2)*sin(x), 
which we can compute as a dot product of [v1,v2] and 
[x^(-1/2*cos(x),x^(-1/2)*sin(x)]. Hence, we convert the vertical matrix V 
into a horizontal vector, do the same for FS, and take the dot product 
using (what else?!) a dot.

A particular solution is 

> convert(%,Vector).convert(FS,Vector);
> solp:=y(x)=simplify(%);

       (-sin(x) + cos(x) x) cos(x)   (cos(x) + sin(x) x) sin(x)
       --------------------------- + --------------------------
                   1/2                           1/2
                  x                             x


                                         1/2
                         solp := y(x) = x

And we add all the homogeneous solutions to get the general nonhomogeneous 
solution.

> nhgsol:=solp+(0=c1*x^(-1/2)*cos(x)+c2*x^(-1/2)*sin(x));

                              1/2   c1 cos(x)   c2 sin(x)
            nhgsol := y(x) = x    + --------- + ---------
                                       1/2         1/2
                                      x           x

Pick an easy one and find the v-primes by Cramer's rule, then do the integration with Maple. Next, repeat that problem using the technique in the variation of parameters worksheet, checking that you get the same v-primes, same v's, and the same particular solution.

Then do all other problems using only the matrix multiplication technique.

Day 16:

On Day 16 we reviewed variation of parameters, especially the implementation of this technique using the efficient matrix inversion technique available in Maple. We also discussed briefly how second order linear constant coefficient differential equations arise in mechanical engineering in spring-mass-dashpot systems and in electrical engineering in L-R-C series circuits.

A series electrical circuit gives a second order linear constant coefficient equation L y" + R y' + 1/C y = f(t). A mass-dashpot-spring system gives a second order linear constant coefficient equation m y" + b y' + k y = f(t). Both situations frequently occur with inputs which are polynomials times exponentials times sines or cosines. Hence the method of undetermined coefficients is frequently used by engineers. (Of course, as always, we will write down the appropriate differential equation on applied problems, let Maple solve it, and work with the solution that she gives.)

Since the characteristic equation in both cases is a quadratic, and since a quadratic with real coefficients has only three kinds of answers, based on the discriminant b^2-4ac, the electrical and spring systems have only three kinds of behavior. If b^2-4ac>0, then we have two distinct real eigenvalues; whereas if b^2-4ac=0, then there are repeated real eigenvalues; and if b^2-4ac<0, then there are complex eigenvalues, with sines and cosines in the homogeneous solutions. Those are the only kinds of solutions that can occur.

The first situation is described as an overdamped system, and is very slow to return to equilibrium. In the third situation, the underdamped case, the system oscillates forever under an exponential curve that goes slowly to zero. The optimal arrangement is the middle case, called critical damping, in which the situation returns to equilibrium as quickly as possible, with at most one crossing of the equilibrium position before it gradually settles back down.

Objects that oscillate have natural frequencies that they prefer. If one excites such systems with a signal, even of small amplitude, which matches the preferred frequency, the system will begin to oscillate at high amplitudes. This phenomenon is called resonance. For both mechanical and electrical systems, if the oscillation builds up to high enough amplitudes, a mechanical or electrical overload can occur, resulting in the destruction of the system.

One of the things mentioned in the chapter is a repeat of something that you probably learned in trigonometry. In the underdamped situation, a solution to a homogeneous equation is is a linear combination c1*sin(wt)+c2*cos(wt). However, for fixed c1 and c2, the linear combination of a sine and cosine can be converted into a linear combination of a single sine (or cosine, if you like) by a silly trick:

c1*sin(wt)+c2*cos(wt) =
sqrt(c1^2 + c2^2) ( c1/(sqrt(c1^2 + c2^2)*sin(wt) + c2/sqrt(c1^2 + c2^2)*cos(wt) ),
which becomes cos(delta)*sin(wt) + sin(delta)*cos(wt),
which is A sin(wt+delta), by the sum of angles formula for sine, where the amplitude A is sqrt(c1^2 + c2^2), and delta = arccos(c1/(sqrt(c1^2 + c2^2)).
The value of delta is called the phase shift.

We will skip Chapter 5 in its entirety, since Chapter 5 is a baby version of the material in Chapter 9 (which we will cover after Chapter 7, the Laplace transform).

Day 17:

On Day 17 we began to study the use of Laplace transforms to solve constant coefficient initial value problems.

Laplace transforms are a very useful tool to solve constant coefficient differential equations, since both homogeneous and non-homogeneous equations are handled the same way. Laplace is also an excellent tool to handle such differential equations when they have discontinuous inputs. Since abrupt changes in the input to electrical and mechanical systems are often encountered in applications, Laplace methods are often used by mechanical, electrical, and chemical engineers.

The Laplace transform of a "real" function f(t) is a function F(s). (Functions of s are always Laplace transforms here, functions of t are always "real", and the uppercase-lowercase same letter convention is often used to show which functions are related in a transform pair.).

The Laplace transform F(s) of a function f(t) is given by multiplying f(t) by the Laplace kernel exp(-st) and integrating with respect to t from t = 0 to t = infinity. As such, a Laplace transform is an improper integral, which may converge or may diverge, depending on the size of s.

Hidden in every Laplace transform is a choice of s to make the integral converge. This fact is indicated in NSS by saying that f(t) is "of exponential order". What that means in Texian is "Is there a constant s that will make the improper integral converge?" The Laplace transform is Int(f(t)/exp(st),t=0..infinity). If the integral is not to pick up an infinite amount of area as we integrate out to infinity, then the integrand f(t)/exp(st) must quickly get very close to zero. For example, if we are trying to take the Laplace transform of exp(49t), then we need exp(49t)/exp(st) to go quickly to zero. How big must we choose s? Anything bigger than 49 will do, so exp(49t) is of exponential order, with s > 49. On the other hand, exp(t^2) is not of exponential order, since no matter how big we choose the constant s, exp(t^2)/exp(st) does not go to zero as t gets large.

There are four ways to get the Laplace transform of a function f(t).

• 1) You can integrate by hand, with the appropriate assumption on s.

• 2) You can integrate by Maple, with the appropriate assumption on s. (For example, to find the Laplace transform of exp(49t) with Maple doing the integration, you must say assume(s > 49); for the limit to be taken. The default for Maple is to put a tilde with a variable for which some assumption has been made. You can and should turn off the tilde, which is ugly and causes problems, by saying interface(showassumed=0); which turns off the tilde.) You can do this globally by going into Tools > Options > Display > Assumed Variables > No annotation.

  One can go backwards from F(s) to f(t) using contour integrals, but that techniques is only taught in Math 407, Complex Variables.

• 3) You can look up the Laplace transform in a table of Laplace transforms, a copy of which is in the inside front cover of your text. Like Highway 6, the table runs both ways: you can start with an f(t) and get F(s), or you can start with an F(s) and go to f(t).

• 4) You can use Maple's Laplace command, in the integral transforms package, after the command with(inttrans); The format is laplace(f(t),t,s); which gives the Laplce transform of a function of t as a function of s. Maple also has a command to go backwards form F(s) to f(t), called invlaplace(F(s),s,t);

There are six properties of the Laplace transform.

• 1) The Laplace transform is a linear operator, since integration is a linear operator.

• 2) Differentiation by t corresponds to multiplication by s (with f(0) subtracted off):
  L( d/dt f) = s L(f) - f(0),
  ie., cover up the derivative, take the Laplace transform of the original function, multiply by s, then subtract off the value of the original function at 0.

  One can apply this rule twice to get the Laplace transform of higher derivatives with respect to t: L( d/dt d/dt f) = s L( d/dt f) - f'(0) = s ( s L(f) - f(0) ) - f'(0).

  The property follows by integration by parts:
  L( d/dt(f) ) = int( exp(-st) d/dt(f) , t=0..infinity) = exp(-st) f(t) at t=infinity - exp(-st) f(t) at t=0 - int( -s exp(-st) f,t=0..infinity)
  which is 0 - f(0) + s L(f), since f(t) is of exponential order, and exp(-st) f(t) goes to zero as t gets large.

  This Property 2 pairs with Property 3.

• 3) Multiplication by t corresponds to differentiation by s (with a minus sign):
  L( t f) = - d/ds L(f),
  ie., cover up the t, take the Laplace transform of the other function, take the derivative with respect to s, then add a minus.

  One can apply this rule more than once to get the Laplace transform of a function multipled by a higher power of t: L( t^2 f) = - d/ds L(t f) = d/ds d/ds L(f).

  This property comes from Leibniz's Integral Rule in Math 251.
  Leibniz's Rule
  L( t f) = int( exp(-st) t f, t=0..infinity) =
  int( -d/ds( exp(-st) ) f, t=0..infinity) =
  -d/ds int( exp(-st) f, t=0..infinity) =
  - d/ds L(f),
  where we used the Leibniz rule to reverse the order of the derivative and the integral.

• 4) Multiplication by exp(at) corresponds to translation right by a:
  L( exp(at) f ) = L(f)(s-a),
  ie., cover up the exp(at), take the Laplace transform of the other function, replace s by s-a.
  L(exp(at) f) = int( exp(-st) exp(at) f, t=0..infinity) = int( exp(-(s-a)) f, t = 0..infinity) = L(f)(s-a)./p>

  Property 4 pairs with another translation rule, Property 5.

• 5) Multiplication of exp(cs) corresponds to translation right by c in t-land, but we need the Heaviside function to complete this rule.

• 6) Multiplication of two Laplace transforms corresponds to a convolution product in t-land. but we need to know about the definition of the convolution product of two functions of t to make sense of this.

Of these, the most important property is the second, since it allows us to change a constant coefficient differential equation into a multiply-by-s equation.

Example: y'(t) = 2 y(t), y(0)=3, an initial value problem. We solve it in four steps:

• 1) Laplace the equation: s L(y) - y(0) = 2 L(y) (NOT a de!)
  

• 2) Subs the inits: s L(y) -3 = 2 L(y)
  

• 3) Solve for the Laplace of the answer: L(y) = 3/(s-2)
  

• 4) Unlaplace: y(t) = 3 exp(2t)

Remark: why does the Laplace method not work well for Cauchy-Euler equations? Note that by rule 3, the Laplace transform of a Cauchy-Euler equation is not a multiply by s equation, but instead, just a different differential equation, this time a differential equation in s.

Of course, in general, the Laplace of the answer is a complicated fraction, not the simple fractions that are found in the Laplace transform tables, so we need a technique from 8th grade algebra to convert complicated fractions to simple fractions. (This is the same problem that you encountered in integrating complicated fractions in Math 152, and the solution is the same: convert the complicated fraction is s to partial fractions.)

We will always use a Maple command to do this easy but computationally onerous calculation. The command is convert(%,parfrac,s);
After we have easy fractions, then we can go backwards through the table to find the answer y(t).

Caveat: convert(%,parfrac,s); only works if % is a rational function in s. If there is an exp(cs) in the object, you need to hold that part in abeyance until after using the parfrac.

The HW problems are in a Maple handout for which the link can be obtained below. You need to begin *memorizing* the short Laplace table on p. 358. (Note that there is a longer one on the inside front cover of NSS.) The next lecture will assume that you know the table.

Day 18:

On Day 18 we reviewed the properties of the Laplace transform and applied them to find the Laplace and inverse Laplace transforms of functions.

The four step procedure to solve constant coefficient initial value 
problems by use of the Laplace transform:
p. 383 #18
> de:=diff(y(t),t$2)-2*diff(y(t),t)-y(t)=exp(2*t)-exp(t);
> inits:=y(0)=1,D(y)(0)=3;

            / 2      \
            |d       |     /d      \
      de := |--- y(t)| - 2 |-- y(t)| - y(t) = exp(2 t) - exp(t)
            |  2     |     \dt     /
            \dt      /


                    inits := y(0) = 1, D(y)(0) = 3

Laplace the differential equation.
> with(inttrans):
> laplace(de,t,s);

   2
  s  laplace(y(t), t, s) - D(y)(0) - s y(0) - 2 s laplace(y(t), t, s)

                                                 1
         + 2 y(0) - laplace(y(t), t, s) = ---------------
                                          (s - 2) (s - 1)

Substitute the initial values.
> subs(inits,%);

   2
  s  laplace(y(t), t, s) - 1 - s - 2 s laplace(y(t), t, s)

                                        1
         - laplace(y(t), t, s) = ---------------
                                 (s - 2) (s - 1)

Solve for the Laplace transform of the solution.
> solve(%,laplace(y(t),t,s));

                              2            3
                          -2 s  - s + 3 + s
                       ------------------------
                        4      3      2
                       s  - 5 s  + 7 s  - s - 2

Find the inverse Laplace transform to get the solution.
> sol:=y(t)=invlaplace(%,s,t);

  sol:=y(t)=
                                    1/2       1/2         1/2
  -exp(2 t) + 1/2 exp(t) (3 cosh(t 2   ) + 3 2    sinh(t 2   ) + 1)

An example of how we might apply the properties of the Laplace transform to find a Laplace transform of a certain f(t) would be p. 365 #9, where we are asked to compute the transform of exp(-t)*t*sin(2*t), which is not in the table.

We recall the Laplace transform of sin(2t) from the table that we have memorized: 2/(s^2+4).

We note that multiplication by t corresponds to differentiation by s, yielding 4s/(s^2+4)^2.

We also note that multiplication by exp(-t) replaces s by s+1, to give 4*(s+1)/((s+1)^2+4)^2, which simplifies to 4*(s+1)/(s^2+2*s+5)^2.

We can check our hand work using Maple's Laplace command:

exp(-t)*t*sin(2*t);
> laplace(%,t,s);
                          exp(-t) t sin(2 t)


                              4 (s + 1)
                           ---------------
                             2           2
                           (s  + 2 s + 5)

We can also use the Laplace properties to find the inverse Laplace transform. For example, p. 374 #5, 1/(s^2+4s+8).

Note that the discriminant b^2-4ac is 4^2-32 = -16, so the quadratic is irreducible over the reals. Hence, we complete the square, to get
1/((s+2)^2+4).

The result reminds us of 2/(s^2+4), which is the Laplace transform of sin(2t). Hence, we write 1/2 sin(2t), allowing us to get the requisite 2 in the numerator. However, the s has been replaced by s+2, because of a multiplication by exp(2t), so the final inverse Laplace transform is 1/2 exp(2t) sin(2t).

We can, of course, check our hand work by Maple's inverse Laplace command:

> 1/(s^2+4*s+8);
> invlaplace(%,s,t);

                                  1
                             ------------
                              2
                             s  + 4 s + 8


                        1/2 exp(-2 t) sin(2 t)

If the quadratic is not irreducible and factors over the reals, then we can use partial fractions to reduce the complicated rational function to simple rational functions of the kind found in the Laplace table that you memorized.

For example, p. 375 #12, (-s-7)/(s^2-s-2).

Note that s^2-s-2 has a discriminant that is positive and a perfect square, so we can factor to get (s+1)(s-2).

We then use partial fractions.

(-s-7)/(s^2-s-2);
> convert(%,parfrac,s);
                                -s - 7
                              ----------
                               2
                              s  - s - 2


                               3       2
                           - ----- + -----
                             s - 2   s + 1

The inverse transform is then clear: -3 exp(2t) + 2 exp(-t). We can check our hand work using Maple's inverse Laplace command.

> (-s-7)/(s^2-s-2);
> invlaplace(%,s,t);
                                -s - 7
                              ----------
                               2
                              s  - s - 2


                        2 exp(-t) - 3 exp(2 t)

Note that when using Maple's integrate command to find the Laplace transform from the definition, one needs to indicate to Maple the assumption on s which is necessary to make the improper integral converge.

For example, p. 359 #5, finding the Laplace transform of exp(2t)cos(3t).

> exp(2*t)*cos(3*t);
> int(%*exp(-s*t),t=0..infinity);

       lim       - (exp(-t (s - 2)) cos(3 t) s - s
  t -> infinity

         - 2 exp(-t (s - 2)) cos(3 t) + 2

                                                    2
         - 3 exp(-t (s - 2)) sin(3 t))/(13 - 4 s + s )

Maple cannot compute the limit and finish the transform, since she does not know the value of s that makes exp(2*t)*cos(3*t) be of exponential order. The limit is left unevaluated.

To make this command work in Maple, you need to assume that s > 2, and tell Maple not to show tildes on the s.

> assume(s>2):
> interface(showassumed=0):
> exp(2*t)*cos(3*t);
> int(%*exp(-s*t),t=0..infinity);

                          exp(2 t) cos(3 t)


                               -2 + s
                            -------------
                                        2
                            13 - 4 s + s

Day 19:

On Day 19 we introduced what NSS calls the unit step function u[0](t) and what Maple calls the Heaviside function, Heaviside(t). Heaviside(t) is defined piecewise and is zero until t = 0, at which time it becomes 1.

Since replacing t by t-1 translates the graph of a function right by one, Heaviside(t-1) changes from being zero (asleep) to one (awake) at t = 1.

One can take any piecewise defined function and write it in terms of Heaviside functions using the (new - old)*Heaviside(t - transition) rule. For example, if we assume, as we always will, that f is zero until t = 0, then

> f:=piecewise(t<3,2,t<5,t,t<7,14-2*t,0);

                        {    2              t < 3
                        {
                        {    t              t < 5
                   f := {
                        { 14 - 2 t          t < 7
                        {
                        {    0            otherwise

> g:=(2-0)*Heaviside(t)+
     (t-2)*Heaviside(t-3)+
     ((14-2*t)-t)*Heaviside(t-5)+
     (0-(14-2*t))*Heaviside(t-7);

  g := 2 Heaviside(t) + (t - 2) Heaviside(-3 + t)

         + (14 - 3 t) Heaviside(-5 + t)

         + (-14 + 2 t) Heaviside(-7 + t)

There is no blue graph shown since it exactly overlaps the red one, 
showing that the two functions are equal.

> plot([f,g],t=0..8,discont=true,color=[red,blue]);

> with(inttrans):
> laplace(f,t,s);
> laplace(g,t,s);

            2 exp(-7 s) - exp(-5 s) (3 + s) + exp(-3 s) (s + 1)
      2/s + ---------------------------------------------------
                                     2
                                    s


            2 exp(-7 s) - exp(-5 s) (3 + s) + exp(-3 s) (s + 1)
      2/s + ---------------------------------------------------
                                     2
                                    s

For example, 1*Heaviside(t)+(t-1)*Heaviside(t-2)+(t^2-t)*Heaviside(t-5) has wakeup times at t = 0, t = 2, and t = 5. In the interval from 5 to infinity, all three Heavisides are awake and are replaced by 1. Hence, the sum for t > 5 is 1 + (t-1) + (t^2-t) = t^2.

Day 20:

Theorem 3 says that the Laplace of exp(at) times f(t) is the Laplace of f(t), which is F(s), but with the s replaced by s - a.

The dual to that is (See also HW p. 365 #31) on p. 387.

The dual, Theorem 8, says that the inverse Laplace of exp(-as) times F(s) is the inverse Laplace of F(s), which is f(t) Heaviside(t), but with the t replaced by t - a.

For example, the inverse Laplace transform of exp(-3s) 1/s^2 is the inverse Laplace of 1/s^2, which is (by the memorized table) t Heaviside(t), but with the t replaced throughout by t-3, eg., (t-3) Heaviside(t-3).

Note that it is always a good idea to include the Heaviside(t) when coming back from Laplace-land, even if no translation by t is contemplated. If you get in the habit of always placing the Heaviside(t) in the inverse Laplace, it will certainly be there if t is replaced by t - a.

The difference between f(t) and f(t) Heaviside(t) reveals what Heaviside "does." By multiplying by Heaviside(t), we force the value of the inverse Laplace to be zero before a certain time. This corresponds to the idea that an object in the laboratory should not begin to oscillate before someone enters the lab to start it oscillating.

We will never have to deal with the inverse Laplace transform of something like exp(3s) 1/s^2, since an exponential with a positive coefficient on s could make the object in the laboratory respond in a non-causal manner, *before* the force was applied to make the object oscillate. It will always be exp(-as).

A good example to work in Maple is the initial value problem y" + y = Heaviside(t-1), y(0) = 0, y'(0)=0, which is a mass-dashpot-spring system, as in p. 213 (12), with m = 1, b = 0, and k = 1, which is initially at the equilibrium position y = 0, and with no initial velocity when you enter the laboratory. As you will note from the graph, the mass remains still until after t = 1, when a 1 Newton force is applied and maintained forever, moving the equilibrium point up by one, about which point the mass then oscillates.

Here is the example:

de:=diff(y(t),t$2)+y(t)=Heaviside(t-1);
> inits:=y(0)=0,D(y)(0)=0;
> sol:=dsolve({de,inits},y(t));
                    / 2      \
                    |d       |
              de := |--- y(t)| + y(t) = Heaviside(t - 1)
                    |  2     |
                    \dt      /


                    inits := y(0) = 0, D(y)(0) = 0


          sol := y(t) = -Heaviside(t - 1) (-1 + cos(t - 1))

> plot(rhs(sol),t=0..10);

The purpose of the Heaviside(t-1) is to keep the (1 - cos(t-1)) from starting to oscillate before t = 1 when the force was first applied.

You should be aware of the Gamma function, mentioned on p. 393. The Maple equivalent is GAMMA(t), with all upper case letters.

As shown on p. 394, the Gamma function has the property that GAMMA(n + 1) = n!. Since the factorial is only defined at positive integers, we can use the GAMMA function to extend the definition to other positive real numbers, including those between integers.

For example, the Laplace transform of sqrt(t) is

> with(inttrans):
> laplace(sqrt(t),t,s);

                                  1/2
                                Pi
                                ------
                                   3/2
                                2 s

Now, according to the memorized table, the Laplace transform of t^n is 
n!/s^(n+1). Since sqrt(t) = t^(1/2), we are not surprised to find the
s^(3/2) in the denominator. On the other hand, we have a problem with the 
numerator, which should have (1/2)!, an undefined quantity.

On the other hand, if we extend the definition of factorial by the GAMMA 
function, (1/2)! should be GAMMA(3/2), which is exactly what we find.

> GAMMA(3/2);

                                  1/2
                                Pi
                                -----
                                  2

We also defined the convolution f(t)*g(t), which is another function of t. (WARNING! The asterisk is NOT Maple's multiplication symbol, it is the mathematical symbol for convolution.)

The Convolution Theorem, Theorem 11, p. 400, (different than the definition of convolution) relates the inverse Laplace of the ordinary product of two Laplace transforms to the convolution product of their inverse Laplace transforms.

To take the inverse Laplace transform of a product of two laplace transforms, we take the convolution product of each individual inverse Laplace transform.

Here is an example:

Taking the inverse Laplace transform using Maple's command,

> invlaplace(1/(s-1),s,t);
> f:=unapply(%,t);
> invlaplace(1/(s-2),s,t);
> g:=unapply(%,t);
                                exp(t)


                           f := t -> exp(t)


                               exp(2 t)


                          g := t -> exp(2 t)

then applying the Convolution Theorem, we get

> f(t-v)*g(v);
> Int(%,v=0..t);
> value(%);

                         exp(t - v) exp(2 v)


                        t
                       /
                      |
                      |   exp(t - v) exp(2 v) dv
                      |
                     /
                       0


                          exp(2 t) - exp(t)

And here is the same result using partial fractions:

> 1/((s-1)*(s-2));
> convert(%,parfrac,s);

                                  1
                           ---------------
                           (s - 1) (s - 2)


                              1       1
                            ----- - -----
                            s - 2   s - 1

which gives the same result if we apply the memorized table.

You should now be able to do all the assigned HW problems through those on p. 405.

For Thursday, we will cover Dirac(t), the derivative of Heaviside(t), which is not really a function, but a distribution. We will see that Dirac(t) is like a strobe light in a 1960's discotheque. When we integrate it against a function f(t) is freezes the value of f at the time when the strobe light flashes. For example,

> Int(exp(-2*t)*Dirac(t-1),t=0..infinity);

                   infinity
                  /
                 |
                 |          exp(-2 t) Dirac(t - 1) dt
                 |
                /
                  0

> value(%);

                               exp(-2)

The Dirac(t) is zero, so the product is zero, except at the single point 
t = 1, when the light strobes exp(-2t), whose value at that instant is 
exp(-2). There is no need to do an antiderivative.

A useful example is the initial value problem y" + y = Dirac(t-1), y(0) = 0, y'(0) = 0. This represents a Chinese gong, of mass m = 1, b = 0, k = 1, which is at rest in its equilibrium posiition at t = 0, when you enter the kung-fu movie. At time t = 1, the gong is given a short, sharp hammer blow (the Dirac impulse), at which time it starts to oscillate. Note the Heaviside(t-1) in the answer, which turns off the start of the oscillation until after the hammer blow.

de:=diff(y(t),t$2)+y(t)=Dirac(t-1);
> inits:=y(0)=0,D(y)(0)=0;
                      / 2      \
                      |d       |
                de := |--- y(t)| + y(t) = Dirac(t - 1)
                      |  2     |
                      \dt      /


                    inits := y(0) = 0, D(y)(0) = 0

The four step Laplace method
> laplace(de,t,s);
> subs(inits,%);
> solve(%,laplace(y(t),t,s));
> sol:=y(t)=invlaplace(%,s,t);

   2
  s  laplace(y(t), t, s) - D(y)(0) - s y(0) + laplace(y(t), t, s) =

        exp(-s)


         2
        s  laplace(y(t), t, s) + laplace(y(t), t, s) = exp(-s)


                               exp(-s)
                               -------
                                2
                               s  + 1


              sol := y(t) = Heaviside(t - 1) sin(t - 1)

> plot(rhs(sol),t=0..10);

Note how the Laplace transform integral of Dirac(t-1) strobes exp(-st) at t = 1, to get exp(-s*1) = exp(-s), which is the Laplace transform of Dirac(t-1), as shown above.

This is the end of Chapter 7. We will then begin Chapter 9, systems of simultaneous first order linear differential equations, where we will again have many reasons to work with Maple's Linear Algebra package.

Chapter 9 will be the last topic in the course.

Day 21:

It is easy to do integrals that involve the Dirac delta function, since no 
antiderivative is required. 

p. 412 #1.
> Int((t^2-1)*Dirac(t),t=-infinity..infinity);

                      infinity
                     /
                    |            2
                    |          (t  - 1) Dirac(t) dt
                    |
                   /
                     -infinity

> value(%);

                                  -1

Since the delta function is zero except when t = 0, the only time that  
t^2 - 1 is illuminated  is at the point  t = 1, when its value is  - 1. 
That is the value of the integral.

To compute the Laplace transform of Dirac(t-2), we use the definition, 
which is an integral. Hence, as noted above, our task is easy:

> Int(exp(-s*t)*Dirac(t-2),t=0..infinity);
> value(%);

                   infinity
                  /
                 |
                 |          exp(-s t) Dirac(t - 2) dt
                 |
                /
                  0


                              exp(-2 s)

The delta function is zero except at  t = 2,  so the only time that 
exp(-s*t) is illuminated is at  t - 2, when its value is  exp(-s*t). That 
is the value of the integral, and the required Laplace transform.

p. 406 #15.

We can solve certain integral equations which contain convolution 
integrals in much the same way that we solve constant coefficient 
differential equations by means of Laplace transforms.

> with(inttrans):
> integral_equation:=y(t)+3*Int(y(v)*sin(t-v),v=0..t)=t;

                                       t
                                      /
                                     |
      integral_equation := y(t) + 3  |   y(v) sin(t - v) dv = t
                                     |
                                    /
                                      0

Note how the Convolution Theorem is used to compute the Laplace transform 
of the convolution integral as a product of laplace(y(t),t,s) and 
1/(s^2+1).

Laplace the integral equation.
> laplace(integral_equation,t,s);

                                3 laplace(y(t), t, s)    1
          laplace(y(t), t, s) + --------------------- = ----
                                        2                 2
                                       s  + 1            s

There are no initial value to substitute, so we go directly to solving for 
the Laplace transform of the answer. 

(By hand we would factor out the Laplace transform, then divide by 
(1 + 3/(s^2+2) = (s^2 + 4)/(s^2 + 1). )

> solve(%,laplace(y(t),t,s));

                                2
                               s  + 1
                             -----------
                              2   2
                             s  (s  + 4)

Use partial fractions to convert to simple fractions.
> convert(%,parfrac,s);

                           1         3
                          ---- + ----------
                             2       2
                          4 s    4 (s  + 4)

Then take the inverse Laplace by the memorized tables, to get

> y(t)=(1/4*t+3/8*sin(2*t))*Heaviside(t);

               y(t) = (t/4 + 3/8 sin(2 t)) Heaviside(t)


p. 398 # 61 works exactly like the mixing problems in Chapter 3. The 
unknown is always q(t). The differential equation is
 
q' = (rate in)(conc in) - (rate out)(conc out). 

The volume is constant at 500 L, since  the rate in is the same as the rate out.

> diff(q(t),t)=piecewise(t<10,12*2/5-12*q(t)/500,t>10,12*3/5-12*q(t)/500);

             d         { 24/5 - 3/125 q(t)        t < 10
             -- q(t) = {
             dt        { 36/5 - 3/125 q(t)        10 < t

We can convert the piecewise function into Heaviside by the (new - old) 
Heaviside(t - transition) scheme, if one is asked to do it by hand, and 
you should check it.
 
Alternatively, by Maple,

> de:=convert(%,Heaviside);
> inits:=q(0)=500*1/5;

           d
     de := -- q(t) = 24/5 + 12/5 Heaviside(-10 + t) - 3/125 q(t)
           dt


                         inits := q(0) = 100

By the four step Laplace procedure,
> with(inttrans):
> laplace(de,t,s);
> subs(inits,%);
> solve(%,laplace(q(t),t,s));
> sol:=q(t)=invlaplace(%,s,t);
> 

  s laplace(q(t), t, s) - q(0) =

             2 + exp(-10 s)
        12/5 -------------- - 3/125 laplace(q(t), t, s)
                   s


  s laplace(q(t), t, s) - 100 =

             2 + exp(-10 s)
        12/5 -------------- - 3/125 laplace(q(t), t, s)
                   s


                    100 (125 s + 6 + 3 exp(-10 s))
                    ------------------------------
                            s (125 s + 3)


                           3 t
  sol := q(t) = -100 exp(- ---) + 200
                           125

                                  /          3 t        \
         + 100 Heaviside(-10 + t) |1 - exp(- --- + 6/25)|
                                  \          125        /

Since the problem asked not for the mass, but for the concentration, we 
need to divide by the volume.

> rhs(sol)/500;

             3 t
  -1/5 exp(- ---) + 2/5
             125

                                  /          3 t        \
         + 1/5 Heaviside(-10 + t) |1 - exp(- --- + 6/25)|
                                  \          125        /

p. 419 #31.

We can solve systems of simultaneous constant coefficient initial value 
problems by the Laplace transform. Laplace transforms each
equation into a "multiply by s" equation, and we then solve the 
simultaneous 8th grade algebra equations simultaneously, finally we use 
the inverse Laplace transform.

> de1:=diff(x(t),t)+y(t)=0;
> de2:=x(t)+diff(y(t),t)=1-Heaviside(t-2);
> inits:=x(0)=0,y(0)=0;

                            /d      \
                     de1 := |-- x(t)| + y(t) = 0
                            \dt     /


                          /d      \
            de2 := x(t) + |-- y(t)| = 1 - Heaviside(t - 2)
                          \dt     /


                     inits := x(0) = 0, y(0) = 0

> laplace({de1,de2},t,s);

                                                        1 - exp(-2 s)
  {laplace(x(t), t, s) + s laplace(y(t), t, s) - y(0) = -------------
                                                              s

        , s laplace(x(t), t, s) - x(0) + laplace(y(t), t, s) = 0}

> subs(inits,%);

  {s laplace(x(t), t, s) + laplace(y(t), t, s) = 0,

                                                      1 - exp(-2 s)
        laplace(x(t), t, s) + s laplace(y(t), t, s) = -------------}
                                                            s

> solve(%,{laplace(x(t),t,s),laplace(y(t),t,s)});

                         -1 + exp(-2 s)
  {laplace(x(t), t, s) = --------------,
                                   2
                          s (-1 + s )

                                -1 + exp(-2 s)
        laplace(y(t), t, s) = - --------------}
                                         2
                                   -1 + s

Maple writes the solution in terms of hyperbolic trig funtions, but we can 
convert them to exponentials.

> invlaplace(%,s,t);
> convert(%,exp):
> sol:=simplify(%);

  {y(t) = sinh(t) - Heaviside(t - 2) sinh(t - 2),

                                                             2
        x(t) = 1 - cosh(t) + 2 Heaviside(t - 2) sinh(t/2 - 1) }


  sol := {y(t) = 1/2 exp(t) - 1/2 exp(-t)

         - 1/2 Heaviside(t - 2) exp(t - 2)

         + 1/2 Heaviside(t - 2) exp(-t + 2), x(t) = 1 - 1/2 exp(t)

         - 1/2 exp(-t) + 1/2 Heaviside(t - 2) exp(t - 2)

         - Heaviside(t - 2) + 1/2 Heaviside(t - 2) exp(-t + 2)}

Checking that both sides are equal in each of the two differential 
equations,

> subs(sol,{de1,de2}):
> simplify(%);

         {0 = 0, 1 - Heaviside(t - 2) = 1 - Heaviside(t - 2)}



p. 395, #4. Use of the Theorem 8 to find Laplace transforms by hand.

Suppose that we wanted to find the Laplace transform  BY HAND  of   t 
Heaviside(t - 1). 

We want our children to look like us, so we rewrite the expression to 
match the (t - 1) inside Heaviside. 

Thus, our expression becomes   

(t - 1) Heaviside(t - 1)  + 1 Heaviside(t - 1).

The first term is  t Heaviside(t), but with  t  replaced by  t - 1 due to 
a multiplication by  exp(-1 s).  Hence, using Theorem 8 on p. 387, we get 
that the Laplace transform of  (t - 1) Heaviside(t - 1) is exp(- s) 1/s^2.  
(The memorized table tells us that the Laplace transform of t is 1/s^2.)    

The second term is  1(t - 1) Heaviside(t - 1), where the 1(t - 1) is not a 
multiplication,  It is the function which is one at each value of t, 
evaluated at (t - 1).  The replacement of t by  (t - 1) is due to a 
multiplication by exp(-1 s).  Hence,  using Theorem 8 on p. 387, we get 
that the Laplace transform of  1 Heaviside(t - 1) = 1(t - 1) Heaviside(t - 
1)  is exp(- s) 1/s.  (The memorized table tells us that the Laplace 
transform of 1 is 1/s.)

We can check with Maple,
 
> exp(-s)*1/s^2+exp(-s)*1/s;
> simplify(%);

                          exp(-s)   exp(-s)
                          ------- + -------
                             2         s
                            s


                           exp(-s) (s + 1)
                           ---------------
                                  2
                                 s

Check:
> t*Heaviside(t-1);
> laplace(%,t,s);

                          t Heaviside(t - 1)


                           exp(-s) (s + 1)
                           ---------------
                                  2
                                 s

An example of a system of simultaneous first order linear equations, solved two ways in Maple.

p. 537 Example 3.

One can use Maple to find formula solutions to systems. For example, we 
can find formulas for x1(t) and x2(t) in terms of t and two parameters C1 
and C2. (If we had initial values for x1(t) and x2(t), then we could solve 
for C1 and C2 to match those initial conditions.)

> restart:
> de1:=diff(x1(t),t)=2*x1(t)-3*x2(t);
> de2:=diff(x2(t),t)=x1(t)-2*x2(t);

                        d
                 de1 := -- x1(t) = 2 x1(t) - 3 x2(t)
                        dt


                         d
                  de2 := -- x2(t) = x1(t) - 2 x2(t)
                         dt

> sol:=dsolve({de1,de2},{x1(t),x2(t)});

  sol := {x1(t) = _C1 exp(t) + _C2 exp(-t),

        x2(t) = 1/3 _C1 exp(t) + _C2 exp(-t)}

We can check by subs and simplify.

> subs(sol,{de1,de2});
> simplify(%);

   d
  {-- (1/3 _C1 exp(t) + _C2 exp(-t)) = 1/3 _C1 exp(t) - _C2 exp(-t),
   dt

        d
        -- (_C1 exp(t) + _C2 exp(-t)) = _C1 exp(t) - _C2 exp(-t)}
        dt


  {1/3 _C1 exp(t) - _C2 exp(-t) = 1/3 _C1 exp(t) - _C2 exp(-t),

        _C1 exp(t) - _C2 exp(-t) = _C1 exp(t) - _C2 exp(-t)}

However, dsolve is seldom allowed, since it does not show me what you have 
learned. Hence, we will often just convert this homogeneous system to a 
single matrix equation, using the coefficient matrix.

> with(LinearAlgebra):
> A:=Matrix(<<2,1>|<-3,-2>>);

                                 [2    -3]
                            A := [       ]
                                 [1    -2]

Single matrix differential equation: X '  = A.X.

> X:=Matrix(<x1(t),x2(t)>);
> map(diff,X,t)=A.X;

                                  [x1(t)]
                             X := [     ]
                                  [x2(t)]


                   [d       ]
                   [-- x1(t)]
                   [dt      ]   [2 x1(t) - 3 x2(t)]
                   [        ] = [                 ]
                   [d       ]   [ x1(t) - 2 x2(t) ]
                   [-- x2(t)]
                   [dt      ]

These equations form a true system and cannot be solve separately.

Our task is to find a new coordinate system in which the system will 
decouple. If we represent vectors in a new basis of eigenvectors for
this system, rather than the default basis of { i , j }, then the matrix 
will diagonalize (all entries off the principle diagonal will become 
zero), and the problem will become easy. 

Although we will eventually learn how to find eigenvalues and eigenvectors 
by hand (at least for 2 by 2 systems), for the nonce we will use Maple's 
Eigenvectors command.

> v,C:=Eigenvectors(A);

                                [ 1]  [3    1]
                        v, C := [  ], [      ]
                                [-1]  [1    1]

The eigenvector  1 i + 1 j  belongs to eigenvalue  - 1, while the 
eigenvector  3 i + 1 j  belongs to the eigenvalue  1.  We "rotate 
coordinates" to use { v1 = i  + j , v2 = 3i + j } as our new coordinate 
axes directions.

This is equivalent to making the substitution  X = C.Y, which is the 
substitution that we will always use, based on the eigenvectors that form 
the columns of C.

> Y:=Matrix(<y1(t),y2(t)>);

                                  [y1(t)]
                             Y := [     ]
                                  [y2(t)]

Here is the substitution, written out.

> X = C.Y;

                     [x1(t)]   [3 y1(t) + y2(t)]
                     [     ] = [               ]
                     [x2(t)]   [ y1(t) + y2(t) ]

As noted in class, if  X = C.Y, then  X ' = C.Y '.

Hence the equation that was X ' = A.X, with dependent variables x1(t) and 
x2(t), becomes a new system with dependent variables y1(t) and y2(t) in 
the new coordinate system. 

The new system is  C.Y ' = A.C.Y, which is, again, a first order linear 
matrix differential equation, but it is not in standard form.

Hence, we divide by C on the left side of each side of the equation to get

Y ' = 1/C.A.C.Y.

> new_de:=map(diff,Y,t)=1/C.A.C.Y;

                             [d       ]
                             [-- y1(t)]
                             [dt      ]   [y1(t) ]
                   new_de := [        ] = [      ]
                             [d       ]   [-y2(t)]
                             [-- y2(t)]
                             [dt      ]

Note that each scalar equation now has only one dependent variable, so the 
system is decoupled. We can now solve each equation one at a time, without 
paying any attention to the other equations. In other words, we are back 
to SINGLE differential equations.

The reason that the system is decoupled is that the matrix  1/C.A.C has 
been DIAGONALIZED, meaning that all elements off the diagonal are zero.
(You may also note that the eigenvalues are *on* the diagonal, in the same
order that Maple placed the eigenvectors as columns of C.)

> 1/C.A.C;

                              [1     0]
                              [       ]
                              [0    -1]

Since we can solve new_de by inspection, we now have solutions for y1(t) 
and y2(t), which we place inside Y.

> subs({y1(t)=C1*exp(t),y2(t)=C2*exp(-t)},Y);

                             [C1 exp(t) ]
                             [          ]
                             [C2 exp(-t)]

Since our substitution is  X = C.Y, we can now find the answers to the 
original system, where x1(t) is on top and x2(t) is on the bottom. (Feel 
free to compare with dsolve, above.)

X = C.Y

> X:=C.%;

                        [3 C1 exp(t) + C2 exp(-t)]
                   X := [                        ]
                        [ C1 exp(t) + C2 exp(-t) ]

We can check our results in matrix form.

Since the original equation was  X ' = A.X, making  X ' - A.X = 0, we 
confirm:

> map(diff,X,t)-A.X;

                                 [0]
                                 [ ]
                                 [0]

The assignment for today is to master the above example. You can check your mastery on p. 541 #1, 2, 3, 4. (Note that I do not mean just to find the eigenvalues and eigenvectors as per the instructions in the book: use *the diagonalization procedure* to completely solve the system.)

This material IS fair game for April 18th. Do not procrastinate in learning it.

Day 22:

On Day 22 we completely solved p. 541 #3 using the diagonalization procedure, remarked on the shortcut that is ONLY available on homogeneous equations, then reworked the same example using diagonalization with an addition to the same problem that made it nonhomogeneous.

Note that we are not going to use undetermined coefficients in this 
problem. Undetermined coefficients will NOT be on the exam. 
Diagonalization will. 

(Undetermined coefficients often does not work very well in systems.)

Here is the data in the problem.

> restart:
> with(LinearAlgebra):
> A:=Matrix(<<6,4>|<1,3>>);
> f:=Matrix(<<-11,-5>>);

                                 [6    1]
                            A := [      ]
                                 [4    3]


                                   [-11]
                              f := [   ]
                                   [ -5]

The data represents the following nonhomogeneous system:  X ' = A X + f

> X:=Matrix(<x1(t),x2(t)>);
> map(diff,X,t)=A.X+f;

                                  [x1(t)]
                             X := [     ]
                                  [x2(t)]


                 [d       ]
                 [-- x1(t)]
                 [dt      ]   [6 x1(t) + x2(t) - 11 ]
                 [        ] = [                     ]
                 [d       ]   [4 x1(t) + 3 x2(t) - 5]
                 [-- x2(t)]
                 [dt      ]

We use Maple to find the eigenvalues and eigenvectors that match the 
coordinate system that the problem prefers.

> v,C:=Eigenvectors(A);

                               [7]  [1    -1/4]
                       v, C := [ ], [         ]
                               [2]  [1     1  ]

The substitution that we will make is  X = C Y.

> Y:=Matrix(<y1(t),y2(t)>);
> X=C.Y;

                                  [y1(t)]
                             Y := [     ]
                                  [y2(t)]


                    [x1(t)]   [y1(t) - 1/4 y2(t)]
                    [     ] = [                 ]
                    [x2(t)]   [  y1(t) + y2(t)  ]


After x1(t) has been replaced by x1(t) = y1(t) - 1/4 y2(t) and 
x2(t) has been replaced by x2(t) = y1(t) + y2(t),
all the x's will be gone and we will have a system of first order linear 
equations in y1(t) and y2(t).

The original problem is  X ' = A X + f.  The substitution is  X = C Y.  
After the substitution, the new system is  C Y ' = A C Y + f.  

When we put this system in standard form by multiplying on the left by 
1/C, we get a changed equation, which is much nicer than the original one. 
Since we have adapted our coordinate system away from the standard 
coordinate system to the one preferred by the problem, the new system will 
decouple and can be solved one equation at a time, rather than having to 
work with two equations at a time.  

Y ' = 1/C A C Y + 1/C f.

Note that because the matrix 1/C.A.C is diagonal - meaning that all the 
elements off the principal diagonal are zero - we only get ONE DEPENDENT 
VARIABLE  PER EQUATION: the system is DECOUPLED.

> map(diff,Y,t)=1/C.A.C.Y+1/C.f;

                    [d       ]
                    [-- y1(t)]
                    [dt      ]   [7 y1(t) - 49/5]
                    [        ] = [              ]
                    [d       ]   [2 y2(t) + 24/5]
                    [-- y2(t)]
                    [dt      ]

While the resulting system is decoupled, so it can be solved one equation 
at a time, in the nonhomogeneous case as contrasted with the homogeneous 
one, we cannot get the answer "by inspection." We have to do some honest 
labor, as in solving a first order linear equation in Math 152:

We can solve the first equation by hand, by putting it in standard form  
y1 '(t) - 7 y1(t) = -49/5, noting that p = -7, so mu = exp( - 7 t), etc., 
then using formula (8) on p. 50.

We can solve the second equation by hand by putting it in standard form
y2 '(t) - 2 y2(t) = 24/5, noting that p = -2, so mu = exp( - 2 t), etc, 
then using formula (8) on p. 50.

> p:=-7;
> int(%,t);
> mu:=exp(%);
> 1/mu*(int(mu*(-49/5),t)+C1):
> y1(t)=simplify(%);
> sol1:=expand(%):
> sol1:=simplify(%);

                               p := -7


                                 -7 t


                           mu := exp(-7 t)


              y1(t) = 1/5 exp(7 t) (7 exp(-7 t) + 5 C1)


                  sol1 := y1(t) = 7/5 + exp(7 t) C1

> p:=-2;
> int(%,t);
> mu:=exp(%);
> 1/mu*(int(mu*(24/5),t)+C2):
> y2(t)=simplify(%);
> sol2:=expand(%):
> sol2:=simplify(%);

                               p := -2


                                 -2 t


                           mu := exp(-2 t)


             y2(t) = -1/5 exp(2 t) (12 exp(-2 t) - 5 C2)


                 sol2 := y2(t) = - 12/5 + exp(2 t) C2

Now that we have the solutions to the easier, decoupled problem, we then 
go backwards through the substitution to find the solutions x1(t) and 
x2(t).

> Y=subs({sol1,sol2},Y);

                   [y1(t)]   [ 7/5 + exp(7 t) C1  ]
                   [     ] = [                    ]
                   [y2(t)]   [- 12/5 + exp(2 t) C2]

The substitution is  X = C Y, so the answers to the original problem are

> original_sols:=X=C.rhs(%);

                     [x1(t)]   [2 + exp(7 t) C1 - 1/4 exp(2 t) C2]
    original_sols := [     ] = [                                 ]
                     [x2(t)]   [ -1 + exp(7 t) C1 + exp(2 t) C2  ]

 
Since we have gone through numerous computations, we would like to check 
our answers. 

If  X ' = A X + f, then it follows that  X' - A X - f = 0, where 0 is not 
a number, but a matrix of numbers. 

Check:
> rhs(original_sols);
> map(diff,%,t)-A.%-f;

                 [2 + exp(7 t) C1 - 1/4 exp(2 t) C2]
                 [                                 ]
                 [ -1 + exp(7 t) C1 + exp(2 t) C2  ]


                                 [0]
                                 [ ]
                                 [0]


Note that since the system is linear and nonhomogeneous, the general 
solution can be written as a sum of a particular solution, plus all the 
homogeneous solutions.  (The homogeneous solutions can be found directly 
from the eigenvalue and an eigenvector belonging to it, as in the easy 
method described in class, but we need some other method to get a 
particular solution.)

> psol:=Matrix(<2,-1>);
> hsol1:=exp(7*t)*Matrix(<1,1>);
> hsol2:=exp(2*t)*Matrix(<-1/4,1>);

                                     [ 2]
                             psol := [  ]
                                     [-1]


                                  [exp(7 t)]
                         hsol1 := [        ]
                                  [exp(7 t)]


                                [-1/4 exp(2 t)]
                       hsol2 := [             ]
                                [  exp(2 t)   ]

> X=psol+C1*hsol1+C2*hsol2;

            [x1(t)]   [2 + exp(7 t) C1 - 1/4 exp(2 t) C2]
            [     ] = [                                 ]
            [x2(t)]   [ -1 + exp(7 t) C1 + exp(2 t) C2  ]

Tuesday's HW Assignment

In addition to using the diagonalization method to find complete solutions to homogeneous systems p. 541 # 1, 2, 3 and 4, you should now be equipped to use diagonalization to find complete solutions to nonhomogeneous systems.

p. 556 # 8, 13, 21, 35.

Again, ignore directions in the book as to how to do these problems; USE DIAGONALIZATION.

Day 23:

Conversion of Higher Order Linear Equations to First Order Linear Systems

For example, p. 336 #1. The third order nonhomogeneous equation
y''' -2 y'' -5 y' + 6y = exp(x) + x^2 can be converted to a system of three first order linear equations.

Original equation, solved for y''':

y''' = -6 y + 5 y' + 2 y'' + exp(x) + x^2

The Rosetta Stone:

z1 = y
z2 = y'
z3 = y''

The nonhomogeneous system:

z1 ' =           z2
z2 ' =                  z3
z3 ' = -6 z1 + 5 z2 + 2 z3 + exp(x) + x^2

The trade is a fair one, since, if we can solve the scalar third order 
equation, then we know y. But if we know y, then we know the answer z1, 
z2, and z3, which is the answer to the system. Conversely, if we know the
answer to the system, then we certainly know z1; and if we know z1 then we
know y, the answer to the scalar third order equation.

Fundamental Solution Matrix

Just as in Chapters 4 and 6 we had a fundamental solution set from which 
all homogeneous solutions could be made, so too in Chapter 9 we have a 
fundamental solution matrix.


We use the shortcut method to find a fundamental solution matrix for a
homogeneous system.



> with(LinearAlgebra):
> A:=Matrix(<<0,1,1>|<1,0,1>|<1,1,0>>);

                               [0    1    1]
                               [           ]
                          A := [1    0    1]
                               [           ]
                               [1    1    0]

> v,C:=Eigenvectors(A);

                            [-1]  [-1    -1    1]
                            [  ]  [             ]
                    v, C := [-1], [ 0     1    1]
                            [  ]  [             ]
                            [ 2]  [ 1     0    1]

We make three homogeneous solutions using the eigenvalues and an 
eigenvector for each of them. (This is the quickie method described in 
class and in 9.5.)

> hsol1:=exp(-1*t)*Matrix(<-1,0,1>);
> hsol2:=exp(-1*t)*Matrix(<-1,1,0>);
> hsol3:=exp(2*t)*Matrix(<1,1,1>);

                                  [-exp(-t)]
                                  [        ]
                         hsol1 := [   0    ]
                                  [        ]
                                  [exp(-t) ]


                                  [-exp(-t)]
                                  [        ]
                         hsol2 := [exp(-t) ]
                                  [        ]
                                  [   0    ]


                                  [exp(2 t)]
                                  [        ]
                         hsol3 := [exp(2 t)]
                                  [        ]
                                  [exp(2 t)]


We make a fundamental solution matrix by laying in the three solutions as 
columns of a matrix F.

> F:=Matrix(<hsol1|hsol2|hsol3>);

                    [-exp(-t)    -exp(-t)    exp(2 t)]
                    [                                ]
               F := [   0        exp(-t)     exp(2 t)]
                    [                                ]
                    [exp(-t)        0        exp(2 t)]


All the homogeneous solutions can be made by taking linear combinations of 
the columns of F.

> C:=Matrix(<C1,C2,C3>);

                                   [C1]
                                   [  ]
                              C := [C2]
                                   [  ]
                                   [C3]


We can make a linear combination  C1 hsol1  +  C2 hsol2 + C3 hsol3  by 
matrix multiplication:

> F.C;

               [-exp(-t) C1 - exp(-t) C2 + exp(2 t) C3]
               [                                      ]
               [       exp(-t) C2 + exp(2 t) C3       ]
               [                                      ]
               [       exp(-t) C1 + exp(2 t) C3       ]


To check that the columns of F are linearly independent, we compute the 
Wronskian determinant.

With scalar equations, we had to take derivatives to make a square matrix; 
but here no derivatives are required, since F is already a square matrix.

Since the determinant is not zero, the solutions are linearly independent.

> Determinant(F);
> simplify(%);

                                   2
                         -3 exp(-t)  exp(2 t)


                                  -3
To check that the columns of F are solutions to the equation X ' = A X, or 
equivalently, X ' - A X = 0, we compute

> map(diff,F,t)-A.F;
                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

Since the first column is all zeros, then the first column of F is a 
solution. Ditto for the othet two solutions. 



 Eigenvectors

Definition: A nonzero vector v is an EIGENVECTOR for a square matrix A if 
A.v is a scalar multiple of itself. If A.v = r*v, then r is said to be 
an EIGENVALUE for A and v is an eigenvector belonging to r.


Example:


> restart:
> with(LinearAlgebra):
> A:=Matrix(<<0,1,1>|<1,0,1>|<1,1,0>>);

                               [0    1    1]
                               [           ]
                          A := [1    0    1]
                               [           ]
                               [1    1    0]

> v,C:=Eigenvectors(A);

                            [-1]  [-1    -1    1]
                            [  ]  [             ]
                    v, C := [-1], [ 0     1    1]
                            [  ]  [             ]
                            [ 2]  [ 1     0    1]

The first column of C is an eigenvector for A belonging to the eigenvalue  
r = - 1.
> A.Column(C,1);
> -1*Column(C,1);

                                 [ 1]
                                 [  ]
                                 [ 0]
                                 [  ]
                                 [-1]


                                 [ 1]
                                 [  ]
                                 [ 0]
                                 [  ]
                                 [-1]

The second column of C is an eigenvector for A belonging to the eigenvalue  
r = -1.
> A.Column(C,2);
> -1*Column(C,2);

                                 [ 1]
                                 [  ]
                                 [-1]
                                 [  ]
                                 [ 0]


                                 [ 1]
                                 [  ]
                                 [-1]
                                 [  ]
                                 [ 0]

The third column of C is an eigenvector for A belonging to the eigenvalue 
r = 2.
> A.Column(C,3);
> 2*Column(C,3);


                                 [2]
                                 [ ]
                                 [2]
                                 [ ]
                                 [2]


                                 [2]
                                 [ ]
                                 [2]
                                 [ ]
                                 [2]


On the other hand, the vector <1,2,3> is NOT an eigenvector for A since  
<5, 4, 3>  is not a scalar multiple of  <1, 2, 3>.

> v:=Matrix(<1,2,3>);
> A.v;

                                    [1]
                                    [ ]
                               v := [2]
                                    [ ]
                                    [3]


                                 [5]
                                 [ ]
                                 [4]
                                 [ ]
                                 [3]



 Since we now know what an eigenvector is, we can go looking for one by 
hand in the case where A is a 2 by 2 matrix.

 Here is an example:


If  v is a nonzero vector which is an eigenvector for A, then  A v = r v  
for some scalar r.

That definition says that  A v = r Id v,  where Id is the multiplicative 
identity matrix with ones down the diagonal and zeros elsewhere. (  r Id v 
= r v, since the identity matrix behaves like the number 1.)

However, if A v = r Id v,  then  A v - r Id v = 0. (That's a matrix of 
zeros, not the number zero.)

Factoring out the  v, we get that  (A - r Id) v = 0, which is what we call 
B below.  

B is just the A matrix with r's subtracted down the diagonal.

> restart:
> with(LinearAlgebra):
> A:=Matrix(<<6,4>|<1,3>>);

                                 [6    1]
                            A := [      ]
                                 [4    3]

> Id:=Matrix(<<1,0>|<0,1>>);

                                  [1    0]
                            Id := [      ]
                                  [0    1]

> B:=A-r*Id;

                             [6 - r      1  ]
                        B := [              ]
                             [  4      3 - r]

Now if B v = 0, then if 1/B exists, we would be forced to accept  
v = 1/B 0 = 0 as our only eigenvector. However, the definition of eigenvector 
given above explicitly excludes v = 0 as a choice for eigenvector.

That leaves us no choice. We must force B NOT to have a multiplicative 
inverse, and that means that we must force B to have zero determinant.  

(Recall that to find the inverse of a two by two matrix B, we switched 
places on the diagonal, switched signs on the off diagonal, and divided by the 
determinant of B. Hence, if we have a zero determinant for B, there is no 
inverse for B.)

Setting the determinant of B to zero is how we get the characteristic 
equation which we use to solve for eigenvalues  r.
 
Here is how we solve for the eigenvalues r = 7 and r = 2.

> char_eqn:=Determinant(B)=0;

                                            2
                    char_eqn := 14 - 9 r + r  = 0

> solve(%,{r});

                           {r = 7}, {r = 2}


Now that we have the eigenvalues r = 7 and r = 2, we substitute each into 
the matrix B and try to solve the system  (A -r Id) v = 0.

Here is the first eigebvector, using r = 7. 

If you look carefully, you will note that what at first blush appeared to 
be a system of two equations in two unknowns is really a system of one 
equation in two unknowns. (The second equation is a multiple of the first 
equation.) This effect is NOT an error, you did it yourself when 
you forced the determinant to be zero.)

Since we only have one equation, we are not restricted to using only  
<0, 0> as an eigenvector.

If we pick v1 = 1, then we can solve to get  v2 = 1. Therefore <1, 1> is 
an eigenvector belonging to r = 7. (If you pick any other value for v1 
besides 1, that choice will also give an eigenvector.)

> subs(r=7,B).Matrix(<v1,v2>)=Matrix(<0,0>);;

                         [ -v1 + v2  ]   [0]
                         [           ] = [ ]
                         [4 v1 - 4 v2]   [0]


Since we only have one equation, we are not restricted to using only  
< 0, 0 > as an eigenvector.

If we pick v1 = - 1/4, then we can solve to get v2 = 1. Therefore 
< - 1/4, 1 > is an eigenvector belonging to r = 7.

> subs(r=2,B).Matrix(<v1,v2>)=Matrix(<0,0>);

                          [4 v1 + v2]   [0]
                          [         ] = [ ]
                          [4 v1 + v2]   [0]

Checking:

> v,C:=Eigenvectors(A);

                               [2]  [-1/4    1]
                       v, C := [ ], [         ]
                               [7]  [ 1      1]


Any scalar multiple of an eigenvector belonging to A is also an 
eigenvector belonging to A.

In particular, had you chosen  v1 = 1 in the section immediately above, 
you would have gotten v2 = -4, producing an eigenvector  < 1, -4 > 
belonging to r = 2. That is also correct. (Note that < 1, -4 > is a scalar 
multiple of Maple's result.)

Again, we verify from the definition that the vectors indicated are 
eigenvectors for A.

> A.Matrix(<-1/4,1>);
> 2*Matrix(<-1/4,1>);
> 

                                [-1/2]
                                [    ]
                                [ 2  ]


                                [-1/2]
                                [    ]
                                [ 2  ]

> A.Matrix(<1,1>);
> 7*Matrix(<1,1>);

                                 [7]
                                 [ ]
                                 [7]


                                 [7]
                                 [ ]
                                 [7]

Note: ANY choice of eigenvectors will work in the matrix C that we use for 
diagonalization. There is nothing special about Maple's choice.



 Thursday's HW Assignment:


 You are now ready to do

 p. 508 # 11, p. 521 #6a, 13 using Maple, 15, 19 
using Maple, 25 using Maple, 29, 34 both by hand and by Maple's 
map(diff,%,t), 39, and
 
 p. 530 #1, 12, 18, 22,  23.
 
 Plus, you should go back and find eigenvalues and eigenvectors of 
p. 541 #1, 2, 3, 4 by hand.

Day 24:

On Day 24 we began to study the final two section in the course: complex 
eigenvalues and eigenvectors, and systems that will not diagonalize.


Here is an example of a Jordan Block, which will not diagonalize.



Here is a Jordan Block. Note that the eigenvalues run along the diagonal, 
with 1's in the first superdiagonal. That is already in the form that is 
closest to diagonal that we can hope to achieve.  (With a diagonal matrix, 
ALL elements off the diagonal are zero; but here we have 1's along the 
first superdiagonal. Every other off diagonal entry is zero.)

We can try to diagonalize, but there are not enough linearly independent 
eigenvectors to make a basis.

> restart:
> with(LinearAlgebra):
> A:=Matrix(<<2,0,0>|<1,2,0>|<0,1,2>>);

                               [2    1    0]
                               [           ]
                          A := [0    2    1]
                               [           ]
                               [0    0    2]

> v,C:=Eigenvectors(A);

                              [2]  [1    0    0]
                              [ ]  [           ]
                      v, C := [2], [0    0    0]
                              [ ]  [           ]
                              [2]  [0    0    0]


If you wanted to see what goes on by hand, note that the eigenvalues are 
2, 2, and 2. When we substract 2's down the diagonal, the equations that 
we need to solve for v1, v2, abd v3 are

> Matrix(<<0,0,0>|<1,0,0>|<0,1,0>>).Matrix(<<v1,v2,v3>>)
     = Matrix(<0,0,0>);

                              [v2]   [0]
                              [  ]   [ ]
                              [v3] = [0]
                              [  ]   [ ]
                              [0 ]   [0]


v2 and v3 must be zero, but v1 can be any value except 0 . (Recall that < 
0, 0, 0 > is NOT counted as an eigenvector.) Choosing v1 = 1, we get 
Maple's answer:  the eigenspace is one dimensional and consists of all 
scalar multiples of < 1, 0, 0 >. In particular, we do not get three 
linearly
independent columns for C that would be necessary for C to be invertible.

The matrix A CANNOT BE DIAGONALIZED. Hence, the system cannot be 
decoupled.

(An astute observer might note that in the case of a Jordan block, the 
system is "close enough" to diagonal in order for us to solve the system.
Note that the last equation can be solved without considering the other 
differential equations, since it contains only one dependent variable.
Then we can substitute that answer in the next to the last equation, which 
will then contain only one dependent variable, etc.  

However, there are many matrices which are NOT Jordan Blocks, but which 
also cannot be diagonalized.)




Here is an example of the use of the matrix exponential to solve a 
nonhomogeneous system whose coefficient matrix A does not diagonalize.

p. 567 #23.

Here is the data for the problem.

> restart:
> with(LinearAlgebra):
> 
X:=Matrix(<x1(t),x2(t),x3(t)>);
A:=Matrix(<<2,-3,9>|<1,-1,3>|<-1,1,-4>>);
f:=Matrix(<0,t,0>);

                                  [x1(t)]
                                  [     ]
                             X := [x2(t)]
                                  [     ]
                                  [x3(t)]


                             [ 2     1    -1]
                             [              ]
                        A := [-3    -1     1]
                             [              ]
                             [ 9     3    -4]


                                    [0]
                                    [ ]
                               f := [t]
                                    [ ]
                                    [0]

Here is the system with the letters put in. 

> map(diff,X,t)=A.X+f;

             [d       ]
             [-- x1(t)]
             [dt      ]
             [        ]   [  2 x1(t) + x2(t) - x3(t)   ]
             [d       ]   [                            ]
             [-- x2(t)] = [-3 x1(t) - x2(t) + x3(t) + t]
             [dt      ]   [                            ]
             [        ]   [9 x1(t) + 3 x2(t) - 4 x3(t) ]
             [d       ]
             [-- x3(t)]
             [dt      ]

Note that A does not diagonalize.

> v,C:=Eigenvectors(A);

                            [-1]  [1/3    0    0]
                            [  ]  [             ]
                    v, C := [-1], [ 0     0    0]
                            [  ]  [             ]
                            [-1]  [ 1     0    0]

C is not invertible.

> Determinant(C);
> 1/C;

                                  0

Error, (in rtable/Power) singular matrix

Now we apply formula (8) on p. 50, where the integrating factor  
mu is given by the matrix exponential.  (The standard form equation is 
X ' -Ax = f, so we use  p = -A, not A.)

> p:=-A;
> map(int,%,t);
> mu:=MatrixExponential(%);

                             [-2    -1     1]
                             [              ]
                        p := [ 3     1    -1]
                             [              ]
                             [-9    -3     4]


                        [-2 t     -t      t ]
                        [                   ]
                        [3 t      t      -t ]
                        [                   ]
                        [-9 t    -3 t    4 t]


  mu :=

        [                           2
        [exp(t) - 3 t exp(t) - 3/2 t  exp(t) , -t exp(t) ,

             2                  ]
        1/2 t  exp(t) + t exp(t)]

        [3 t exp(t) , exp(t) , -t exp(t)]

        [                   2
        [-9 t exp(t) - 9/2 t  exp(t) , -3 t exp(t) ,

                      2                    ]
        exp(t) + 3/2 t  exp(t) + 3 t exp(t)]

The formula on p. 50 gives

> mu.f;
> map(int,%,t)+Matrix(<c1,c2,c3>);
> 1/mu.%:
> X:=map(simplify,%);

                            [   2        ]
                            [ -t  exp(t) ]
                            [            ]
                            [  t exp(t)  ]
                            [            ]
                            [    2       ]
                            [-3 t  exp(t)]


                   [              2              ]
                   [ -(2 - 2 t + t ) exp(t) + c1 ]
                   [                             ]
                   [    (-1 + t) exp(t) + c2     ]
                   [                             ]
                   [               2             ]
                   [-3 (2 - 2 t + t ) exp(t) + c3]


  X :=

        [                                   2
        [-2 + t + exp(-t) c1 - 3/2 exp(-t) t  c1 + 3 exp(-t) t c1

                                       2                  ]
         + exp(-t) t c2 + 1/2 exp(-t) t  c3 - exp(-t) t c3]

        [t - 3 exp(-t) t c1 - 1 + exp(-t) c2 + exp(-t) t c3]

        [              2
        [-9/2 exp(-t) t  c1 + 3 t + 9 exp(-t) t c1 + 3 exp(-t) t c2

                                                          2   ]
         - 6 + exp(-t) c3 - 3 exp(-t) t c3 + 3/2 exp(-t) t  c3]


Checking our answer, using matrix multiplication: if  X ' = A.X = f, then 
X ' - A.X - f = 0. 

> map(diff,X,t)-A.X-f:
> simplify(%);

                                 [0]
                                 [ ]
                                 [0]
                                 [ ]
                                 [0]

It is interesting to note that the matrix exponential gives us a 
fundamental solution matrix. (WARNING:  A, not - A. Note that mu was the 
matrix exponential of  -A*t, but we took 1/mu in the formula on p. 50, 
which is the matrix exponential of A*t.)

> F:=MatrixExponential(A*t);

  F :=

        [                             2
        [exp(-t) + 3 t exp(-t) - 3/2 t  exp(-t) , t exp(-t) ,

             2                    ]
        1/2 t  exp(-t) - t exp(-t)]

        [-3 t exp(-t) , exp(-t) , t exp(-t)]

        [                   2
        [9 t exp(-t) - 9/2 t  exp(-t) , 3 t exp(-t) ,

                                     2        ]
        exp(-t) - 3 t exp(-t) + 3/2 t  exp(-t)]

> map(diff,F,t)-A.F;

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

> Determinant(F):
> simplify(%);

                              exp(-3 t)





 You are invited to look up a more complete formula in the Lab Manual which
uses definite integrals, rather than adding c's, and which handles 
nonhomogeneous initial value problems with a single integral. p. 25?


Here is an example of a homogeneous system with complex eigenvalues.

p. 549 #3.

Since Maple assumes that ALL variables are complex, while you assume that 
t is real, we need to communicate to Maple what you are assuming, then we 
turn off the tildes that are applied to assumed variables.

> restart:
> with(LinearAlgebra):
> assume(t,real):
> interface(showassumed=0):
> A:=Matrix(<<1,0,0>|<2,1,-1>|<-1,1,1>>);

                              [1     2    -1]
                              [             ]
                         A := [0     1     1]
                              [             ]
                              [0    -1     1]


Note that since the data from the original problem contains real numbers, 
then if complex eivenvalues and eigenvectors occur, they must occur
as complex conjugate pairs.

> v,C:=Eigenvectors(A);

                       [1 + I]  [-2 + I    -2 - I    1]
                       [     ]  [                     ]
               v, C := [1 - I], [  -I        I       0]
                       [     ]  [                     ]
                       [  1  ]  [  1         1       0]

We get two solution vectors by taking the real and imaginary parts of 
*either* complex solution vector. The other solution vector is real, and 
that is laid in, too.

> exp((1+I)*t);
> evalc(%);
> %*Matrix(<-2+I,-I,1>);
> map(expand,%);
> hsol1:=map(Re,%);
> hsol2:=map(Im,%%);

                            exp((1 + I) t)


                   exp(t) cos(t) + exp(t) sin(t) I


                            [(-2 + I) %1]
                            [           ]
                            [   -I %1   ]
                            [           ]
                            [    %1     ]

                  %1 := exp(t) cos(t) + exp(t) sin(t) I


        [-2 exp(t) cos(t) + exp(t) cos(t) I - 2 I exp(t) sin(t)

         - exp(t) sin(t)]

        [-exp(t) cos(t) I + exp(t) sin(t)]

        [exp(t) cos(t) + exp(t) sin(t) I]


                      [-2 exp(t) cos(t) - exp(t) sin(t)]
                      [                                ]
             hsol1 := [         exp(t) sin(t)          ]
                      [                                ]
                      [         exp(t) cos(t)          ]


                       [exp(t) cos(t) - 2 exp(t) sin(t)]
                       [                               ]
              hsol2 := [        -exp(t) cos(t)         ]
                       [                               ]
                       [         exp(t) sin(t)         ]

> hsol3:=exp(1*t)*Matrix(<1,0,0>);

                                   [exp(t)]
                                   [      ]
                          hsol3 := [  0   ]
                                   [      ]
                                   [  0   ]

We can form a fundamental solution matrix by laying in the homogeneous 
solutions as columns of F.

> F:=Matrix(<<hsol1|hsol2|hsol3>>);

  F :=

        [-2 exp(t) cos(t) - exp(t) sin(t) ,

        exp(t) cos(t) - 2 exp(t) sin(t) , exp(t)]

        [exp(t) sin(t) , -exp(t) cos(t) , 0]

        [exp(t) cos(t) , exp(t) sin(t) , 0]

> map(diff,F,t)-A.F;

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

> Determinant(F):
> simplify(%);

                               exp(3 t)

We can get all solutions from a fundamental solution matrix.

> X=F.Matrix(<c1,c2,c3>);

  X =

        [(-2 exp(t) cos(t) - exp(t) sin(t)) c1

         + (exp(t) cos(t) - 2 exp(t) sin(t)) c2 + exp(t) c3]

        [exp(t) sin(t) c1 - exp(t) cos(t) c2]

        [exp(t) cos(t) c1 + exp(t) sin(t) c2]

It is perhaps worthwhile to remark that we can also use the matrix 
exponential to find still a different fundamental solution matrix. 

(Yes, you need to learn both ways.)

> F1:=MatrixExponential(A*t);

  F1 :=

        [exp(t) , 2 exp(t) sin(t) - exp(t) cos(t) + exp(t) ,

        -exp(t) sin(t) - 2 exp(t) cos(t) + 2 exp(t)]

        [0 , exp(t) cos(t) , exp(t) sin(t)]

        [0 , -exp(t) sin(t) , exp(t) cos(t)]

> map(diff,F1,t)-A.F1;

                            [0    0    0]
                            [           ]
                            [0    0    0]
                            [           ]
                            [0    0    0]

> Determinant(F1):
> simplify(%);

                               exp(3 t)



 

 Using Maple and the examples above, you should now be equipped to work 
on
 
p. 549 # 1,7,9,13,17,21

p. 555 # 8,13,21,35

p. 566 # 24.


Day 25:

 On Day 25 we examined a predator-prey system of simultaneous nonlinear 
differential equations that population biologists use.
 

A Predator-Prey example.

Note that the equations are nonlinear since we have a product of dependent 
variables. (The number of interactions per day is proportional to the 
product of the number of rabbits and the number of foxes.)

x(t) is the number of rabbits at time t years; y(t) is the number of foxes 
at time t years.  The loop shows that the sizes of the populations
follows an approximate eight year cycle.

> restart:
> de1:=diff(x(t),t)=x(t)*(1-1/2*y(t));
> de2:=diff(y(t),t)=y(t)*(-3/4+x(t)/4);
> inits:=x(0)=3,y(0)=0.5;

                        d
                 de1 := -- x(t) = x(t) (1 - 1/2 y(t))
                        dt


                      d
               de2 := -- y(t) = y(t) (- 3/4 + 1/4 x(t))
                      dt


                    inits := x(0) = 3, y(0) = 0.5

> 
sol:=dsolve({de1,de2,inits},{x(t),y(t)},numeric,method=classical[rk4],stepsize=0.01);

               sol := proc(x_classical)  ...  end proc

> with(plots):
> odeplot(sol,[x(t),y(t)],t=0..8);


 You should note that with systems, as with single equations, nonlinearity 
implies that formula solutions are seldom possible. We use a two 
dimensional version of the rk4 numeric approcimation. (We could also have 
used the forward Euler method. In the Euler method, we follow "arrows." 
In Math 251, what do we call the "arrow" [x'(t),y'(t)]?)

 The general instructions are that you do not need to learn the 
undetermined coefficient and "variation of parameter" methods.

 The undetermined coefficient method does not work well with systems, since
we need to multiply the initial forms of solutions by x^s, where s is 
chosen to be the smallest s that assures that no undetermined coefficients 
will die. Since the value for s may not be the same for both equations,
we cannot always satisfy both requirements. The HW problem that suggests
undetermined coefficients shows you the mess that we get.

 As for the "variation of parameter" method discussed on p. 553 (11), the
mysterious X(t) in formula (11) is just any fundamental solution matrix.
Since you are now past masters of the art of making and verifying 
fundamental matrices, you could easily use (11) to solve problems.

 However, before you decide to memorize yet another meaningless formula,
note that p. 553 (11) and p. 50 (8) are exactly the same object:


p.  50  (8)   x(t) = 1/mu * ( int( mu*rhs,     t) + C ) 
p. 553 (11)   x(t) = X(t) . ( int( 1/X(t).rhs, t) + C )



 since, as noted in class, mu = exp(-A*t), so 1/mu = exp(A*t) = X(t) is
a fundamental solution matrix.

 Similarly, you might want to compare p. 554 (13) with the Lab Manual's p.
26, #5.


 General Advice:


 Keep up to date. Make sure that you do two to three hours outside of class
for each hour in class. Study daily, not just once a week. It will stay 
with you longer. All quizzes and exams are cumulative.

 Be sure to bring questions for the start of class on the problems that you 
have attempted and the material that you have read. Your job is to make 
sure that any questions that you have get answered.
 

Those problems that involve extensive computation should be attempted in 
Maple, as well as by hand. Theory problems usually do not require Maple 
since the computations are usually minimal.

Download sample worksheets:
HW Solutions->right click week1.mw->Save Link Target
As->Save

Worksheet is now downloaded, open in Maple.


 What Maple Worksheets Should I Look at Before the Next Class?

10. c) Use of DEplot for direction fields.

5. Euler's Method via dsolve/numeric.

 Variables Separable, First Order Linear worksheet.

Click  here for a Bernoulli example in Maple.

 
A  good exam overview :  click here.


Here is the library for future reference:


10. c) Use of DEplot for direction fields.

NSS p.31
 
7. Phase line (stable and unstable equilibria).

8. Phase line (nodes).

NSS, p. 35

5. Euler's Method via dsolve/numeric.

Variables Separable, First  Order Linear  worksheet.

The Bernoulli worksheet immediately below in this page.

A  good exam overview :  click here, 
click here, click here.

Click  here for a Bernoulli example in Maple.

Click  here for a Reduction of Order example in Maple.

Click  here for a Variation of Parameters example in Maple.

Click  here for a Nonhomogenous System example in Maple.

Click  here for a System with Complex Eigenvalues done in Maple.

Click  here for suggested Laplace HW 
problems, solutions, and an introductory worksheet.

Click  here for Matrix Exponential.