Part III. First Order Ordinary Differential Equations

Section 1. Entering, Solving, Plotting

7.* Enter the differential equation y' + y = cos(t) with the name DE. Obtain the solution satisfying the generic initial condition y(t0) = y0. Name it soln.

> DE := diff(y(t),t) + y(t) = cos(t);
soln := dsolve( {DE, y(t0)=y0} );

DE := (diff(y(t), t))+y(t) = cos(t)

soln := y(t) = 1/2*cos(t)+1/2*sin(t)-exp(-t)*(-y0+1/2*cos(t0)+1/2*sin(t0))/exp(-t0)

a. Use the entry  eval( soln, t=t0)  to verify that soln satisfies the initial condition.

> eval(soln,t=t0);

y(t0) = y0

b. Use unapply to make rhs(soln) into a function of t, t0, and y0 named phi (see page 35 in Part III, Section 1 of the manual).

> phi := unapply(rhs(soln),t,t0,y0);

phi := proc (t, t0, y0) options operator, arrow; 1/2*cos(t)+1/2*sin(t)-exp(-t)*(-y0+1/2*cos(t0)+1/2*sin(t0))/exp(-t0) end proc

c. Use phi to plot some solutions starting at points evenly spaced on the y axis.

> plot( [ [t,phi(t,0,y0),t=0..10] $ y0=-2..2] );

[Plot]

d. Use phi to plot some solutions starting at points evenly spaced on the t axis.

> plot( [ [t,phi(t,t0,0),t=t0..10] $ t0=-2..2] );

[Plot]

e. Use phi to plot the solutions starting at points evenly spaced around the unit circle in the style of the two plots on page 35 of the manual. That is, one picture runs time forward, another runs time backward.

> Points := [cos(k*Pi/6),sin(k*Pi/6)] $ k=0..11:
plot( [ [t,phi(t,cos(k*Pi/6),sin(k*Pi/6)),t=cos(k*Pi/6)..2]

        $ k=0..11, [ Points ] ], style=[line$12,point],

        color=[red$12,black], view=[-2..2,-2..2] );

[Plot]

> plot( [ [t,phi(t,cos(k*Pi/6),sin(k*Pi/6)),t=-2..cos(k*Pi/6)]
        $ k=0..11, [ Points ] ], style=[line$12,point],

        color=[red$12,black], view=[-2..2,-2..2] );

[Plot]

Section 2. Working with Solutions: Modeling

3.* The following problem is adapted from Ledder, Chapter 1, Section 1.

Suppose you borrow $12,000 to buy a car. The loan is to be paid in 60 equal monthly installments at an interest rate of 5% per year.

a. Assume the payments are actually made continuously at whatever rate is needed to pay off the loan in 60 months. Determine the continuous rate per month that would be required.

> DE := diff(P(t),t) = 0.05*P(t) - Q;
soln := dsolve( {DE,P(0)=12000} );

DE := diff(P(t), t) = 0.5e-1*P(t)-Q

soln := P(t) = 20*Q+exp(1/20*t)*(12000-20*Q)

> eqn := subs(P(t)=0,t=5,soln);
Q := solve(eqn,Q);

eqn := 0 = 20*Q+exp(1/4)*(12000-20*Q)

Q := 600*exp(1/4)/(-1+exp(1/4))

> Monthly_Rate = evalf(Q/12)*dollars;

Monthly_Rate = 226.0405830*dollars

b. Compare the answer to part a to the answer if 60 equal monthly payments are made at a constant annual interest rate of 5% applied to the outstanding balance. In other words, the first payment, due one month after the loan is made, would be

12000*(exp(Float(5, -2)/12)-1)+P[1]

where P[1] is the amount that is put towards reducing the principal (the $12000) in the first month. The second payment is

(12000-P[1])*(exp(Float(5, -2)/12)-1)+P[2]

where P[2] is the amount put towards reducing the principal in the second month. The nth payment is

(12000-(sum(P[k], k = 1 .. n-1)))*(exp(Float(5, -2)/12)-1)+P[n] .

Hint for b. The second payment is supposed to equal to the first. Therefore, the following equation must be satisfied.

> eqn := 12000*(exp(0.05/12)-1)+P[1] =
      (12000 - P[1])*(exp(0.05/12)-1) + P[2];

eqn := 50.10431+P[1] = 50.10430800-0.4175359e-2*P[1]+P[2]

This determines P[2] in terms of P[1].

> P[2] := solve(eqn,P[2]);

P[2] := 0.2000000000e-5+1.004175359*P[1]

Create a for..do loop that calculates P[3], P[4], ... , P[60] in terms of P[1], then find P[1] using the fact that

sum(P[k], k = 1 .. 60) = 12000

> for n from 3 to 60
 do

   eqn := (12000-add(P[k],k=1..n-2))*(exp(0.05/12) - 1) + P[n-1] =

          (12000-add(P[k],k=1..n-1))*(exp(0.05/12) - 1) + P[n]:

   P[n] := solve(eqn,P[n]):

 end do:

unassign('n');

> add(P[k],k=1..60) = 12000:
P[1] := solve(%,P[1]);

P[1] := 176.4078443

> Check;
for n from 1 to 60 by 19

 do

  payment,n,equals,(12000-add(P[k],k=1..n-1))*(exp(0.05/12) - 1) +

  P[n];

 end do;

unassign('n');

Check

payment, 1, equals, 226.5121523

payment, 20, equals, 226.5121541

payment, 39, equals, 226.5121538

payment, 58, equals, 226.5121537

