Second and Higher Order Differential Equations

This tutorial contains many Mathematica scripts. You, as the user, are free to use all codes for your needs, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately. Any comments and/or contributions for this tutorial are welcome; you can send your remarks to <Vladimir_Dobrushkin@brown.edu>

Return to computing page
Return to computing page for the second course APMA0330
Return to Mathematica tutorial page for the first course
Return to Mathematica tutorial page for the second course
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340

Content

Part IV. Second and hifger order differential equations

Fundamental set of solutions. Wronskian
General solution
Reduction of order
Non-homogeneous equations.
Lagrange's method
Method of undetermined coefficients

Operator methods (not sure yet)
Numerical Solutions of second order ODEs
Applications

 

 

1.4.1. Fundamental Set of Solutions

 

To check linearly independence of two functions, we have two options. First, two functions are linearly independent if and only if one of them is a constant multiple of another. Another option is to calculate the Wronskian if you know that these two functions are solutions of the same differential equation.

Example. Given two functions \( \{ x , 3x \} ,\) determine whether they are linearly dependent.

y={x,3 x}
w={y,D[y,x]}
Det[w]
Out[15]= 0

I. Wronskian

The Mathematica command Wronskian[{y1 ,y2, ...},x] gives the Wronskian determinant for the functions \( y1 (x), y2 (x) , \ldots ,\) depending on x.

 

 

 

 

 

1.4.2. General Solutions of Homogeneous Equations

I. Verification

Example 4.4.1. We verify that the function \( y(x) = C_1 \,\cos 2x + C_2 \,\sin 2x + x^3 -3x/2 \) is a solution of the differential equation \( y'' + 4\,y = 4\,x^3 .\)

