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]
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 ]
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]
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]
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]=
Out[105]=
In[106]:=
Simplify[y[t]/.soln]
Out[106]=
In[107]:=
Qrate = Solve[ (y[t]/.soln/.t->5)==0, Q]
Out[107]=
In[108]:=
MonthlyRate = (Q/.Qrate)/12
Out[108]=
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
where P[1] is the amount that is put towards reducing the principal (the %12,000) in the first month. The second payment is
where P[2] is the amount that is put towards reducing the principal in the second month. The nth payment is
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]=
This determines P[2] in terms of P[1].
In[2]:=
Solve[ equation, P[2] ]
Out[2]=
In[3]:=
P[2] = P[2]/.%[[1]]//FullSimplify
Out[3]=
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
= 12000
In[13]:=
Solve[Sum[P[k], {k,60}] == 12000, P[1]]
Out[13]=
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]=
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]
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} ] ]
b. Obtain the general solution formula.
In[36]:=
DSolve[ y'[t]==Sin[y[t]], y[t], t]
Out[36]=
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]=
In[39]:=
y[t]/.soln/.t->1
N[%]
Out[39]=
Out[40]=
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 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]=
In[58]:=
g = DSolve[ {DE,y[0]==1}, y, t][[1,1,2]];
g[t]
Out[59]=
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" ]
Created by Mathematica (December 8, 2004)