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 general solution.

a.   Plot some solutions starting at points evenly spaced on the y axis.   
b.  Use plot some solutions starting at points evenly spaced on the t axis.  
c.   Plot 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.

Hint. Use Table to generate a list of solutions, then plot them. See the example below.

In[87]:=

DE = y'[t] + y[t] == Cos[t];
solna = Table[ DSolve[{DE,y[0]==y0},y[t],t], {y0,-2,2}];
Show[ ListPlot[ Table[{0,k},{k,-2,2}], PlotStyle->PointSize[0.02]],
      Plot[ Evaluate[Table[y[t]/.solna[[k]],{k,5}]], {t,0,12}], AspectRatio→2/3]

[Graphics:../HTMLFiles/Solutions_61.gif]

Solution to b.

In[82]:=

solnb = Table[ DSolve[{DE,y[t0]==0},y[t],t], {t0,-2,2}];
Show[ ListPlot[ Table[{k,0},{k,-2,2}], PlotStyle->PointSize[0.02]],
      Table[Plot[ y[t]/.solnb[[k]], {t,k-3,12}],{k,5}], AspectRatio→1/3 ]

[Graphics:../HTMLFiles/Solutions_62.gif]

Solutions to c.

In[76]:=

solnc1 = Table[ DSolve[{DE,y[Cos[k*Pi/6]]==Sin[k*Pi/6]},y[t],t], {k,12}];
Show[ ListPlot[ Table[{Cos[k*Pi/6],Sin[k*Pi/6]},{k,12}],
                PlotStyle->PointSize[0.02]],
      Table[Plot[ y[t]/.solnc1[[k]], {t,Cos[k*Pi/6],2}],{k,12}] ,
      PlotRange->{{-2,2},{-2,2}}, AspectRatio→1/1]

[Graphics:../HTMLFiles/Solutions_63.gif]

In[78]:=

solnc2 = Table[ DSolve[{DE,y[Cos[k*Pi/6]]==Sin[k*Pi/6]},y[t],t], {k,12}];
Show[ ListPlot[ Table[{Cos[k*Pi/6],Sin[k*Pi/6]},{k,12}],
                PlotStyle->PointSize[0.02]],
      Table[Plot[ y[t]/.solnc1[[k]], {t,-2,Cos[k*Pi/6]}],{k,12}] ,
      PlotRange->{{-2,2},{-2,2}}, AspectRatio→1/1]

[Graphics:../HTMLFiles/Solutions_64.gif]

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.   

Let y(t) be the amount owed at time t and Q be the yearly rate of pay back.

In[104]:=

DE = y'[t] == 0.05 y[t] - Q
soln = DSolve[ {DE, y[0]==12000} , y[t], t]

Out[104]=

y^′[t]  -Q + 0.05 y[t]

Out[105]=

{{y[t] 2.71828^(-0.05 t) (12000. 7.38906^(0.05 t) - 20. 7.38906^(0.05 t) Q + 20. ^(0.05 t) Q)}}

In[106]:=

Simplify[y[t]/.soln]

Out[106]=

{^(0.05 t) (12000. - 20. Q) + 20. Q}

In[107]:=

Qrate = Solve[ (y[t]/.soln/.t->5)==0, Q]

Out[107]=

{{Q2712.49}}

In[108]:=

MonthlyRate = (Q/.Qrate)/12

Out[108]=

{226.041}

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 (e^0.05/12 - 1) + P[1]

     where P[1] is the amount that is put towards reducing the principal (the %12,000) in the first month. The second payment is

(12000 - P[1]) (e^0.05/12 - 1) + P[2]

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

(12000 - Underoverscript[∑, k = 1, arg3] P[k]) (e^0.05/12 - 1) + P[n]

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

In[1]:=

equation = 12000(Exp[0.05/12]-1)+P[1] ==
           (12000-P[1])*(Exp[0.05/12]-1) + P[2]

Out[1]=

50.1043 + P[1] 0.00417536 (12000 - P[1]) + P[2]

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

In[2]:=

Solve[ equation, P[2] ]

Out[2]=

{{P[2]  -1. (-50.1043 + 0.00417536 (12000. - 1. P[1]) - 1. P[1])}}

In[3]:=

P[2] = P[2]/.%[[1]]//FullSimplify

Out[3]=

0. + 1.00418 P[1]

So, we do this over and over again to define P[3], P[4], ... , P[60] in terms of P[1]. See the Help page for the Do function.

In[6]:=