P[1] is the first payment towards the principal. The 60th payment is almost all towards the principal.

> Last_Principal_payment=P[60];

Last_Principal_payment = 225.5703165

Section 3. Slope Fields: DEplot

2.* Use DEplot to make a nice-looking slope field for the autonomous equation y' = sin(y). (Use the window t=-6..6, y=-6..6.)

> with(DEtools):

> DE := diff(y(t),t) = sin(y(t));
DEplot( DE, y(t), t=-6..6, y=-6..6, arrows=line, dirgrid=[17,17] );

DE := diff(y(t), t) = sin(y(t))

[Plot]

a. Put some solution curves into the plot by adding { [y(0)=k] $ k=-6..6 }, linecolor=blue . Comment on the relationship between one curve and the next.

> DEplot( DE, y(t), t=-6..6, y=-6..6, arrows=line, dirgrid=[17,17],
        { [y(0)=k] $ k=-6..6 }, linecolor=blue);

[Plot]

Each curve is a horizontal translate of the other.

b. Obtain the general solution using dsolve.

> dsolve( DE );

y(t) = arctan(2*exp(t)*_C1/(1+exp(2*t)*_C1^2), (-exp(2*t)*_C1^2+1)/(1+exp(2*t)*_C1^2))

c. Use dsolve to obtain the solution satisfying the initial condition y(0) = 1. What is the value of this solution when t = 1? Get the exact value and an approximation. If the solution is called "soln" use

eval(soln,t=1); evalf(%);

> soln := dsolve( {DE,y(0)=1} );

soln := y(t) = arctan(2*exp(t)*sin(1)*(-1+cos(1))/(-sin(1)^2-2*(exp(t))^2+2*(exp(t))^2*cos(1)+(exp(t))^2*sin(1)^2), -(2*(exp(t))^2*cos(1)+sin(1)^2+(exp(t))^2*sin(1)^2-2*(exp(t))^2)/(-sin(1)^2-2*(exp(t...

> eval(soln,t=1); evalf(%);

y(1) = -arctan(2*exp(1)*sin(1)*(-1+cos(1))/(2*(exp(1))^2*cos(1)+sin(1)^2+(exp(1))^2*sin(1)^2-2*(exp(1))^2))+Pi

y(1) = 1.956294970

d. Comment on the long-term behavior of solutions to this differential equation.

The slope field suggests that as t approaches infinity, all solutions have a horizontal asymptote. That is, the solutions have a finite limit (depending on the initial value) as t appoaches infinity. The arctan function in the solution formula confirms this observation. One can also appeal to the fact that the values y = k*Pi are stationary points on the phase line. All solutions are either attracted to k*Pi or repelled from k*Pi for some k.

Section 4. Approximate Solutions

7. Modify the ModEuler procedure to make a procedure that implements the classical Runge-Kutta algorithm (see Ledder, Section 2.6). Test in by applying it to the problems described in Exercises 2 - 6. (Copy and Paste).

> RK := proc(f,t0,y0,h,N)
 local n, k, T, Y;

   T[0] := t0;

   Y[0] := y0;

   for n from 0 to N-1

     do

       T[n+1] := T[n] + h:

       k[1] := f(T[n],Y[n]):

       k[2] := f(T[n] + 0.5*h,Y[n] + 0.5*k[1]*h):

       k[3] := f(T[n] + 0.5*h,Y[n] + 0.5*k[2]*h):

       k[4] := f(T[n+1],Y[n] + k[3]*h):

       Y[n+1] := Y[n] + 1/6*(k[1]+2*k[2]+2*k[3]+k[4])*h:

     end do:

   [ [T[m],Y[m]] $ m=0..N ]

end proc;

RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....RK := proc (f, t0, y0, h, N) local n, k, T, Y; T[0] := t0; Y[0] := y0; for n from 0 to N-1 do T[n+1] := T[n]+h; k[1] := f(T[n], Y[n]); k[2] := f(T[n]+.5*h, Y[n]+.5*k[1]*h); k[3] := f(T[n]+.5*h, Y[n]+....

Test on Exercise 5.

> DE := diff(y(t),t) = cos(t)*y(t);
f := (t,y) -> cos(t)*y;

DE := diff(y(t), t) = cos(t)*y(t)

f := proc (t, y) options operator, arrow; cos(t)*y end proc

> soln := dsolve( {DE, y(0)=1} );
g := unapply(rhs(soln),t);

soln := y(t) = exp(sin(t))

g := `@`(exp, proc (t) options operator, arrow; sin(t) end proc)

> L1 := RK(f,0,1,0.5,20):
plot( [L1,g(t),L1], t=0..10, y=0..3, style=[line$2,point],

      color=[red,blue,black]);

[Plot]

> L2 := RK(f,0,1,0.25,40):
plot( [L2,g(t),L2], t=0..10, y=0..3, style=[line$2,point],

      color=[red,blue,black]);

[Plot]

> Matrix([ [t,0.5*k$k=0..5],
        ['h=0.5',evalf[6]('L1[k,2]'$k=1..6)],

        ['h=0.25',evalf[6]('L2[2*k+1,2]'$k=0..5)],

        ['g(t)',evalf[6](g(0.5*k)$k=0..5)] ] );

Matrix([[t, 0., .5, 1.0, 1.5, 2.0, 2.5], [h = .5, 1., 1.61486, 2.31919, 2.71076, 2.48190, 1.81876], [h = .25, 1., 1.61513, 2.31974, 2.71144, 2.48254, 1.81931], [g(t), 1., 1.61515, 2.31978, 2.71148, 2....