y[x_] = c1 Cos[2 x] + c2 Sin[2 x] + x^3 - 3*x/2;
Simplify[y''[x] + 4 y[x] == 4 x^3]
Out[2]= True
ics = {y[0] == 1, y'[0] == -1};
c1c2 = Solve[ics, {c1, c2}];
SolRule[x_] = y[x] /. c1c2[[1]]
SolRule''[x] + 4 SolRule[x] == 4 x^3 && SolRule[0] == 1 && SolRule'[0] == -1
Simplify[%]
Out[7]= True

To check linearly dependence:

Clear[y];
y={Cos[2 x], Sin[2 x]}
w={y,D[y,x]}
Det[w]
Out[11]= 2 Cos[2 x]^2 + 2 Sin[2 x]^2
y={x,3 x}
w={y,D[y,x]}
Det[w]
Out[15]= 0

Example 4.4.2. We consider the linear differential operator:

 

II. Linear Differential Operators

Clear[x, y];
L[x_,y_] =y''[x] -y'[x]- 6y[x]
CharPoly[lambda_] =Coefficient[L[x,Exp[lambda #]&], Exp[lambda x]]
roots =lambda/.Solve[CharPoly[lambda]==0,lambda]
Out[15]= {-2, 3}
Soln = Map[Exp[# x] &, roots]
AllSoln[x_]=Soln.{c1,c2}
Out[17]= c1 E^(-2 x) + c2 E^(3 x)
Simplify[L[x, AllSoln] == 0] (* check the answer *)
Out[18]= True
DSolve[y''[x] - y'[x] - 6 y[x] == 0, y[x], x]

Higher Order Equations

L[x_, y_] = y'''[x] - 2 y''[x] - 5 y'[x] + 6 y[x]
char[lambda_] =Coefficient[L[x,Function[t,Exp[lambda t]]],Exp[lambda x]]
roots = r /. Solve[char[r] == 0, r]
solns = Map[Function[k, Exp[k x]], roots]
y[x_] = solns.{c2,c1,c3} (* or *)
Clear[x,y];
L[x_, y_] = y'''[x] - 2 y''[x] - 5 y'[x] + 6 y[x]
DSolve[L[x, y] == 0, y[x], x]
y[x_] = Expand[y[x] /. %[[1]] ]
basis = Table[Coefficient[y[x], C[i]], {i, 1, 3}]
Factor[Coefficient[L[x, Function[t, Exp[r t]]], Exp[r x]]]
W[x_] = NestList[Function[t, D[t, x]], basis, 2]
Det[W[x]]
Out[7]= 30 E^(2 x) (* Wronskian *)

Composition of Linear Differential Operators:

L2[x_, y_] := a[2] D[y[x], {x, 2}] + a[1] D[y[x], {x, 1}] + a[0] y[x]
L3[x_, y_] := b[3] D[y[x], {x, 3}] + b[2] D[y[x], {x, 2}] + b[1] D[y[x], {x, 1}] + b[0] y[x]
L2L3 = Collect[L2[x, Function[z, L3[z, y]]], y]

 

III. Real Distinct Roots

 

 

 

 

IV. Complex Roots

ExpToTrig[Exp[I 2 x]]
Out[9]= Cos[2 x] + I Sin[2 x]
y={E^(a x) Cos[2x], E^(a x) Sin[2x]};
w = {y, D[y, x]};
Det[w] (* or *)
Simplify[Det[w], Trig -> True]
Out[22]= 2 E^(2 a x)
which is not zero. So these two functions are linearly independent.
DSolve[y''[x]-4 y'[x]+5 y[x]==0,y[x],x]
Out[1]= {{y[x] -> E^(2 x) C[2] Cos[x] + E^(2 x) C[1] Sin[x]}}
soln[x_] = Expand[y[x]/.%[[1]]]
Out[2]= E^(2 x) C[2] Cos[x] + E^(2 x) C[1] Sin[x]
y1[x_] = Coefficient[soln[x],C[1]]
Out[3]= E^(2 x) Sin[x]
y2[x_] = Coefficient[soln[x], C[2]]
Out[4]= E^(2 x) Cos[x]
Wronskian[{y1[x], y2[x]}, x] (* built-in Wronskian *)
Out[5]= -E^(4 x)
DSolve[{y''[x]-4 y'[x]+5 y[x]==0,y[0]==1},y[x],x]
Out[6]= {y[x] -> E^(2 x) (Cos[x] + C[1] Sin[x])}}

IVP:

DSolve[{y''[x] - 4 y'[x] + 5 y[x] == 0, y[0] == 1}, y[x], x];
soln[x_] = Expand[y[x] /. %[[1]]] /. C[1] -> c1
curves = Table[soln[x], {c1, -2, 2}];
Plot[Evaluate[curves], {x, -2, 2}, PlotRange -> {-10, 10}]
Out[11]= E^(2 x) Cos[x] + c1 E^(2 x) Sin[x]
Out[13]=

DSolve[{y''[x] - 4 y'[x] + 5 y[x] == 0, y'[0] == 1}, y[x], x];
soln[x_] = Expand[y[x] /. %[[1]]] /. C[1] -> c2
curves = Table[soln[x], {c2, -2, 2}];
Plot[Evaluate[curves], {x, -0.5, 2}, PlotRange -> {-10, 10}]
Out[15]= 1/2 E^(2 x) Cos[x] - 1/2 c2 E^(2 x) Cos[x] + c2 E^(2 x) Sin[x]
Out[17]=

Example.
Consider the second order differnetial equation \( y'' +y=0, \) which has two linearly independent solutions: \( \cos x, \sin x .\)

DSolve[y''[x] + y[x] == 0, y[x], x]
Out[14]= {{y[x] -> C[1] Cos[x] + C[2] Sin[x]}}
toplot1 =
Table[C[2] Cos[x] - C[1] Sin[x] /. {C[2] -> i, C[1] -> 0}, {i, -2.5, 2.5, 0.5}];
toplot2 =
Table[C[2] Cos[x] - C[1] Sin[x] /. {C[2] -> 0, C[1] -> i}, {i, -2.5, 2.5, 0.5}];
plot1 = Plot[Evaluate[toplot1], {x, 0, 4*Pi},
DisplayFunction -> Identity]

plot2 = Plot[Evaluate[toplot2], {x, 0, 4*Pi}, DisplayFunction -> Identity]

 

 

V. Boundary Value Problem


L[x_, y_] = y'''[x] - 2 y''[x] - 5 y'[x] + 6 y[x]
yb[x_] = y[x] /. DSolve[L[x, y] == 36 x, y[x], x][[1]]
const = Solve[{yb[0] == 2, yb[1] == 1, yb[2] == -1}, {C[1], C[2], C[3]}]
z[x_] = Simplify[yb[x] /. const[[1]]]
Out[5]= 5 + (E^(-1 + 3 x) (10 - 3 E - 18 E^2 + 10 E^3))/((-1 + E)^2 (1 + E) (1 + E +
E^2 + E^3 + E^4)) + ( E^(-1 + x) (-10 + 18 E^2 + 3 E^3 - 10 E^5))/((-1 + E)^2 (1 + E) (1 +
E + E^2)) + (E^(4 - 2 x) (-18 + 10 E + 10 E^3 - 3 E^4))/(
1 - E^3 - E^5 + E^8) + 6 x
{p1, p2, p3} = N[{C[1], C[2], C[3]} /. const[[1]]]
zb[x_] = yb[x] /. {C[1] -> p1, C[2] -> p2, C[3] -> p3}
Plot[zb[x], {x, -1, 2.5}]

 

1.4.3. Reduction of Order

I. Double Roots

DSolve[y''[x]-4 y'[x]+4 y[x]==0,y[x],x]
soln[x_] = Expand[y[x]/.%[[1]]]
y1[x_] = Coefficient[soln[x],C[1]]
y2[x_] = Coefficient[soln[x], C[2]]
basis = {y1[x], y2[x]}
WronskianDet = Det[{basis, D[basis, x]}]
Out[22]= E^(4 x)

soln = DSolve[{x''[t] + 6 x'[t] + 9 x[t] == 0, x[0] == -1,
x'[0] == 5}, x[t], t]
s2[t_] = x[t] /. soln[[1]]
Plot[s2[t], {t, 0.2, 3}]

Out[11]= {{x[t] -> E^(-3 t) (-1 + 2 t)}}
Out[12]= E^(-3 t) (-1 + 2 t)
Out[13]=


When does the extremum occur?
Solve[s2'[t] == 0, t]
out[48]= {{t -> 5/6}}

Where is the extremum?
s2[t /. %]
Out[49]= {2/(3 E^(5/2))}
N[%]
Out[50]= {0.0547233}

 

 

II. Reduction of Order

 

1.4.4. Variation of Parameters/Lagrange

I. Case Sensitivity

 

Suppose the fundamental set of solutions is \( x^3 \) and \( 1/x, \) and the right-hand side function is
\( (x-3)^2 /x , \) so we need to solve the non-homogeneous equation

\[ L[y]=(x-3)^2/x . \]
Let \( y_1=x^3 , \quad y_2=1/x \) be the fundamental set of solutions of the second order homogeneous differential equation
\[ y'' +p(x)y' +q(x)y=0 . \]

Then

\[ p(x) = - w'/w ; \qquad q(x) =(y'_1 y''_2 -y''_1 y'_2)/w , \] where \( w(x) \) is the Wronskian of \( y_1 , y_2 . \)

Now we determine the linear differential operator \( L[y] \) having \( y_1 (x) \ \mbox{and} \ y_2 (x) \) as solutions of \( L[y]=0. \)

y1[x_] = x^3
y2[x_] = 1/x
w[x_] = Wronskian[{y1[x], y2[x]}, x]
Out[3]= -4x
p[x_] = -w'[x]/w[x]
Out[4]= -(1/x)
q[x_] = (y1'[x] y2''[x] - y1''[x] y2'[x])/w[x]
Out[5]= -(3/x^2)
L[x_, y_] = y''[x] + p[x]*y'[x] + q[x]*y[x] (* oerator *)
Out[6]= -((3 y[x])/x^2) - Derivative[1][y][x]/x + (y^\[Prime]\[Prime])[x]

f[x_] = (x - 3)^2 /x;
basis = {x^3, 1/x}
W[x_] = Table[D[basis, {x, i}], {i, 0, 1}]// Flatten
Out[9]= {x^3, 1/x}, {3 x^2, -(1/x^2)}
F[x_] = {0, f[x]}
Out[10]= {0, (-3 + x)^2/x}
v[x_] = Integrate[Inverse[W[x]].F[x], x]
Out[11]= {1/4 (-(9/(2 x^2)) + 6/x + Log[x]),
1/4 (-((9 x^2)/2) + 2 x^3 - x^4/4)}
yp[x_] = Simplify[basis.v[x]] (* a particular solution *)
Out[12]= 1/16 x (-36 + 32 x - x^2 + 4 x^2 Log[x])
Simplify[L[x, yp] == f[x]] (* check solution *)
Out[13]= True

 

 

 

 

1.4.5. Method of Undetermined Coefficients

 

Clear[x,y];
L[x_, y_] = y'''[x] - 2 y''[x] - 5 y'[x] + 6 y[x]
rhs = x^3
L[x,y] == rhs
assume[x_] = a0 + a1*x + a2*x^2 + a3*x^3
result = Collect[L[x, assume], x]
eqns = CoefficientList[result, x] == CoefficientList[rhs, x]
yp[x_] = assume[x] /. Solve[eqns][[1]]
Out[40]= 209/216 + (37 x)/36 + (5 x^2)/12 + x^3/6

 

 

 

1.4.6. Operator Method

I. Case Sensitivity

 

 

 

 

1.4.7. Numerical Solutions of Second Order ODEs

 

Example 1.4.3:
Rayleigh's equation is a non-linear equation

\[ x'' + \mu \left( x^2 -1\right) x +x =0 , \]

that arises in the study of the motion of a violin string. We plot three solutions subject to three different initial conditions:

numsol = NDSolve[{x''[t] + ((1/3)*(x'[t])^2 - 1)*x'[t] + x[t] == 0, x[0] == 1,
x'[0] == 0}, x[t], {t, 0, 15}]
Out[1]= {{x[t] -> InterpolatingFunction[{{0., 15.}}],<>],[t]}}
numsol1 = NDSolve[{x''[t] + ((1/3)*(x'[t])^2 - 1)*x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], {t, 0, 15}]
Out[2]= {{x[t] -> InterpolatingFunction[{{0., 15.}}],<>],[t]}}
numsol2 = NDSolve[{x''[t] + ((1/3)*(x'[t])^2 - 1)*x'[t] + x[t] == 0, x[0] == 0.1, x'[0] == 0}, x[t], {t, 0, 15}]
Out[3]= {{x[t] -> InterpolatingFunction[{{0.,15.}}, <>][t]}}
To see three graphs, we type:
numsol3 = NDSolve[{x''[t] + ((1/3)*(x'[t])^2 - 1)*x'[t] + x[t] == 0, x[0] == 0, x'[0] == 2}, x[t], {t, 0, 15}]
Plot[Evaluate[x[t] /. {numsol, numsol2, numsol3}], {t, 0, 15},
PlotStyle -> {GrayLevel[0.1], Dashed, Thick}]

To identify each graph, it is helpful to use PlotLegend option:

Plot[Evaluate[x[t] /. {numsol, numsol2, numsol3}], {t, 0, 15}, PlotStyle -> {GrayLevel[0.1], Dashed, Thick},
PlotLegends -> Automatic]
Plot[Evaluate[x[t] /. {numsol, numsol2, numsol3}], {t, 0, 15}, PlotStyle -> {GrayLevel[0.1], Dashed, Thick},
PlotLegends -> {"IC: x(0)=1, x'(0)=0", "IC: x(0)=0.1, x'(0)=0", "IC: x(0)=0, x'(0)=2"}]

 

I. Pendulum Equation

There seems no doubt that the first person who investigated and established the mathematical theory and properties of the pendulum was

Christiaan Huygens (this spelling of his name is taken from the title of his 1658 book Horologium Oscillatorium), the Dutch

philosopher. The main point of Huygens discovery was that the curve in which a bob hang by a string of insensible weight must move in order to be isochronous in its vibrations, is not a circle, but a cycloid.

For a point pendulum supported by massless, inextensible cord of length \( \ell , \) the equation of motion for oscillations in a vacuum is

\[ \ddot{\theta} + \left( g/\ell \right) \sin \theta =0 , \]
where \( \ddot{\theta} = {\text d}^2 \theta / {\text d}t^2 \) and g is gravitational acceleration.

The resistance of the air acts on both the pendulum ball and the pendulum wire. It causes the amplitude to decrease with time and increases the period of oscillation slightly. The drag force is proportional to the velocity for values of the Reynolds number of the order 1 or less. For values of the Reynolds number of order \( 10^3 \) to \( 10^5 , \) the force is proportional to the square of the velocity. The maximum Reynolds number based on diameter for the ball is 1100, where the quadratic force law should apply, while the maximum value based on diameter of wire is about 6, where the linear force law should prevail.

Since the damping force is neither linear nor quadratic, but rather a combination of the two, it makes sense to establish a damping function which contains both effects simultaneously. Therefore, the pendulum equation becomes

\[ \ddot{\theta} + 2\alpha \,\dot{\theta} + \epsilon \,\mbox{sign} \left( \dot{\theta} \right) \dot{\theta}^2 + \omega_0^2 \sin\theta =0 , \]

where \( \alpha = \kappa /(2m) ,\quad \epsilon = c\ell /m, \quad \omega_0 = \left( g/\ell \right)^{1/2} ,\) and \( \mbox{sign}(x) = \begin{cases} 1, & \ x>0, \\ -1, & \ x<0. \end{cases} \)

In order to find out the relationship between the angle of the pendulum and time t, a good way is to plot \[Theta] against t.
In order to make the plot, first solve the differential equation with the command NDSolve.

This can be done by assigning constant values to \[Alpha], \[Epsilon], and \[Omega] when using the NDSolve command, get three
different functions in respect to three different \[Alpha] values, and use the Plot command to plot the three functions in one graph

numsol1 =
  NDSolve[{y''[x] == -10*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsol2 =
  NDSolve[{y''[x] == -1*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsol3 =
  NDSolve[{y''[x] == -0.1*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[y[x] /. {numsol1, numsol2, numsol3}], {x, 0, 100}, PlotRange -> All, PlotStyle -> {GrayLevel[0.1], Dashed, Thick}, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=10", "\[Alpha]=1", "\[Alpha]=0.1"}]

Another method is to set up a function that takes \[Theta]0, \[Omega]0,\[Alpha], \[Epsilon] and tmax as input and out put a function of \[Theta].

DampedDrivenPendulumFunction2[{\[Theta]0_, \[Omega]0_}, \[Alpha]_, \ \[Epsilon]_, tmax_] := Module[{\[Theta]}, \[Theta] /. First[ NDSolve[{ \[Theta]''[t] == -Sin[\[Theta][t]] - 2 \[Alpha] \[Theta]'[ t] - \[Epsilon] Sign[\[Theta]'[t]] (\[Theta]'[t])^2, \[Theta][0] == \[Theta]0, \[Theta]'[0] == \[Omega]0}, \[Theta], {t, 0, tmax}
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 0.1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 10, 0, 100][x]}, {x, 0, 100}, PlotLegends -> "Expressions", PlotRange -> Full, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]
The third method presented here is very different from the first and the second: we create a loop to add solutions into a list.
\[Epsilon] = 0; g = 9.8; \[Omega]0 = 1;
listofsolutions = {};
listofvaluesforalpha = {.1, 1, 10};
Numerically solve the corresponding second order differential equation, leaving alpha as a variable:
pendulumeq = \[Theta]''[t] + 2*\[Alpha] \[Theta]'[ t] + \[Epsilon] Sign[\[Theta]'[ t]]*((\[Theta]'[t])^2) + (\[Omega]0^2)*Sin[\[Theta][t]] == 0;
Get distinct equations for the different values of \[Alpha] and add them to the listofsolutions:
counter = 1;
While[counter <= Length[listofvaluesforalpha], AppendTo[listofsolutions, NDSolve[{Evaluate[ pendulumeq /. {\[Alpha] -> listofvaluesforalpha[[counter]]}], \[Theta][ 0] == .7, \[Theta]'[0] == 0}, \[Theta][t], {t, 0, 100}]]; counter++]
Plot[Evaluate[theta[t] /. listofsolutions], {t, 0, 50}, PlotRange -> All]

The three different methods all get the same result, from which we can see that the largest damping coefficient does not stop the pendulum the fastest. So the qualitative difference among the three diffrent coefficients is that when \[Alpha] is 0.1, the pendulum oscilates (cross 0). In other cases, the pendulum doesn' t oscilate (doesn' t cross 0).

The critical dammping value is the damping coefficient that stop the \ pendulum the fastest.

In order to find out the critical damping value of alpha, first we \ need to define the meaning of stop, because the pendulum will end up \ always be in very small oscilation.

The first method defined - 0.035 < \[Theta] < 0.035 as stop. By plotting different functions of \[Theta] with different alpha \ values (0.6, 0.7, 0.8, 0.9, 1.0) in one graph, you can see that \[Alpha] = 0.7 starts to be in the range [-0.035, 0.035] the fastest.

numsola1 =
  NDSolve[{y''[x] == -0.6*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolb1 =
  NDSolve[{y''[x] == -0.7*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolc1 =
  NDSolve[{y''[x] == -0.8*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsold1 =
  NDSolve[{y''[x] == -0.9*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolf1 =
  NDSolve[{y''[x] == -1.0*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[ y[x] /. {numsola1, numsolb1, numsolc1, numsold1, numsolf1}], {x, 3, 9}, PlotStyle -> {GrayLevel[0.2], Dashed, Thick}, PlotRange -> All, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=0.6", "\[Alpha]=0.7", "\[Alpha]=0.8", "\[Alpha]=0.9", "\[Alpha]=1.0"}]
To find out more accurate value of \[Alpha], plot several different \[Alpha] values close to 0.7:
numsola =
  NDSolve[{y''[x] == -0.74*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolb =
  NDSolve[{y''[x] == -0.73*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolc =
  NDSolve[{y''[x] == -0.72*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsold =
  NDSolve[{y''[x] == -0.71*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolf =
  NDSolve[{y''[x] == -0.70*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolg =
  NDSolve[{y''[x] == -0.69*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsolh =
  NDSolve[{y''[x] == -0.68*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
numsoli =
  NDSolve[{y''[x] == -0.67*2*y'[x] - Sin[y[x]], y[0] == 0.7, y'[0] == 0}, y[x], {x, 0, 1000}]
Plot[Evaluate[ y[x] /. {numsola, numsolb, numsolc, numsold, numsolf, numsolg, numsolh, numsoli}], {x, 3, 8}, PlotStyle -> {GrayLevel[0.1], Dashed, Thick}, PlotRange -> All, AxesLabel -> {"t", "\[Theta]"}, PlotLegends -> {"\[Alpha]=0.74", "\[Alpha]=0.73", "\[Alpha]=0.72", "\[Alpha]=0.71", "\[Alpha]=0.70", "\[Alpha]=0.69", "\[Alpha]=0.68", "\[Alpha]=0.67"}]
Then \[Alpha] = 0.69 should rouphly be the critical damping value.
Pi/2*0.005
Out[25]= 0.00785398
findmid[x_, y_] := (x + y)/2
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 1, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 0.7, 0, 100][x]}, {x, 0, 30}, PlotLegends -> "Expressions", PlotRange -> {-0.007853981633974483`, 0.007853981633974483`}, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]
From this graph, when \[Alpha] = 1, the pendulum settles in the stop range faster than when \[Alpha] = 0.7. This means that 0.7 is too small, and 1 may be too large. Then try the mid value 0.75.
findmid[1, 0.7]
Plot[{DampedDrivenPendulumFunction2[{0.7, 0}, 0.85, 0, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 1, 2, 100][x], DampedDrivenPendulumFunction2[{0.7, 0}, 0.7, 0, 100][x]}, {x, 0, 30}, PlotLegends -> "Expressions", PlotRange -> {-0.007853981633974483`, 0.007853981633974483`}, PlotStyle -> {Thick, Automatic, Dashed}, AspectRatio -> 1.0]

When \[Alpha] = 0.85, the pendulum settles in the range faster than when \[Alpha] = 1. So the critical value may be smaller than 0.85. Try the mid value of 0.85 and 0.7. Finally, the answer got from this method is 0.81.

The third method defines stop as not passing 0 in the time span 0 < t \ < 500, uses a function "piano" to check if the pendulum oscilates.

\[CurlyEpsilon] = 0; g = 9.8;
\[CapitalOmega]0 = 1;
listofsolutions = {};
listofvaluesfor\[Alpha] = Range[0, 2, .0005];
pendulumeq = \[Theta]''[t] + 2*\[Alpha]*\[Theta]'[t] + \[CurlyEpsilon]* Sign[\[Theta]'[t]]*((\[Theta]'[t])^2) + (\[CapitalOmega]0^2)* Sin[\[Theta][t]] == 0;
counter = 1;
While[counter <= Length[listofvaluesfor\[Alpha]], AppendTo[listofsolutions, NDSolve[{Evaluate[ pendulumeq /. {\[Alpha] -> listofvaluesfor\[Alpha][[counter]]}], \[Theta][ 0] == .7, \[Theta]'[0] == 0}, \[Theta][t], {t, 0, 1000}]]; counter++];
evaluatedlistofsolutions = Evaluate[\[Theta][t] /. listofsolutions];
Piano is a function that checks to see if the pendulum ever swings past \[Theta] = 0 for each of the functions in the list.If the answer is no, then it means that the pendulum no longer oscillates for that value of \[Alpha].
piano[key_] :=
  Catch[If[FindMinValue[ evaluatedlistofsolutions[[key]], {t, 0, 500}] > 0, Throw[key], Throw[906]]];
counter2 = 1;
While[piano[counter2] == 906, counter2++];
Print["The Value of \[Alpha] at which the function \[Theta][t] stops \ oscillating is ", listofvaluesfor\[Alpha][[counter2]]]
The Value of \[Alpha] at which the function \[Theta][t] stops oscillating is 0.8985
Plot[evaluatedlistofsolutions[[piano[counter2]]], {t, 0, 20}, PlotRange -> All]
Now we assume that alpha is not zero. .
DampedDrivenPendulumFunction2[{\[Theta]0_, \[Omega]0_}, \[Alpha]_, \ \[Epsilon]_, tmax_] :=
  Module[{\[Theta]}, \[Theta] /. First[ NDSolve[{ \[Theta]''[t] == -Sin[\[Theta][t]] - 2 \[Alpha] \[Theta]'[ t] - \[Epsilon] Sign[\[Theta]'[t]] (\[Theta]'[t])^2, \[Theta][0] == \[Theta]0, \[Theta]'[0] == \[Omega]0}, \[Theta], {t, 0, tmax}] ]]
g5 = DampedDrivenPendulumFunction2[{0.7, 0}, 0.1, 1, 100]
g6 = DampedDrivenPendulumFunction2[{0.7, 0}, 10, 1, 100]
g7 = DampedDrivenPendulumFunction2[{0.7, 0}, 1, 1, 100]
Plot[{g5[x], g6[x], g7[x]}, {x, 0, 100}, PlotLegends -> "Expressions", PlotRange -> Full, PlotStyle -> {Thick, Automatic, Dashed}]
When \[Epsilon] is 1, the qualitative difference among different \ alpha when apha is 1, 0.1, or 10 is the same. But when alpha is 0.1, it oscilates with less times compared to the case when epsilon is 0.

 

 

II. Spring Equation

Consider the following spring equation with aging:

\[ \ddot{x} + 4\,e^{-t/4}\,x(t) =0 \qquad x(0) =1, \quad \dot{x} (0)=2 . \]

We solve it and plot in the following sequence of steps.

sol = DSolve[{x''[t] + 4*Exp[-t/4]*x[t] == 0, x[0] == 1, x'[0] == 2},
x, t] // InputForm
Out[1]=
{{x -> Function[{t}, (BesselJ[0, 16*Sqrt[E^(-t/4)]]*
BesselY[0, 16] - BesselJ[0, 16]*BesselY[0,
16*Sqrt[E^(-t/4)]] + BesselJ[1, 16]*
BesselY[0, 16*Sqrt[E^(-t/4)]] -
BesselJ[0, 16*Sqrt[E^(-t/4)]]*BesselY[1, 16])/
(BesselJ[1, 16]*BesselY[0, 16] - BesselJ[0, 16]*
BesselY[1, 16])]}}
Plot[x[t] /. sol[[1]], {t, 0, 30}]
Wrap sol =... in a bracket for clarity to obtain the same result:
(sol = DSolve[{x''[t] + 4*Exp[-t/4]*x[t] == 0, x[0] == 1, x'[0] == 2},
x, t]) // InputForm

f[t_] = x[t] /. sol[[1]]
f[0.5]

Out[5]= {1.40733}
f[0.5][[1]]
Out[6]= 1.40733
This gives 1.40733 without the curly brackets.

Another way to remove curly brackets:

f[t_] = x[t] /. sol[[1]][[1]]
Out[7]=
(BesselJ[0, 16 Sqrt[E^(-t/4)]] BesselY[0, 16] -
BesselJ[0, 16] BesselY[0, 16 Sqrt[E^(-t/4)]] +
BesselJ[1, 16] BesselY[0, 16 Sqrt[E^(-t/4)]] -
BesselJ[0, 16 Sqrt[E^(-t/4)]] BesselY[1, 16])/(BesselJ[1,
16] BesselY[0, 16] - BesselJ[0, 16] BesselY[1, 16])
f[0.5]
Out[8]= 1.40733

 

 

 

 


1.4.9. Applications

 

We present some applications of differential equations of order greater than one.

Example 4.9.1: Flow Problems
A 500 liter container initially contains 10 kg of salt. A brine mixture of 100 gramms of salt per liter is entering the container at 6 liter per minute. The well-mixed contents
are being discharged from the tank at the rate of 6 liters per minute.

Express the amount of salt in the container as a function of time.

Salt is coming at the rate: 6*(0.1)=0.6 kg/min
d yin /dt =0.6 ; d yout/dt = 6x/500

dx/dt = 0.6 -6x/500 ; x(0)=10

DSolve[{x'[t]==6/10-6*x[t]/500, x[0]==10},x[t],t]
Out[15]= {{x[t] -> 10 E^(-3 t/250) (-4 + 5 E^(3 t/250))}}

Suppose that the rate of discharge is reduced to 5 liters per minute.

DSolve[{x'[t]==6/10-5*x[t]/(500+t), x[0]==10},x[t],t]
SolRule[t_] = Apart[x[t] /. First[%]]
Out[1]= {{x[t] -> (1/( 10 (500 + t)^5))(3125000000000000 + 187500000000000 t + 937500000000 t^2 + 2500000000 t^3 + 3750000 t^4 + 3000 t^5 + t^6)}}
Out[2]= 50 + t/10 - 1250000000000000/(500 + t)^5
To check the answer:
Simplify[SolRule'[t] == 6/10 - 5*SolRule[t]/(500 + t)]
Together[SolRule'[t] == Together[6/10 - 5*SolRule[t]/(500 + t)]]
SolRule[0]
Out[3]= True
Out[4]= True
Out[5]= 10

Example 4.9.2: Pure oscillations give the solutions of the following differential equation subject to the given initial conditions (displacement is 1, but velocity is zero):

\[ \ddot{x} (t) + 25\,x(t) =0 , \qquad x(0) =1, \quad \dot{x}(0) =0. \]

Clear[x,y,t,z];
soln = DSolve[{x''[t] + 25 x[t] == 0, x[0] == 1, x'[0] == 0}, x[t], t]
Plot[x[t] /. soln, {t, -1, 2.5}]
s[t_] = x[t]/.soln[[1]]
Plot[s[t],{t,0,3},AxesLabel->{"t","Displacement"}]

Example 4.9.3: Consider the initial value problem

\[ \ddot{x} + 17\,\dot{x} + 16\,x(t) =0, \qquad x(0) =1, \quad \dot{x} (0) =5. \]
soln = DSolve[{x''[t] + 17 x'[t] + 16 x[t] == 0, x[0] == 1, x'[0] == 5}, x[t], t];
s1[t_] = x[t] /. soln[[1]]
Out[32]= 1/5 E^(-16 t) (-2 + 7 E^(15 t))
Plot[s1[t], {t, 0, 3}, AxesLabel -> {"t", "Displacement"}]

When does the maximum excursion occur?

FindRoot[s1'[t] == 0, {t, 0}]
Out[9]= {t -> 0.101322}
What is the maximum excursion?
s1[t /. %]
Out[10]= 1.18603

Example 4.9.4: Now we turn to a differential equation whose characteristic polynomial has a double root:

\[ \ddot{x} + 6\,\dot{x} + 9\,x =0, \qquad x(0) =-1, \quad \dot{x}(0) =5. \]

soln = DSolve[{x''[t] + 6 x'[t] + 9 x[t] == 0, x[0] == -1, x'[0] == 5}, x[t], t]
s2[t_] = x[t] /. soln[[1]]
Plot[s2[t], {t, 0.2, 3}, PlotStyle->{Thick,Magenta}]

Out[11]= {{x[t] -> E^(-3 t) (-1 + 2 t)}}
Out[12]= E^(-3 t) (-1 + 2 t)
Out[13]=

When does the extremum occur?
Solve[s2'[t] == 0, t]
out[48]= {{t -> 5/6}}
Where is the extremum?
s2[t /. %]
out[49]= {2/(3 E^(5/2))}
N[%]
out[50]= {0.0547233}

Example 4.9.5: We consider another initial value problem

\[ \ddot{x} + 2\,\dot{x} + 37\,x =0, \qquad x(0) =-1, \quad \dot{x} (0) =5. \]
Since the characteristic polynomial has complex roots \( \lambda = -1 \pm 6{\bf j} ,\) the solutins decay with oscillations.
soln = DSolve[{x''[t] + 2 x'[t] + 37 x[t] == 0, x[0] == -1,
x'[0] == 5}, x[t], t]
Out[2]= -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])}}
s3[t_] = x[t] /. soln[[1]]

Clear[x,t];
soln = DSolve[{x''[t] + 2 x'[t] + 37 x[t] == 0, x[0] == -1, x'[0] == 5}, x[t], t]
s3[t_] = x[t] /. soln[[1]]
Plot[s3[t], {t, 0, 3}, PlotRange -> {-1, 1}]

Out[19]= {{x[t] -> -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])}}
Out[20]= -(1/3) E^-t (3 Cos[6 t] - 2 Sin[6 t])
Out[21]=

Now we define the amplitude of these oscillations:
amplitude = s3[t] /. c3_ Exp[c4_] (c1_ Cos[a_] + c2_ Sin[a_]) -> c3*Sqrt[c1^2 + c2^2]
Out[19]= -(Sqrt[13]/3)
Plot[{amplitude*Exp[-t], -amplitude*Exp[-t], s3[t]}, {t, 0, 3}, PlotStyle -> {Blue, Blue, Black}]

Example 4.9.6: Forced oscillations are modelled by the following nonhomogeneous equation:

\[ \ddot{x} + 10\,\dot{x} + 41\,x =25\,\sin 4t, \qquad x(0) =-1, \quad \dot{x}(0) =5. \]

soln = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 25 Sin[4 t], x[0] == -1, x'[0] == 5}, x[t], t]
Out[1]= {{x[t] ->
1/712 E^(-5 t) (-392 Cos[4 t] - 445 E^(5 t) Cos[4 t] +
125 E^(5 t) Cos[4 t] Cos[8 t] + 200 Sin[4 t] -
200 E^(5 t) Cos[8 t] Sin[4 t] + 200 E^(5 t) Cos[4 t] Sin[8 t] +
125 E^(5 t) Sin[4 t] Sin[8 t])}}

We extract the transient solution and the steady-state solution. First, we define the solution function:

ss[t_] = x[t] /. soln[[1]]
Out[2]= 5/712 E^(-5 t) (64 Cos[4 t] - 89 E^(5 t) Cos[4 t] + 25 E^(5 t) Cos[4 t] Cos[8 t] + 40 Sin[4 t] -
40 E^(5 t) Cos[8 t] Sin[4 t] + 40 E^(5 t) Cos[4 t] Sin[8 t] + 25 E^(5 t) Sin[4 t] Sin[8 t])
Then we find the steady-state solution by finding the general solution of the given differential equation
gensol = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 25 Sin[4 t]}, x[t], t]
Out[3]= {{x[t] -> E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] - 5/712 (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] +
40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])}}
gs[t_] = x[t] /. gensol[[1]]
Out[4]= E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] - 5/712 (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] +
40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])
and then vanishing arbitrary constants:
steadyState[t_] = gs[t] /. {C[1] -> 0, C[2] -> 0}
Out[5]= -(5/712) (89 Cos[4 t] - 25 Cos[4 t] Cos[8 t] + 40 Cos[8 t] Sin[4 t] - 40 Cos[4 t] Sin[8 t] - 25 Sin[4 t] Sin[8 t])
Then the transient solution will be their diference:
transient[t_] = Simplify[ss[t] - steadyState[t]]
Out[6]= 1/89 E^(-5 t) (-49 Cos[4 t] + 25 Sin[4 t])
Then we check the initial conditions:
steadyState[0]
Out[7]= -(40/89)
D[steadyState[t], t] /. t -> 0
Out[8]= 100/89
transient[0]
Out[9]= -(49/89)
D[transient[t], t] /. t -> 0
Out[10]= 345/89
ss[0]
Out[11]= -1
D[ss[t], t] /. t -> 0
Out[12]= 5

Now we consider a similar problem but with another input function:

Clear[soln,L,x]
L[t_,x_] = x''[t] + 10 x'[t] + 41 x[t];
solnRule = DSolve[{L[t, x] == 8 *Exp[-5*t] Cos[2 t], x[0] == -1, x'[0] == 5}, x[t], t]
Out[3]= {{x[t] ->
1/6 E^(-5 t) (-13 Cos[4 t] + 6 Cos[t]^2 Cos[4 t] +
Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])}}
s[t_] = x[t] /. solnRule[[1]]
Out[4]= 1/6 E^(-5 t) (-13 Cos[4 t] + 6 Cos[t]^2 Cos[4 t] +
Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])
Then we find the general solution
gensol = DSolve[{x''[t] + 10 x'[t] + 41 x[t] == 8*Exp[-5*t] Cos[2 t]}, x[t], t]
Out[5]= {{x[t] ->
E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] +
1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] +
3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])}}
gs[t_] = x[t] /. gensol[[1]]
Out[6]= E^(-5 t) C[2] Cos[4 t] + E^(-5 t) C[1] Sin[4 t] +
1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] +
3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])
Now the steady-state solution is
steadyState[t_] = gs[t] /. {C[1] -> 0, C[2] -> 0}
Out[6]= 1/6 E^(-5 t) (6 Cos[t]^2 Cos[4 t] + Cos[4 t] Cos[6 t] + 3 Sin[2 t] Sin[4 t] + Sin[4 t] Sin[6 t])

and the transient solution becomes

transient[t_] = Simplify[s[t] - steadyState[t]]
Out[7]= -(13/6) E^(-5 t) Cos[4 t]
Now we plot two transient solutions:
Plot[{1/89 E^(-5 t) (-49 Cos[4 t] + 25 Sin[4 t]) , transient[t]},{t,0,3}, PlotStyle->{Thick,{Blue,Pink}}]

We also plot two solutions

Plot[{ss[t], s[t]}, {t, 0, 2}, PlotStyle -> {Thick, {Cyan, Purple}}]

Transient solutions could be plotted ina similar way.

Plot[x[t] - steadyState[t], {t, 0, 1.5}, PlotRange -> {-1/2, 0.2}]

Example 4.9.7: Undamped Spring with External Forcing.

Enter the coefficients and the differential equation:

F0 = 1; omega = 2; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
Solve the initial value problem:
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->
1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])}}
Copy the output from the preceding step to define the solution function soln[t]:
soln[t_] = y[t] /. %[[1]]
Out[4]= 1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])
Plot the solution function.
Plot[soln[t], {t, 0, 10}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Pink}]
Define the derivative function as v[t], and plot the trajectory in the phase plane. (You may need to enlarge the plot to see details.)

Example 4.9.8: Beats

We reconsider Example 1.4.7 for periodic input, when resonance is observed.
Enter the coefficients and the differential equation.

F0 = 1; omega = 4; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
Solve the initial value problem:
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->
1/90 (-10 Cos[5 t] + 9 Cos[t] Cos[5 t] + Cos[5 t] Cos[9 t] +
9 Sin[t] Sin[5 t] + Sin[5 t] Sin[9 t])}}
Copy the output from the preceding step to define the solution function y[t]:

y[t_] = y[t] /. %[[1]]
Out[4]= 1/90 (-10 Cos[5 t] + 9 Cos[t] Cos[5 t] + Cos[5 t] Cos[9 t] +
9 Sin[t] Sin[5 t] + Sin[5 t] Sin[9 t])
Plot the solution function:
Plot[y[t], {t, 0, 15}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Green}]

Define the derivative function as v[t], and plot the trajectory in the phase plane:

v[t_] = D[y[t], t]
Out[6]= 1/90 (36 Cos[5 t] Sin[t] + 50 Sin[5 t] - 36 Cos[t] Sin[5 t] +
4 Cos[9 t] Sin[5 t] - 4 Cos[5 t] Sin[9 t])
ParametricPlot[{y[t], v[t]}, {t, 0, 15}, AspectRatio -> 1, AxesLabel -> {"y[t]", "v[t]"}, PlotStyle->{Thick,Gray}]

Example 4.9.9: Resonance

We repeat the commands (after quitting the kernel) here for defining the coefficients and the differential equation---now with \[omega] very close to k:

F0 = 1; omega = 4.9; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
Out[3]= {{y[t] ->-1.0101 Cos[5. t] + 1. Cos[0.1 t] Cos[5. t] +
0.010101 Cos[5. t] Cos[9.9 t] + 1. Sin[0.1 t] Sin[5. t] +
0.010101 Sin[5. t] Sin[9.9 t] }}
y[t_] = y[t] /. %[[1]]
Out[4]= -1.0101 Cos[5. t] + 1.0101 Cos[4.9 t] Cos[5. t]^2 +
1.0101 Cos[4.9 t] Sin[5. t]^2
Plot[y[t], {t, 0, 15}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Cyan}]
Now we find the velocity
v[t_] = D[y[t], t]
Out[6]= 0. - 4.94949 Cos[5. t]^2 Sin[4.9 t] + 5.05051 Sin[5. t] -
4.94949 Sin[4.9 t] Sin[5. t]^2
and plot v versus y:
ParametricPlot[{y[t], v[t]}, {t, 0, 15}, AspectRatio -> 1, AxesLabel -> {"y[t]", "v[t]"}, PlotStyle->{Thick,Orange}]

F0 = 1; omega = 2; k = 5;
ode = y''[t] + k^2*y[t] == F0*Cos[omega*t];
DSolve[{ode, y[0] == 0, y'[0] == 0}, y[t], t]
y[t_] = y[t] /. %[[1]]
Out[3]= {{y[t] ->
1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])}}
Out[4]= 1/210 (-10 Cos[5 t] + 7 Cos[3 t] Cos[5 t] + 3 Cos[5 t] Cos[7 t] +
7 Sin[3 t] Sin[5 t] + 3 Sin[5 t] Sin[7 t])
Plot[y[t], {t, 0, 10}, AxesLabel -> {"t", "y[t]"}, PlotStyle->{Thick,Dashed}]

v[t_] = D[y[t], t]
Out[6]= 1/210 (14 Cos[5 t] Sin[3 t] + 50 Sin[5 t] - 14 Cos[3 t] Sin[5 t] + 6 Cos[7 t] Sin[5 t] - 6 Cos[5 t] Sin[7 t])
ParametricPlot[{y[t], v[t]}, {t, 0, 10}, AspectRatio -> 1, PlotStyle->{Thick,Brown}, AxesLabel -> {"y[t]", "v[t]"}, PlotLabel -> "v[t] vesus y[t]"]

 

 


 

 

1.4.8. References

Martha L. Abell, James P. Braselton
Differential Equations With Mathematica,
Academic Press, 2004
ISBN 0120415623, 978012041562

Alfred Gray, Mike Mezzino, and Mark Pinsky
Introduction to Ordinary Differential Equations with Mathematica: An Integrated Multimedia Approach
Springer
ISBN-10: 0387944818 | ISBN-13: 978-0387944814

Brian R. Hunt, Ronald L. Lipsman, John E. Osborn, Donald A. Outing, Jonathan Rosenberg
Differential Equations with Mathematica, 3rd Edition
ISBN: 978-0-471-77316-0

Clay C. Ross
Differential Equations
An introduction with Mathematica
Springer-Verlag, 1995
ISBN: 0-387-94301-3