Do[
eqn[n] = (12000 - Sum[P[k],{k,n-2}])(Exp[0.05/12]-1)+P[n-1]==
         (12000 - Sum[P[k],{k,n-1}])(Exp[0.05/12]-1)+P[n];
P[n] = FullSimplify[P[n]/.Solve[eqn[n],P[n]][[1]]],
{n,3,60} ]

Find P[1] using the fact that

Underoverscript[∑, k = 1, arg3] P[k]= 12000

In[13]:=

Solve[Sum[P[k], {k,60}] == 12000, P[1]]

Out[13]=

{{P[1] 176.408}}

The first payment will be this amount plus the interest on 12000 (for one month).

In[14]:=

12000(Exp[0.05/12]-1) + P[1]/.%

Out[14]=

{226.512}

Section 3. Slope Fields

In[15]:=

<<Graphics`PlotField`

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

In[30]:=

u = 1; v = Sin[y];
vf = PlotVectorField[ {u,v}/Sqrt[u^2+v^2], {t,-6,6}, {y,-6,6}, Axes->True]

[Graphics:../HTMLFiles/Solutions_79.gif]

a.   Put some solution curves into the plot using DSolve and Table. Comment  on the relationship between one curve and the next.  Instead of DSolve the solution curves are obtained using Mathematica's numeric solver function, NDSolve. This function is discussed in the next section.

In[34]:=

solns = Table[ NDSolve[{y'[t]==Sin[y[t]],y[0]==k},y[t],{t,-6,6}], {k,-6,6} ];
Show[ vf, Table[ Plot[y[t]/.solns[[k]], {t,-6,6}, PlotStyle->RGBColor[1,0,0]], {k,13} ] ]

[Graphics:../HTMLFiles/Solutions_80.gif]

b.   Obtain the general solution formula.

In[36]:=

DSolve[ y'[t]==Sin[y[t]], y[t], t]

Out[36]=

{{y[t] 2 ArcCot[^(-t - C[1])]}}

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.  

In[38]:=

soln = DSolve[ {y'[t]==Sin[y[t]], y[0]==1}, y[t], t]

Out[38]=

{{y[t] 2 ArcCot[^(-t) Cot[1/2]]}}

In[39]:=

y[t]/.soln/.t->1
N[%]

Out[39]=

{2 ArcCot[Cot[1/2]/]}

Out[40]=

{1.95629}

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

The direction 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 approaches infinity. The ArcCot function in the solution formula confirms this observation. One can also appeal to the fact that the values y = k π are stationary points on the phase line. All solutions are either attracted to  k π or repelled from k π  for  some integer k.

Section 4. Approximate Solutions

7.* Modify the ModEuler procedure to make a modular function 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). Is the error cut in half in problem 4?

In[51]:=

RK[f_,t0_,y0_,h_,N_] :=
Module[
       {T,K1,K2,K3,K4,Y},
        T[0] = t0; Y[0] = y0;
        T[n_]  :=  T[n] = T[n-1] + h;
        K1[n_] :=  K1[n] = f[T[n-1],Y[n-1]];
        K2[n_] :=  K2[n] = f[T[n-1] + 0.5h,Y[n-1] + 0.5K1[n]*h];
        K3[n_] :=  K3[n] = f[T[n-1] + 0.5h,Y[n-1] + 0.5K2[n]*h];
        K4[n_] :=  K4[n] = f[T[n],Y[n-1] + K3[n]*h];
        Y[n_]  :=  Y[n] = Y[n-1] + 1/6*(K1[n]+2K2[n]+2K3[n]+K4[n])h;
        Table[ {T[n],Y[n]}, {n,0,N}]
      ]

Test the algorithm on the IVP in the manual (Exercise 5).

In[48]:=

DE = y'[t] == Cos[t]y[t]
f[t_,y_] := Cos[t]y

Out[48]=

y^′[t] Cos[t] y[t]

In[58]:=

g = DSolve[ {DE,y[0]==1}, y, t][[1,1,2]];
g[t]

Out[59]=

^Sin[t]

In[68]:=

approx = RK[f,0,1,0.5,20];
pts = ListPlot[ approx, PlotStyle->PointSize[0.01], PlotRange->{0,3} ];
lines = ListPlot[ approx, PlotJoined->True, PlotStyle->RGBColor[1,0,0],  
                 PlotRange->{0,3} ];
Show[ pts, lines, Plot[ g[t], {t,0,10} ],
      AspectRatio->1/3, PlotLabel->"h = 0.5, Runge-Kutta" ]

[Graphics:../HTMLFiles/Solutions_88.gif]


Created by Mathematica  (December 8, 2004)