Differential equations

Classification of differential equations

Differential equations (DEs) form the basis of physics. Every physical process evolving in time, within classical of quantum mechanics,  is described by a DE. Also many time independent physical situations are describable in terms of DEs. Examples are problems of electrostatics and magnetostatics, theory of elasticity, stationary states in quantum mechanics.

Differential equations are equations for functions that contain their derivatives, such as

(x[t])'==f[t,x[t]].

DE containing derivatives with respect to only one variable, as above, are called ordinary differential equations (ODEs). If there are derivatives of multivariate functions over different variables in the equation, it is called partial differential equation (PDE). In this chapter we study ODEs.

The highest order of derivatives in a DE defines the order of the DE. Above is the first-order DE. Higher-order DE can be rewritten in a form of systems of first-order differential equations by considering lower-order derivatives as additional functions. For instance, the second-order DE

(x[t])''==f[t,(x[t])',x[t]]

can be rewritten as the system of two first-order DE

(v[t])'==f[t,v[t],x[t]]
(x[t])'==v[t].

Combined order of a system of DEs is the sum of the orders of the individual equations (2 in the example above). Rewriting a differential equation as a system of first-order DEs does not change the combined order.

DEs can be resolved or unresolved with respect to the highest derivative. Above are resolved DEs. Unresolved DEs can be writen as F( t,(x[t])',x[t])==0 etc. In some cases these equations can be resolved with respect to the highest derivative (such as ((((x[t])'))^2) == f[t, x[t]]), and in some cases not. It is difficult to find a physical model described by a DE non-resolvable with respect to the highest derivative.

DEs can be linear and nonlinear with respect to derivatives and the function. An example of a linear ODE is

(x[t])''+b (x[t])'+c x[t]==f[t].

Here the function on the right is called "inhomogeneous term". If the latter is zero, the linear ODE is callled "homogeneous". ODEs with constant coefficients (as a and b above) can be solved analytically (see below). ODEs with variable coefficients (functions of the independent variable) usually cannot be solved analytically. In some cases such ODEs define a special function as their solution, if these equations are important. An important property of homogeneous linear equations is that if (Subscript[x, 1])[t] and (Subscript[x, 2])[t] are solutions of the equation, then also (Subscript[C, 1])(Subscript[x, 1])[t]+(Subscript[C, 2])(Subscript[x, 2])[t] is a solution. The solution of an inhomogeneous linear DE can be represented as a particular solution of this inhomogeneous equation plus a general solution of the corresponding homogeneous equation that contains integration constants.

Vector differential equations

(r[t])'==f[t,r[t]]

etc. are in fact systems of DEs for components

(x[t])'==(Subscript[f, x])[t,x[t],y[t],z[t]]
(y[t])'==(Subscript[f, y])[t,x[t],y[t],z[t]]
(z[t])'==(Subscript[f, z])[t,x[t],y[t],z[t]].

In Mathematica, one writes systems of equations in a vector form by putting equations in a List.

Boundary conditions for differential equations

General solution of a system of DEs of a combined order N depends of N integration constants (Subscript[C, n]). This situation is similar with that of indefinite integrals.  The integration constants have to be defined from the boundary conditions (BC) at the left and right end of the interval, where the solution is considered. There should be exactly N boundary conditions to find N integration constants. 

If the independent variable is time, BC at the left end are called initial conditions. If the re are only initial conditions, the problem is called evolutional and in general is solved simpler than a problem with BC at both ends.

As an example of an evolutional problem consider a particle moving in one dimension under the influence of a constant force

m (x[t])''==F

where x[t] is particle's coordinate and m is its mass. First, DE is reduced to the system of first-order DEs

(v[t])'==(F/m)
(x[t])'==v[t].

for the velocity and coordinate. Then the first equation is integrated directly, left and right parts

∫((dv[t])/(dt))dt=∫dv[t]=v[t]=∫(F/m)dt=(F/m)t+(Subscript[C, 1])

Then the resulting equation for x[t] is also integrated directly

∫((dx[t])/(dt))dt=∫dx[t]=x[t]=∫((F/m)t+(Subscript[C, 1]))dt=(F/m)((t^2)/2)+(Subscript[C, 1])t+(Subscript[C, 2])

One can check that the solution

x[t]=(F/m)((t^2)/2)+(Subscript[C, 1])t+(Subscript[C, 2])

satisfies the initial second-order DE. In fact, solution of this structure could be guessed. Now the boundary conditions are found for the initial conditions

v[0]==(Subscript[v, 0]),     x[0]==(Subscript[x, 0]).

This yields

(Subscript[C, 1])==(Subscript[v, 0]),          (Subscript[C, 2])==(Subscript[x, 0])

and the well-known formulas for the motion with a constant acceleration a=F/m

v[t]=at+(Subscript[v, 0])
x[t]=((a (t^2))/2)+(Subscript[v, 0])t+(Subscript[x, 0]).

Solution of differential equations with Mathematica

The Mathematica command for analytical solution of DEs is DSolve that takes both DEs and systems of DEs as its argument. In the latter case, individual equations of a system have to be organized as a List. For instance, equation of motion in one dimension for a body under the influence of a constant force is solved as

DSolve[m x''[t]==F,x[t],t]

It can be used to define the solution of the DE as a function of time

xt[t_]=x[t]/.DSolve[m x''[t]==F,x,t][[1]]

Let us solve this equation again formulating it as a system of two first-order DEs

sol=DSolve[{m v'[t]==F,x'[t]==v[t]},{x[t],v[t]},t]

Differential equations can be solved together with boundary conditions:

sol=DSolve[{m v'[t]==F,x'[t]==v[t],x[0]==x0,v[0]==v0},{x[t],v[t]},t] //Expand

Again, the time-dependent solution is defined as

xt[t_]=x[t]/.sol[[1]]
vt[t_]=v[t]/.sol[[1]]

Let us now solve the equation of projectile motion in the x-y plane in the vector form. To make the vector DE solvable by Mathematica, we first have to convert the equation for the vector (list) to a list of equations with Thread.

FVec={-m g,0};       RVec[t_]={z[t],x[t]};
DSolve[Thread[m RVec''[t]==FVec],RVec[t],t]

Numerical solution of DEs in Mathematica is done by NDSolve that we will consider below

Analytical solution of differential equations

Only DEs of special kinds can be solved analytically. No general analytical method of solution of DEs exists. Below we will consider some cases in which DEs can be solved analytically

Separation of variables

If a first-order DE has a form

(x[t])'==g[x]f[t]

one can integrate this DE directly after division by g[x]

∫(1/(g[x]))((dx[t])/(dt))dt=∫((dx)/(g[x]))==∫f[t]dt.

One can see that on the left there is only x and on the right there is only t. One can see it even clearer if one formally rewrites the initial equation in the form

((dx)/(g[x]))=f[t]dt

before the integration. This is why this situation is called "separation of variables".

Viscous friction

As an example of separation of variables consider a particle moving in one dimension under the influence of a viscous friction

(Subscript[F, fr])=-(Subscript[k, 1]) v

where (Subscript[k, 1]) is a constant. For a sphere of radius R moving slowly enough in a liquid of viscosity η (so that the motion of the liquid is laminar) the friction coefficient has the form

(Subscript[k, 1])=(Piecewise[{{6π η R,, solid sphere}, {4π η R,, air bubble.}}])

The first formula is also applicable to a solid sphere moving in the air slowly enough. If the shape of the body is non-spherical, then the numerical factor is modified. (But for a body of an arbitrary shape there is also a viscous torque that causes body's rotation. We will neglect this effect.)  The Newton's second law for the velocity

(v[t])'==(F/m)==-α v,        α≡((Subscript[k, 1])/m)

separates as

((dv)/v)==-αdt.

Integration yields

∫((dv)/v)=Log[v]==- α t+C

and finally, with the help of the initial condition v[0]==(Subscript[v, 0]),

v=(e^(- α t))(e^C)=(Subscript[v, 0])(e^(- α t)).

Thus the velocity of a body in a viscous liquid exponentially decreases with time. Now the DE for the coordinate

(x[t])'==v[t]=(Subscript[v, 0])(e^(- α t))

can be immediately integrated resulting in

x[t]=∫(Subscript[v, 0])(e^(- α t))dt=-((Subscript[v, 0])/α)(e^(- α t))+C

From the initial condition x[0]==(Subscript[x, 0]) one can eliminate C and obtain

x[t]=(Subscript[x, 0])+((Subscript[v, 0])/α)(1-(e^(- α t)))

One can see that the particle moves only a finite distance (Subscript[v, 0])/α in the liquid.

Media resistance at high speeds

Resistance (friction) acting on a body fast moving in a liquid or in the air does not depend on the viscosity and is turbulent. The corresponding friction force is impossible to obtain exactly but qualitative arguments and dimensional estimations yield

(Subscript[F, fr])=-(Subscript[k, 2]) (v^2),

where

(Subscript[k, 2])=(1/2)ρSγ

ρ being the density of the liquid, S being the cross-section of the body, and γ being a numerical factor. This friction force can be understood as the Bernoulli pressure P=(1/2)((ρv)^2) multiplied by the cross-section S, plus a numerical factor γ that is difficult to find: (Subscript[F, fr])=γPS.

In this case the DE for the velocity is nonlinear

(v[t])'==-(Subscript[α, 2]) (v^2),        (Subscript[α, 2])≡((Subscript[k, 2])/m)

but still it can be solved by the separation of variables

((dv)/(v^2))==-(Subscript[α, 2])dt

and then

∫((dv)/(v^2))=((-1)/v)==- (Subscript[α, 2]) t+C

that is, with the help of the boundary condition,

v[t]=(1/((Subscript[α, 2]) t+1/(Subscript[v, 0])))

It is noteworthy that at long times the velocity (Subscript[v,  ])does not depend on the initial velocity (Subscript[v, 0]), so that the system forgets the initial condition. This solution, as well as the expression for the friction force, is valid for v>0.

Time dependence of the coordinate x[t] is obtained by the integration of the velocity

x[t]=∫v[t]dt=∫(1/((Subscript[α, 2]) t+1/(Subscript[v, 0])))dt=(1/(Subscript[α, 2]))Log[(Subscript[α, 2]) t+(1/(Subscript[v, 0]))]+C

From the initial condition x[0]==(Subscript[x, 0]) one can eliminate C and obtain

x[t]=(Subscript[x, 0])+(1/(Subscript[α, 2]))Log[(((Subscript[α, 2]) t+(1/(Subscript[v, 0])))/(1/(Subscript[v, 0])))]=(Subscript[x, 0])+(1/(Subscript[α, 2]))Log[1+(Subscript[v, 0])(Subscript[α, 2]) t].

Here x[t] increases infinitely but logarithmically slow at long times. However, at long times the speed becomes small and one has to switch to the viscous friction force above. The crossover between the two different expressions for the friction force should be at the speed at which both friction forces are comparable.

Using an integration factor

In some cases one can find a function to multiply by the differential equation, after which the left part containing the unknown function becomes a derivative. In this case one can integrate the equation one time and solve or simplify the problem. Finding an integration factor, if it is possible, is a matter of luck and skill.

Energy conservation

The most prominent example of using an integrating factor in physics is probably the derivation of the energy conservation law from the equation of motion of conservative systems. Consider a point body under the action of a conservative (potential) force

m r''==F[r]=-((∂U[r])/(∂r)),

Dot-multiplying this DE by the integrating factor r' and using the chain rule for differentiation, one can bring this equation into the form of time derivatives on the right and on the left

m r''·r'=m(1/2)((((r')^2)))'==-((∂U[r])/(∂r))·r'=-U'[r]

Integrating this equation over time and substituting r'=v, one obtains the energy conservation law

((m (v^2))/2)+U[r]==E,

where the conserved total energy E arises as an integration constant. For a motion in one dimension the energy conservation law can be used to completely solve the problem. Solving the equation

((m ((x')^2))/2)+U[x]==E,

for x', one obtains the first-order DE

x'==±(Sqrt[((2(E-U[x]))/m)]).

In this DE the variables separate and it can be integrated as

∫(1/(Sqrt[E-U[x]]))dx ==±(Sqrt[(2/m)])t,

so that the problem reduces to the evaluation of the integral on the left. After that one has to find x[t] from a nonlinear equation.

In multi-dimensional cases, the energy conservation law alone is insufficient to define the motion because it is one equation containing more then one unknown function: x[t], y[t], etc., an underdetermined problem.

As an illustration consider a pendulum of mass m and length L described by the angle ϕ (measured from the vertical down direction) as the dynamical variable and having the energy

E=((m (v^2))/2)+mgh=(((mL^2)(ϕ^(′2)))/2)+mgL(1-Cos[ϕ]).

The equation of motion of the pendulum is

ϕ''+(Subscript[ω, 0]^2)Sin[ϕ]==0,       (Subscript[ω, 0])≡(Sqrt[(g/L)]).

Instead of solving this equation, one can use its one-time-integrated form following from the energy conservation law E=const:

ϕ'==±(Sqrt[((2(E-U[ϕ]))/(mL^2))])=±(Sqrt[((2(E-mgL(1-Cos[ϕ])))/(mL^2))]).

This equation can be integrated by separation of variables

∫(1/(Sqrt[E-mgL(1-Cos[ϕ])]))dϕ ==±(Sqrt[(2/(mL^2))])t.

The pendulum has turning points where its kinetic energy turns to zero and the sign in front of the square root changes sign. Let us call the corresponding angle (Subscript[ϕ, max]). It depends on E and can be found from the equation

E==mgL(1-Cos[(Subscript[ϕ, max])].

In the sequel we will use (Subscript[ϕ, max]) as the measure of pendulum's energy. A quarter of the pendulum's period T corresponds to its motion from ϕ=0 to ϕ=(Subscript[ϕ, max]). Thus the period can be found from

(∫_0^(Subscript[ϕ, max]))(1/(Sqrt[mgL(Cos[ϕ]-Cos[(Subscript[ϕ, max])])]))dϕ=(Sqrt[(2/(mL^2))])(T/4),

that is,

T=(2^(3/2))(Sqrt[(L/g)])(∫_0^(Subscript[ϕ, max]))(1/(Sqrt[Cos[ϕ]-Cos[(Subscript[ϕ, max])]]))dϕ.

that is,

T[ϕmax_]=(2^(3/2))(Sqrt[(L/g)])Integrate[(1/(Sqrt[Cos[ϕ]-Cos[ϕmax]])),{ϕ,0,ϕmax},Assumptions->Cos[ϕmax]<1]

This is a pretty complicated expression and Mathematica cannot simplify it. The plot

Plot[(1/(2π))T[ϕmax]/.{L->1,g->1},{ϕmax,0,π}]

shows that for (Subscript[ϕ, max])≪1 the period approaches the known result

T=2π(Sqrt[(L/g)]).

Indeed,

Limit[T[ϕmax],{ϕmax->0}][[1]]

However, Mathematica produces a wrong series expansion of the result at (Subscript[ϕ, max])≪1:

Series[T[ϕmax],{ϕmax,0,3}]

There are also glitches in plotting. All this underscores the importance of analytical work. One can put the result in a much better shape rewriting it as

T=2(Sqrt[(L/g)])(∫_0^(Subscript[ϕ, max]))(1/(Sqrt[((Sin[((Subscript[ϕ, max])/2)])^2)-((Sin[(ϕ/2)])^2)]))dϕ

using the new integration variable ξ defined by

Sin[ξ]=((Sin[(ϕ/2)])/(Sin[((Subscript[ϕ, max])/2)])).

The integral is transformed as follows

T=2(Sqrt[(L/g)])(∫_0^(π/2))(1/(Sqrt[((Sin[((Subscript[ϕ, max])/2)])^2)-((Sin[((Subscript[ϕ, max])/2)])^2)((Sin[ξ])^2)]))((dϕ)/(dξ))dξ=2(Sqrt[(L/g)])(∫_0^(π/2))(1/(Sin[((Subscript[ϕ, max])/2)]Cos[ξ]))((2Cos[ξ]Sin[((Subscript[ϕ, max])/2)])/(Cos[(ϕ/2)]))dξ

Now with the help of

Cos[(ϕ/2)]=(Sqrt[1-((Sin[(ϕ/2)])^2)])=(Sqrt[1-((Sin[((Subscript[ϕ, max])/2)])^2)((Sin[ξ])^2)])

one finally obtains

T=4(Sqrt[(L/g)])(∫_0^(π/2))((dξ)/(Sqrt[1-((Sin[((Subscript[ϕ, max])/2)])^2)((Sin[ξ])^2)]))=4(Sqrt[(L/g)])K[((Sin[((Subscript[ϕ, max])/2)])^2)],

where

K[m]=(∫_0^(π/2))((dξ)/(Sqrt[1-m ((Sin[ξ])^2)]))

is the elliptic intergal of the first kind. In this more transparent representation, one can easily see what is the period at small amplitudes by expanding the integrand in small (Subscript[ϕ, max]). Using (((1 - x))^(-1/2))≃1+x/2  one easily obtains

T≃4(Sqrt[(L/g)])(∫_0^(π/2))(1+(1/2)((Sin[((Subscript[ϕ, max])/2)])^2)((Sin[ξ])^2))dξ=2π(Sqrt[(L/g)])(1+(1/4)((Sin[((Subscript[ϕ, max])/2)])^2))≃2π(Sqrt[(L/g)])(1+(1/16)(Subscript[ϕ, max]^2)).

Let us plot this result

T[ϕmax_]=4(Sqrt[(L/g)])EllipticK[((Sin[((ϕmax)/2)])^2)]
Plot[(1/(2π))T[ϕmax]/.{L->1,g->1},{ϕmax,0,π}]

the same as above but a more accurate plot.

Mathematica yields

T[ϕmax_]=4(Sqrt[(L/g)])(∫_0^(π/2))((\[DifferentialD]ξ)/(Sqrt[1-((Sin[((ϕmax)/2)])^2)((Sin[ξ])^2)]))
Plot[(1/(2π))T[ϕmax]/.{L->1,g->1},{ϕmax,0,π}]

again a different but equivalent analytical expression. The series expansion of the result has the form

Series[T[(Subscript[ϕ, max])],{(Subscript[ϕ, max]),0,3}]

now a correct expression.

Viscous friction with a constant force

As another example let us consider a particle moving in one dimension under the influence of the viscous friction and a time-dependent force F[t]. The equation for the velocity has the form

(v[t])'==-α v[t]+f[t],      f[t]≡((F[t])/m)

with the positive direction of the x-axis along the force F. If the force is time-independent, the variables separate,

((dv)/(-α v+f))==dt

and the DE can be solved in this way. However, for the time-dependent force separation of variables does not occur. In tis case, fortunately, one can find an integration factor that solves the problem. It is convenient to move the v-dependent term to the left,

(v[t])'+α v[t]==f.

Now, multiplication of this DE by the integration factor (e^(α t)) makes the left part a derivative:

(e^(α t))((v[t])'+α v[t])=(∂_t)((e^(α t))v[t])==(e^(α t))f[t],

so that the equation can be integrated as

(e^(α t))v[t]=∫(e^(α t))f[t]dt.

In particular, for a time-independent force such as gravity the integral can be immediately done

(e^(α t))v[t]=(1/α)(e^(α t))f+C.

The solution is

v[t]=(f/α)+(e^(-α t))C=(Subscript[v, 0])(e^(-α t))+(f/α)(1-(e^(-α t))),

where the integration constant C was obtained from the boundary condition. Asymptotically at large times, the particle forgets its initial velocity and takes the stationary velocity value (Subscript[v, ∞])=(f/α). This result can be obtained directly from the original DE setting (v[t])'⇒ 0. After v[t] has been found, one can straighforwardly find v[t] by integration.

Non-evolutionary boundary conditions: Shooting

If boundary conditions are set at both sides of the interval of the independent variable, there is a little difference with the case of BCs on one side if there is a general analytical solution of the DE. As before, integration constants are determined from the BC after the general solution is obtained. If, however the DE can only be solved numerically, there is a big difference. DEs of the evolution type (all BCs on the left) are solved moving from the left to right along the trajectory that is uniquely defined by the BCs on the left. If BCs are set at both ends of the interval, one cannot immediately find the trajectory because BCs on the left do not contain the whole information about it. There is a freedom in the trajectory that is used to satisfy the BCs on the right end. Numerically one can calculate trial trajectories and then choose those that satisfy the BCs on the right. Of course, this procedure is more complicated and time consuming than solving evolution equations. Because the first idea of such type of a problem is hitting a target at the right end of the interval, this class of problems and the method of their solution outlined above is called shooting.

BCs on both ends: Free straight motion

As an example, suppose a body moves straight with no forces applied and was t=(Subscript[t, 0]) at x=(Subscript[x, 0]). At t=(Subscript[t, 1])>(Subscript[t, 0]) the body should be at x=(Subscript[x, 1]). Find the motion of the body. The solution of the differential equation

x''[t]==0

is

x[t]=(Subscript[C, 1])t+(Subscript[C, 2])

and the boundary conditions are

x[(Subscript[t, 0])]==(Subscript[x, 0]),          x[(Subscript[t, 1])]==(Subscript[x, 1])

This yields the integration constants

(Subscript[C, 1])=(((Subscript[x, 1])-(Subscript[x, 0]))/((Subscript[t, 1])-(Subscript[t, 0]))),         
(Subscript[C, 2])=(Subscript[x, 0])-(Subscript[C, 1])(Subscript[t, 0])=(((Subscript[x, 0])((Subscript[t, 1])-(Subscript[t, 0]))-((Subscript[x, 1])-(Subscript[x, 0]))(Subscript[t, 0]))/((Subscript[t, 1])-(Subscript[t, 0])))=(((Subscript[x, 0])(Subscript[t, 1])-(Subscript[x, 1])(Subscript[t, 0]))/((Subscript[t, 1])-(Subscript[t, 0]))),

Obviously (Subscript[C, 1])=v, the constant velocity. This velocity is not given as the inital condition but follows from both boundary conditions.

More general BC: Vertical motion with gravity

In many cases the right boundary of the interval is not set but there is a condition that plays a role of a BC. As an example, a stone it thrown vertically up and it has to hit a target at height h with a speed (Subscript[v, f]) (final velocity). The general solution of the equation of motion

z''[t]==-g

has the form

z[t]=-(1/2)(gt^2)+(Subscript[C, 1])t+(Subscript[C, 2])
v[t]=-gt+(Subscript[C, 1])

and the BCs read

z[0]==0,         z[(Subscript[t, f])]==h,          v[(Subscript[t, f])]=(Subscript[v, f]).

There are three boundary conditions for a second-order DE but this is not a contradiction because (Subscript[t, 1]) is an additional unknown. One obtains

(Subscript[C, 2])==0

and the system of equations

-(1/2)(Subscript[gt, f]^2)+(Subscript[C, 1])(Subscript[t, f])==h,          -(Subscript[gt, f])+(Subscript[C, 1])==(Subscript[v, f]).

Mathematica yields

Solve[{-(1/2)g (Subscript[t, f]^2)+(Subscript[C, 1])(Subscript[t, f])==h,-g (Subscript[t, f])+(Subscript[C, 1])==(Subscript[v, f])},{(Subscript[t, f]),(Subscript[C, 1])}]

Here the solution with (Subscript[t, f])>0 is physically relevant and one obtains, finally

z[t]=-(1/2)(gt^2)+(Sqrt[2 g h+(Subscript[v, f]^2)])t
v[t]=-gt+(Sqrt[2 g h+(Subscript[v, f]^2)]).

The main result of this calculation is that the initial shooting velocity has to be (Subscript[v, 0])=(Sqrt[(Subscript[v, f]^2)+2 g h])>(Subscript[v, f]).

Projectile shooting with a viscous friction

Let us consider shooting of a projectile that has to hit a target at a horizontal distance L taking into account the air resistance described (for simplicity) as a viscous friction. The initial speed of the projectile is (Subscript[v, 0]). The equations of motion are linear and thus they split into a system of uncoupled equations for the horizontal and vertical motion

(Subscript[v, x])'==-α (Subscript[v, x]),               x'==(Subscript[v, x])
(Subscript[v, z])'==-α (Subscript[v, z])-g,         z'==(Subscript[v, z])

Initial conditions are

x[0]==0,                         z[0]==0
(Subscript[v, x])[0]== (Subscript[v, 0])Cos[θ],      (Subscript[v, z])[0]== (Subscript[v, 0])Sin[θ],

where θ is the shooting angle. The final conditions are

x[(Subscript[t, f])]==L,           z[(Subscript[t, f])]==0

Solution of the system of DEs above contains four integration constants, and also there are unknown (Subscript[t, f]) and θ. For this we have six boundary conditions, so that the problem is well determined. Solving equations for the velocity components as above, one obtains

(Subscript[v, x])[t]=(Subscript[v, 0])Cos[θ](e^(-α t))
(Subscript[v, z])[t]=(Subscript[v, 0])Sin[θ](e^(-α t))-(g/α)(1-(e^(-α t)))

Integrating these equations yield

x[t]=∫(Subscript[v, 0])Cos[θ](e^(-α t))dt=-(Subscript[v, 0])Cos[θ]((e^(-α t))/α)+(Subscript[C, x])
z[t]=∫((Subscript[v, 0])Sin[θ](e^(-α t))-(g/α)(1-(e^(-α t))))dt=-(Subscript[v, 0])Sin[θ]((e^(-α t))/α)-(g/α)(t+((e^(-α t))/α))+(Subscript[C, z])

Finding the integration constants from the initial conditions one obtains

x[t]=(Subscript[v, 0])Cos[θ]((1-(e^(-α t)))/α)
z[t]=(Subscript[v, 0])Sin[θ]((1-(e^(-α t)))/α)-(g/α)(t-((1-(e^(-α t)))/α))=((Subscript[v, 0])Sin[θ]+(g/α))((1-(e^(-α t)))/α)-(g/α)t

It remains now to find (Subscript[t, f]) and θ from the final conditions above, that is

(Subscript[v, 0])Cos[θ]((1-(e^(-α (Subscript[t, f]))))/α)==L
((Subscript[v, 0])Sin[θ]+(g/α))((1-(e^(-α (Subscript[t, f]))))/α)-(g/α)(Subscript[t, f])==0.

Eliminating the (Subscript[t, f]) combination with the help of the first equation, one obtains

((α/g)(Subscript[v, 0])Sin[θ]+1)(L/((Subscript[v, 0])Cos[θ]))==(Subscript[t, f]).

Substituting this result into the first equation, one obtains the transcedental equation for the shooting angle

(Subscript[v, 0])Cos[θ](1/α)(1-Exp[-α((α/g)(Subscript[v, 0])Sin[θ]+1)(L/((Subscript[v, 0])Cos[θ]))])-L==0

In reality, the air resistance is small, thus this equation can be simplified by expansion at small α

Normal[Series[(Subscript[v, 0])Cos[θ](1/α)(1-Exp[-α((α/g)(Subscript[v, 0])Sin[θ]+1)(L/((Subscript[v, 0])Cos[θ]))])-L,{α,0,2}]]

Simplifying this, one obtains

Simplify[%((2 g  Cos[θ])/(α (Subscript[v, 0])L))]

thus for small α the equation has the form

Sin[2θ]-((g L)/(Subscript[v, 0]^2))-((αL)/((Subscript[v, 0])Cos[θ]))(Sin[2θ]-((g L)/(3(Subscript[v, 0]^2))))==0

Setting α=0 one recovers the well known result

Sin[2θ]=((g L)/(Subscript[v, 0]^2))

that yields two values of θ corresponding to the high and low trajectories. The target can be hit only if (Subscript[v, 0]^2)≥gL. The final form of the solution is

θ=(Subscript[ϑ, 1])=(1/2)ArcSin[((g L)/(Subscript[v, 0]^2))]<(π/4),         θ=(Subscript[ϑ, 2])=(1/2)(π-ArcSin[((g L)/(Subscript[v, 0]^2))])>(π/4)

Let us plot the trajectories with and without the air resistance

Parameters={g->1.,L->1,(Subscript[v, 0])->1.1,α->0.05};
(Subscript[ϑ, 1])=(1/2)ArcSin[((g L)/(Subscript[v, 0]^2))]/.Parameters ;               (Subscript[ϑ, 1]) (( 180)/π)
(Subscript[ϑ, 2])=(1/2)(π-ArcSin[((g L)/(Subscript[v, 0]^2))])/.Parameters ;     (Subscript[ϑ, 2])(( 180)/π)
x0[t_,θ_]:=(Subscript[v, 0])Cos[θ]t;                   z0[t_,θ_]:=-(1/2)g (t^2)+(Subscript[v, 0])Sin[θ]t;      (* Without air resistance *)
x[t_,θ_]:=(Subscript[v, 0])Cos[θ]((1-(E^(-α t)))/α);         z[t_,θ_]:=((Subscript[v, 0])Sin[θ]+(g/α))((1-(E^(-α t)))/α)-(g/α)t; (* With air resistance *)

Show[
ParametricPlot[{x0[t,(Subscript[ϑ, 2])],z0[t,(Subscript[ϑ, 2])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 2])]))/.Parameters},PlotStyle->{Red,Dashed}],
ParametricPlot[{x0[t,(Subscript[ϑ, 1])],z0[t,(Subscript[ϑ, 1])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 1])]))/.Parameters},PlotStyle->{Blue,Dashed}],
ParametricPlot[{x[t,(Subscript[ϑ, 2])],z[t,(Subscript[ϑ, 2])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 2])]))/.Parameters},PlotStyle->{Red}],
ParametricPlot[{x[t,(Subscript[ϑ, 1])],z[t,(Subscript[ϑ, 1])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 1])]))/.Parameters},PlotStyle->{Blue}]
]

One can see that the high trajectory is more sensitive to the air resistance because the time and exposure to the resistance is longer. In fact, the model of the viscous friction is not good for real shooting because in assumes small speeds whereas the actual speeds are high. A better mathematical model includes the (v^2) damping but in this case the equations of motion are nonlinear and cannot be solved analytically.

Let us now obtain the correction to the shooting angle due to the air resistance. Search for the solution in the form with corrections linear in small α:

θ=(Subscript[θ, i])=(Subscript[ϑ, i])+(Subscript[δθ, i]),       i=1,2,      |(Subscript[δθ, i])|≪1

Inserting this into the equation, using

Sin[2θ]=Sin[2((Subscript[ϑ, i])+(Subscript[δθ, i]))]≃Sin[2(Subscript[ϑ, i])]+2Cos[2(Subscript[ϑ, i])](Subscript[δθ, i])

and keeping only linear-α terms in the equation, one obtains

2Cos[2(Subscript[ϑ, i])](Subscript[δθ, i])-((αL)/((Subscript[v, 0])Cos[(Subscript[ϑ, i])]))(Sin[2(Subscript[ϑ, i])]-((g L)/(3(Subscript[v, 0]^2))))==0

or,

(Subscript[δθ, i])=((αL)/((Subscript[v, 0])Cos[2(Subscript[ϑ, i])]Cos[(Subscript[ϑ, i])]))((g L)/(3(Subscript[v, 0]^2))).

This formula gives the corrections to the no-resustance formula for the air resistance. For the low trajectory  Cos[2(Subscript[ϑ, 1])]>0 and one has to aim higher that is intuitively clear. For the high trajectory Cos[2(Subscript[ϑ, 2])]<0 and one has to aim lower. For the maximal-distance shot, (Subscript[ϑ, 1])=(Subscript[ϑ, 2])=π/4, one has Cos[2(Subscript[ϑ, 2])]=0 and the formula above fails. In this case one cannot hit the target at all because of the air resistance.

Let us plot the trajectories with and without the air resistance correction

Parameters={g->1.,L->1,(Subscript[v, 0])->1.1,α->0.05};
(Subscript[ϑ, 1])=(1/2)ArcSin[((g L)/(Subscript[v, 0]^2))]/.Parameters ;                                                  (* uncorrected *)     
              (Subscript[θ, 1])=(Subscript[ϑ, 1])+((α L)/((Subscript[v, 0])Cos[2(Subscript[ϑ, 1])]Cos[(Subscript[ϑ, 1])]))((g L)/(3(Subscript[v, 0]^2)))/.Parameters ;    (* corrected *)
(Subscript[ϑ, 2])=(1/2)(π-ArcSin[((g L)/(Subscript[v, 0]^2))])/.Parameters ;                                          (* uncorrected *)
              (Subscript[θ, 2])=(Subscript[ϑ, 2])+((α L)/((Subscript[v, 0])Cos[2(Subscript[ϑ, 2])]Cos[(Subscript[ϑ, 2])]))((g L)/(3(Subscript[v, 0]^2)))/.Parameters ;       (* corrected *)
x[t_,θ_]:=(Subscript[v, 0])Cos[θ]((1-(E^(-α t)))/α);            z[t_,θ_]:=((Subscript[v, 0])Sin[θ]+(g/α))((1-(E^(-α t)))/α)-(g/α)t;    (* trajectory *)

Show[
ParametricPlot[{x[t,(Subscript[ϑ, 2])],z[t,(Subscript[ϑ, 2])]}/.Parameters,{t,0,((1.2L)/((Subscript[v, 0])Cos[(Subscript[ϑ, 2])]))/.Parameters},PlotStyle->{Red}], (* uncorrected *)
ParametricPlot[{x[t,(Subscript[ϑ, 1])],z[t,(Subscript[ϑ, 1])]}/.Parameters,{t,0,((1.2L)/((Subscript[v, 0])Cos[(Subscript[ϑ, 1])]))/.Parameters},PlotStyle->{Blue}], (* uncorrected *)
ParametricPlot[{x[t,(Subscript[θ, 2])],z[t,(Subscript[θ, 2])]}/.Parameters,{t,0,((1.1L)/((Subscript[v, 0])Cos[(Subscript[θ, 2])]))/.Parameters},PlotStyle->{Thick,Red}], (* corrected *)
ParametricPlot[{x[t,(Subscript[θ, 1])],z[t,(Subscript[θ, 1])]}/.Parameters,{t,0,((1.1L)/((Subscript[v, 0])Cos[(Subscript[θ, 1])]))/.Parameters},PlotStyle->{Thick,Blue}] (* corrected *)
]

The target is almost hit. Note that the correcion due to the air resistance has been calculated approximately!

Linear DEs with constant coefficients

Linear DEs (and systems of DEs) with constant coefficients is probably the most important class of DEs. The solution of homogeneous equations can be searched in the form

x[t] =((Ce)^(λt))

that reduces the DE to an algebraic equation.

Viscous friction

For instance, the solution of the equation of motion of a body with a viscous friction

v'[t]+α v[t]==0

can be searched in the form  v[t]=((Ce)^(λt)). Substituting this into the equation, one obtains the equation for λ

((Cλe)^(λt))+((αCe)^(λt))=((Ce)^(λt))(λ+α)==0 .

Since we are looking for a nontrivial solution with C≠0, this equation yields λ=-α. Of course, this simple equation can be solved by other methods as well.

Harmonic oscillator

Let us consider the equation of one-dimensional motion of a body on a spring providing the restoring force F=-k x

mx''[t]==F[x]=-kx[t]

that can be rewritten in the generic form of the oscillator equation

x''[t]+(Subscript[ω, 0]^2)x[t]==0  ,           (Subscript[ω, 0])≡(Sqrt[(k/m)]).

The solution can be searched in the form x[t] =((Ce)^(λt)) that yields the quadratic equation for λ

(λ^2)+(Subscript[ω, 0]^2)==0 .

The solution of this equation is imaginary, λ=±i (Subscript[ω, 0]), and the general solution of the DE has the form

x[t] =(SubPlus[C])(e^((Subscript[iω, 0])t))+(SubMinus[C])(e^(-(Subscript[iω, 0])t)),

where (Subscript[C, ±]) are integration constants. Since the solution x[t] must be real, one should express complex exponentials via Sin and Cos

(e^(±(Subscript[iω, 0])t))=Cos[(Subscript[ω, 0])t]±i Sin[(Subscript[ω, 0])t].

In this way one obtains the solution as a linear combination  of sinusoidal functions

x[t] =(Subscript[C, 1])Cos[(Subscript[ω, 0])t]+(Subscript[C, 2])Sin[(Subscript[ω, 0])t],

where (Subscript[C, 1,2]) are other integration constants. One can check that this solution satisfies the initial DE. One could avoid using complex exponentials and work with sinusoidal functions from the very beginning. However, presence of two different functions, Sin and Cos, complicates calculations.

The velocity of the body is given by

v[t] =x'[t] =-(Subscript[C, 1])(Subscript[ω, 0])Sin[(Subscript[ω, 0])t]+(Subscript[C, 2])(Subscript[ω, 0])Cos[(Subscript[ω, 0])t].

With the help of the initial conditions

v[0]==(Subscript[v, 0]),     x[0]==(Subscript[x, 0])

one can find the integration constants

(Subscript[C, 1])=(Subscript[x, 0]),     (Subscript[C, 2])=((Subscript[v, 0])/(Subscript[ω, 0]))

and the final form of the solution

x[t] =(Subscript[x, 0])Cos[(Subscript[ω, 0])t]+((Subscript[v, 0])/(Subscript[ω, 0]))Sin[(Subscript[ω, 0])t].

Systems of linear DEs with constant coefficients

Similar method works for systems of linear equations with constant coefficients, however, with some linear algebra.

Damped harmonic oscillator

Let us consider again a harmonic oscillator but this time we represent the problem as a system of two first-order DEs. To avoid repetition, we add viscous friction that is responsible for damping of the oscilations. The system of equation for a damped oscillator has the form

x''[t]+αx'[t]+(Subscript[ω, 0]^2)x[t]==0  ,           α≡((Subscript[k, 1])/m),           (Subscript[ω, 0])≡(Sqrt[(k/m)]).

or

v'[t]+αv[t]+(Subscript[ω, 0]^2)x[t]==0  
x'[t]-v[t]==0.

We search for a solution in the form

v[t]=(Subscript[C, v])(e^(-λt)),       x[t]=(Subscript[C, x])(e^(-λt)).

Substituting this into the DEs, calculating derivatives, and cancelling (e^(λt)), one obtains the system of linear equations for the coefficients (Subscript[C, v]) and (Subscript[C, x])

(α-λ)(Subscript[C, v])+(Subscript[ω, 0]^2)(Subscript[C, x])==0
-(Subscript[C, v])-(Subscript[λC, x])==0.

Since the inhomogeneous (right) part of this system of equations is zero, the solution will be (Subscript[C, x])=(Subscript[C, y])=0 if the two equations are independent. Only if the two equations are equivalent (and thus one of them can be dropped) there is a solution with nonzero (Subscript[C, v]) and (Subscript[C, x]) that we are looking for. Mathematically dependence of linear equations in a system of linear equations is guaranteed by the zero determinant of this system, in our case

Det[(({{α-λ, (Subscript[ω, 0]^2)}, {-1, -λ}}))=-λ(α-λ)+(Subscript[ω, 0]^2)==0.

The values of λ that satisfy this condition are called "Eigenvalues" and the problem to find them is called "eigenvalue problem". The quadratic equation for λ generalizes those above for the undamped harmonic oscillator and the body in a viscous liquid. The solution for λ is

(Subscript[λ, 1,2])=(α/2)±(Sqrt[((((α/2)))^2)-(Subscript[ω, 0]^2)]).

For α>(Sqrt[2])(Subscript[ω, 0]) both (Subscript[λ, 1]) and (Subscript[λ, 2]) are real and positive, and the solution has the general form

v[t]=(Subscript[C, v1])(e^(-(Subscript[λ, 1])t))+(Subscript[C, v2])(e^(-(Subscript[λ, 2])t)),          x[t]=(Subscript[C, x1])(e^(-(Subscript[λ, 1])t))+(Subscript[C, x2])(e^(-(Subscript[λ, 2])t))

Here, however, the constants are not independent. One can express (Subscript[C, v1]) via (Subscript[C, x1]) and (Subscript[C, v2]) via (Subscript[C, x2]) with the help of the system of linear equations for these coefficients. Even easier, one can obtain v[t] as the derivative of x[t]

x[t]=(Subscript[C, 1])(e^(-(Subscript[λ, 1])t))+(Subscript[C, 2])(e^(-(Subscript[λ, 2])t)),      v[t]=-(Subscript[λ, 1])(Subscript[C, 1])(e^(-(Subscript[λ, 1])t))-(Subscript[λ, 2])(Subscript[C, 2])(e^(-(Subscript[λ, 2])t)).

This is called "aperiodic regime" for a damped harmonic oscillator, there are no oscillations in this regime. For α<(Sqrt[2])(Subscript[ω, 0]) both (Subscript[λ, 1]) and (Subscript[λ, 2]) are complex with a positive real part and can be written as

(Subscript[λ, 1,2])=(α/2)±i(Sqrt[(Subscript[ω, 0]^2)-((((α/2)))^2)]).

Using

(e^(-(Subscript[λ, 1,2])t))=Exp[-(α/2)t](Cos[(Sqrt[(Subscript[ω, 0]^2)-((((α/2)))^2)])t]±i Sin[(Sqrt[(Subscript[ω, 0]^2)-((((α/2)))^2)])t]),

one can rewrite the solution in the form

x[t]=(Subscript[C, 1])Exp[-(α/2)t]Cos[(Sqrt[(Subscript[ω, 0]^2)-((((α/2)))^2)])t]+(Subscript[C, 2])Exp[-(α/2)t]Sin[(Sqrt[(Subscript[ω, 0]^2)-((((α/2)))^2)])t],

and the expression  for v[t] can be obtained by differentiation. Finally, (Subscript[C, 1]) and (Subscript[C, 2]) can be obtained from the initial conditions. In this periodic regime, the motion is oscillating but damped, and the frequency of the oscillations is reduced with respect to (Subscript[ω, 0]).

Systems of linear DEs with constant coefficients and exponential or sinusoidal inhomogeneous term: Resonance

The method of searching for the solution in the exponential form discussed above can be extended for the case when the inhomogeneous term is exponential or sinusoidal. As an example consider damped harmonic oscillator subject to a sinusoidal force,

x''[t]+αx'[t]+(Subscript[ω, 0]^2)x[t]==f[t]=(Subscript[f, 0])Cos[ωt]  ,                f[t]≡((F[t])/m)

Apart for the general solution of the corresponding homogeneous equation found above that tends to zero at large times, there is a particular solution of this system of equations that is the stationary response of the oscillator to the external force F[t]. This stationary response can be searched for in the form x[t] = A Cos[ω t+(Subscript[ϕ, 0])], where A is the amplitude and (Subscript[ϕ, 0]) is the phase shift.  However, this solution requires too many trigonometric operations and is cumbersome. A more elegant solution makes use of complex exponentials

(e^(iωt))=Cos[ωt]+i Sin[ωt]

and the linearity argument. By the latter, the particular solution of the equation

x''[t]+αx'[t]+(Subscript[ω, 0]^2)x[t]==f[t]=(Subscript[f, 0])(e^(iωt))  ,

is the sum of the solution of the initial DE and the additional DE

x''[t]+αx'[t]+(Subscript[ω, 0]^2)x[t]==f[t]=(Subscript[f, 0])i Sin[ωt] .

The solution of the former is real, whereas that of the latter is imaginary. Thus to obtain the physically meaningful part of the solution, it is sufficient to take the real part of the solution of the equation with (Subscript[f, 0])(e^(iωt)) on the right. This equation solves much better than the equation with (Subscript[f, 0])Cos[ωt]. Searching for the solution in the form

x[t]=A (e^(iωt)),

substituting this into the DE, using (i^2)=-1 and canceling (e^(iωt)), one obtains the equation for A

(-(ω^2)+iωα+(Subscript[ω, 0]^2))A==(Subscript[f, 0])

that yields

A==((Subscript[f, 0])/((Subscript[ω, 0]^2)-(ω^2)+iωα))   .

Now the physical part of the solution is given by

x[t]=Re[A (e^(iωt))].

Separating real and imaginary parts, one transforms this solution to

x[t]=Re[(Re[A]+iIm[A])( Cos[ωt]+i Sin[ωt])]=Re[Re[A]Cos[ωt]-Im[A]Sin[ωt]+i(...)]=Re[A]Cos[ωt]-Im[A]Sin[ωt].

One can see that there is a phase shift between the force and the response to it because of the second term above. One can rewrite this expression as

x[t]=(Sqrt[((Re[A])^2)+((Im[A])^2)])Cos[ω t+(Subscript[ϕ, 0])]=|A|Cos[ω t+(Subscript[ϕ, 0])],

where (Subscript[ϕ, 0]) is the phase shift and | A | is the absolute value of the complex number A,

|A|=(Sqrt[((Re[A])^2)+((Im[A])^2)])=(Sqrt[(SuperStar[AA])])=(Sqrt[((Subscript[f, 0])/((Subscript[ω, 0]^2)-(ω^2)+iωα)) ((Subscript[f, 0])/((Subscript[ω, 0]^2)-(ω^2)-iωα)) ])=((Subscript[f, 0])/(Sqrt[((((Subscript[ω, 0]^2)-(ω^2)))^2)+(ω^2)(α^2)])).

One can see that at small damping, α≪(Subscript[ω, 0]), the amplitude of the response of the oscillator to the sinusoidal force peaks at resonance, ω=(Subscript[ω, 0]).

Variation of constants for linear differential equations

If a general solution of a homogeneous linear DE or system of DEs is known, one can find the solution of the corresponding inhomogeneous problem with the method called "variation of constants".

Consider the first-order linear DE

x'[t]+P[t]x[t]+Q[t]==0

The homogeneous DE with Q[t]=0 separates

((dx)/x)==-P[t]dt

and solves to

x[t]=C Exp[-∫P[t]dt].

Now one can search the solution of the original inhomogeneous DE in the form

x[t]=C[t] Exp[-∫P[t]dt]

that gives the name to the method. Substituting this into the DE, after cancelletions one obtains

C'[t] Exp[-∫P[t]dt]+Q[t]==0,

a simple equation for C[t] that can be integrated to

C[t]=C-∫Q[t]Exp[∫P[t]dt]dt,

where C is another integration constant. Finally for x[t] one obtains

x[t]=C Exp[-∫P[t]dt]-Exp[-∫P[t]dt]∫Q[t]Exp[∫P[t]dt]dt,

a general solution of the homogeneous problem plus a particular solution of the inhomogeneous problem. Above we have solved the equation of motion of a particle with viscous friction with a constant force by finding an integration factor for the DE. This DE also could be solved by the method of variation of constants. For this particular problem the two methods are equivalent. However, the method of integration factor works for some nonlinear equations whereas the method of variation of constants can be generalized for high-order DEs and systems of DEs.

Matrix solution of systems of ODEs with constant coefficients

Every ODE or system of ODEs with constant coefficients can be written in the form of the first-order matrix ODE

X'[t]+A.X[t]+F[t]==0,

where A is a constant matrix and X[t] and F[t] are vectors (columns). Similarly to the solution in the preceding section, one can write the solution of this equation in the form

X[t]=Exp[-At].C-Exp[-At].∫Exp[At].F[t]dt

that contains matrix exponentials. Here C is an integration constant vector. The order of factors in this expression is important. This method does not work if A is a time-dependent matrix. Mathematica can calculate such expressions analytically for small matrices A and numerically for large matrices A.

Numerical solution of differential equations

As mentioned above, only a small fraction of differential equations can be integrated analytically. Already the first-order DE of the type

x'[t]==f[t,x[t]]

with a general function f does not have an analytical solution. Some DEs can be solved analytically but their solution is cumbersome. This is why it is important to develop numerical methods of solution of DEs. The basic method uses presentation of the solution as a power series around some point (Subscript[t, 0])

x[t]=a+b(t-(Subscript[t, 0]))+c(((t-(Subscript[t, 0])))^2)+...

The coefficients in this formula can be found from the series expansion

Series[x'[t]-f[t,x[t]],{t,t0,1}]

One obtains

a=x[(Subscript[t, 0])]
b=x'[t0]=f[t0,x[t0]]
c=(1/2)x''[t0]=(1/2) ((∂f)/(∂t))[t0,x[t0]]+(1/2)f[t0,x[t0]]((∂f)/(∂x))[t0,x[t0]]

etc. Keeping the powers up to the first one yields the simplest method of numerical solution of DEs, the Euler method. Starting at the initial point (Subscript[t, 0]) and using the boundary condition x[(Subscript[t, 0])]= a one makes a step h and calculates x[(Subscript[t, 1])] with (Subscript[t, 1])=(Subscript[t, 0])+h by the approximate formula

x[(Subscript[t, 1])]=x[(Subscript[t, 0])]+f[(Subscript[t, 0]),x[(Subscript[t, 0])]]h.

After calculating x[(Subscript[t, 1])] one renames (Subscript[t, 1])⇒(Subscript[t, 0]) and repeats the procedure many times to cover the whole time interval. Since the Euler method uses the lowest-order power expansion, it is not very accurate and one has to use a very small step h to obtain accurate results. Using the second-order power expansion, one obtains the elementary iteration in the form

x[(Subscript[t, 1])]=x[(Subscript[t, 0])]+f[(Subscript[t, 0]),x[(Subscript[t, 0])]]h+(1/2)(((∂f)/(∂t))[t0,x[t0]]+f[t0,x[t0]]((∂f)/(∂x))[t0,x[t0]])(h^2)

that is more accurate than the Euler method. There are very many different methods of numerical solution of DEs of different orders that are implemented in code libraries and, in particular, in Mathematica. Some of these methods use derivatives of f while others avoid using derivatives, taking the necessary information from previous iterations to represent derivatives in the form of differences. There are numerical methods for higher-order DEs and systems of DEs, both resolved and unresolved with respect to the highest derivative.

Let us illustrate the Euler and quadratic methods in comparison with the practically numerically exact solution by Mathematica's standard solver for a particular equation for which Mathematica cannot find an analytical solution.

(* Mathematica's default numerical method (high order, automatic precision) *)
f[t_,xx_]=Sin[xx]+Sin[t];
Equ=x'[t]==f[t,x[t]];
IniCond=x[0]==0;
tMax=20;
Sol=NDSolve[{Equ,IniCond},x,{t,0,tMax}]
xt[t_]:=x[t]/.Sol[[1]]
BestSolution=Plot[xt[t],{t,0,tMax},PlotRange->All]

One can see that the solution shows a stationary asymptotic behavior looking like a sinusoidal function oscillating around π. However, parametric plotting of {x[t], x'[t]} shows that this is not a sinusoidal function as the resulting shape is not exactly a circle

ParametricPlot[{xt[t]-π,f[t,xt[t]]},{t,0.6tMax,tMax},PlotRange->All]

DSolve[Equ,x[t],t]

(* Euler method *)
f[t_,xx_]=Sin[xx]+Sin[t];
t0=0;  tMax=20;  NSteps=50;       h=(tMax/NSteps);     
x0=0.;     
xtList={{t0,x0}};
Do[
t1=t0+h;       x1=x0+f[t0,x0]h;
xtList=Append[xtList,{t1,x1}];
t0=t1;             x0=x1
,{k,1,NSteps}];

EulerSolution=ListPlot[xtList,PlotStyle->Gray]

(* Quadratic method *)
f[t_,xx_]=Sin[xx]+Sin[t];
fDert[t_,xx_]=(∂_t)f[t,xx]; 
fDerx[t_,xx_]=(∂_xx)f[t,xx];
t0=0;  tMax=20;  NSteps=50;       h=(tMax/NSteps);     
x0=0.;     
xtList={{t0,x0}};
Do[
t1=t0+h;       x1=x0+f[t0,x0]h+(1/2)(fDert[t0,x0]+f[x0,t0]fDerx[t0,x0])(h^2);
xtList=Append[xtList,{t1,x1}];
t0=t1;             x0=x1
,{k,1,NSteps}];

QuadraticSolution=ListPlot[xtList,PlotStyle->Black]

Show[EulerSolution,QuadraticSolution,BestSolution]

Let us now try the two methods on a simplest DE of the type

x'[t]==(t^2),       x[0]=0

that has the solution

x[t]=((t^3)/3).

tMax=20;
Exact=Plot[((t^3)/3),{t,0,tMax}]

(* Euler method *)
f[t_,xx_]=(t^2);
t0=0;  tMax=20;  NSteps=10;       h=(tMax/NSteps);     
x0=0.;     
xtList={{t0,x0}};
Do[
t1=t0+h;       x1=x0+f[t0,x0]h;
xtList=Append[xtList,{t1,x1}];
t0=t1;             x0=x1
,{k,1,NSteps}];

EulerSolution=ListPlot[xtList,PlotStyle->Gray,PlotRange->All]

(* Quadratic method *)
f[t_,xx_]=(t^2);
fDert[t_,xx_]=(∂_t)f[t,xx]; 
fDerx[t_,xx_]=(∂_xx)f[t,xx];
t0=0;  tMax=20;  NSteps=10;       h=(tMax/NSteps);     
x0=0.;     
xtList={{t0,x0}};
Do[
t1=t0+h;       x1=x0+f[t0,x0]h+(1/2)(fDert[t0,x0]+f[x0,t0]fDerx[t0,x0])(h^2);
xtList=Append[xtList,{t1,x1}];
t0=t1;             x0=x1
,{k,1,NSteps}];

QuadraticSolution=ListPlot[xtList,PlotStyle->Black,PlotRange->All]

Show[Exact,EulerSolution,QuadraticSolution]

Note that for this equation the cubic numerical method would give an exact result. Since there is no stationary asymptotic behavior, numerical errors accumulate and numerical trajectories deviate more and more from the exact trajectory. Excercise: Plot relative deviations of the approximate solutions from the exact solution. Do they increase with time?

Numerically solvable differential equations with non-evolutionary boundary conditions: More shooting

Let us consider the shooting problem with more realistic quadratic friction

(Subscript[v, x])'==-(Subscript[α, 2]) (Subscript[v, x])(Sqrt[(Subscript[v, x]^2)+(Subscript[v, z]^2)]),               x'==(Subscript[v, x])
(Subscript[v, z])'==-(Subscript[α, 2]) (Subscript[v, z])(Sqrt[(Subscript[v, x]^2)+(Subscript[v, z]^2)])-g,         z'==(Subscript[v, z]),

initial conditions

x[0]==0,                         z[0]==0
(Subscript[v, x])[0]== (Subscript[v, 0])Cos[θ],      (Subscript[v, z])[0]== (Subscript[v, 0])Sin[θ],

where θ is the shooting angle, and the final conditions

x[(Subscript[t, f])]==L,           z[(Subscript[t, f])]==0

The system of DEs is nonlinear and cannot be solved analytically in general. Also one cannot immediately find the trajectory numerically because θ is unknown. The method of solution of such problems is to find a numerical solution of this problem with an arbitrary θ and then adjust θ to hit the target. So first we define the trajectory numerically as a function of θ by first defining θ-dependent solution of the DEs

Parameters={L->1,g->1,(Subscript[v, 0])->1.3,(Subscript[α, 2])->0.5};
Sol[θ_]:=NDSolve[{
vx'[t]==-(Subscript[α, 2])vx[t](Sqrt[((vx[t])^2)+((vz[t])^2)]),             x'[t]==vx[t],
vz'[t]==-(Subscript[α, 2])vz[t](Sqrt[((vx[t])^2)+((vz[t])^2)])-g,        z'[t]==vz[t],
x[0]==0,                         z[0]==0,
vx[0]== (Subscript[v, 0])Cos[θ],      vz[0]== (Subscript[v, 0])Sin[θ]
}/.Parameters,{vx,x,vz,z},{t,0,((2(Subscript[v, 0])Sin[θ])/g)/.Parameters}]

and then defining the trajectories

xt[t_,θ_]:=x[t]/.Sol[θ][[1]];
zt[t_,θ_]:=z[t]/.Sol[θ][[1]];

Check that it works

xt[1,π/4]
zt[1,π/4]

In plotting these functions one cannot use t because this interferes with t used by NDSolve.

θθ=π/4;
ParametricPlot[{xt[tt,θθ],zt[tt,θθ]},{tt,0,((2(Subscript[v, 0])Sin[θθ])/g)/.Parameters}]

Now we define the final (landing) time as a function of the angle

tf[θ_]:=t/.FindRoot[zt[t,θ]==0,{t,((2(Subscript[v, 0])Sin[θ])/g)/.Parameters}]

Check this:

θθ=π/6;
zt[tf[θθ],θθ]

OK. One can plot xt[tf[θ], θ] vs θ and solve the equation xt[tf[θ], θ]==L graphically:

Plot[{xt[tf[θ],θ],L/.Parameters},{θ,0,π/2}]

However, one cannot find the aiming angles with FindRoot (as well as one cannot find the maximum of xt[tf[θ], θ] with FindMaximum):

(Subscript[ϑ, 1])=(1/2)ArcSin[((g L)/(Subscript[v, 0]^2))]/.Parameters  ;     
(Subscript[ϑ, 2])=(1/2)(π-ArcSin[((g L)/(Subscript[v, 0]^2))])/.Parameters ;
(Subscript[θ, 1])=θ/.FindRoot[xt[tf[θ],θ]==L/.Parameters,{θ,(Subscript[ϑ, 1])}]
(Subscript[θ, 2])=θ/.FindRoot[xt[tf[θ],θ]==L/.Parameters,{θ,(Subscript[ϑ, 2])}]

(Complaints..)

As the function xt[tf[θ], θ] itself is good, one can bypass Mathematica bugs by making its interpolated copy. In spite of error messages, this procedure works.

xtInterpolated=FunctionInterpolation[xt[tf[θ],θ],{θ,0,π/2}]

Plot[{xtInterpolated[θ],L/.Parameters},{θ,0,π/2}]

Now one can find the maximal horizontal distance and the corresponding aiming angle for the parameters defined above

FindMaximum[xtInterpolated[θ],{θ,(π/4)}]
(Subscript[θ, max,degree])=(θ/.FindMaximum[xtInterpolated[θ],{θ,(π/4)}][[ 2]])(180/π)

as well as the two aiming angles for the given distance

(Subscript[ϑ, 1])=(1/2)ArcSin[((g L)/(Subscript[v, 0]^2))]/.Parameters  ;     
(Subscript[ϑ, 2])=(1/2)(π-ArcSin[((g L)/(Subscript[v, 0]^2))])/.Parameters ;
(Subscript[θ, 1])=θ/.FindRoot[xtInterpolated[θ]==L/.Parameters,{θ,(Subscript[ϑ, 1])}]
(Subscript[θ, 2])=θ/.FindRoot[xtInterpolated[θ]==L/.Parameters,{θ,(Subscript[ϑ, 2])}]

Here are the uncorrected and exact aiming angles, respectively, in degrees

(Subscript[ϑ, 1,Degree])=(Subscript[ϑ, 1])(180/π)
(Subscript[θ, 1,Degree])=(Subscript[θ, 1])(180/π)

(Subscript[ϑ, 2,Degree])=(Subscript[ϑ, 2])(180/π)
(Subscript[θ, 2,Degree])=(Subscript[θ, 2])(180/π)

Let us now plot the trajectories with corrected and uncorrected aiming angles, as well as trajectories without the air resictance

x0[t_,θ_]:=(Subscript[v, 0])Cos[θ]t;        z0[t_,θ_]:=-(1/2)g (t^2)+(Subscript[v, 0])Sin[θ]t;      (* Without air resistance *)
Show[
ParametricPlot[{x0[t,(Subscript[ϑ, 2])],z0[t,(Subscript[ϑ, 2])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 2])]))/.Parameters},PlotStyle->{Red,Dashed}],
ParametricPlot[{x0[t,(Subscript[ϑ, 1])],z0[t,(Subscript[ϑ, 1])]}/.Parameters,{t,0,(L/((Subscript[v, 0])Cos[(Subscript[ϑ, 1])]))/.Parameters},PlotStyle->{Blue,Dashed}],
ParametricPlot[{xt[tt,(Subscript[θ, 2])],zt[tt,(Subscript[θ, 2])]},{tt,0,tf[(Subscript[θ, 2])]},PlotStyle->{Thick,Red}], (* corrected *)
ParametricPlot[{xt[tt,(Subscript[θ, 1])],zt[tt,(Subscript[θ, 1])]},{tt,0,tf[(Subscript[θ, 1])]/.Parameters},PlotStyle->{Thick,Blue}], (* corrected *)
ParametricPlot[{xt[tt,(Subscript[ϑ, 2])],zt[tt,(Subscript[ϑ, 2])]},{tt,0,tf[(Subscript[ϑ, 2])]},PlotStyle->{Red}], (* uncorrected *)
ParametricPlot[{xt[tt,(Subscript[ϑ, 1])],zt[tt,(Subscript[ϑ, 1])]},{tt,0,tf[(Subscript[ϑ, 1])]},PlotStyle->{Blue}] (* uncorrected *)
]

Stability

Solutions of DEs and systems of DEs in some cases approach stationary regimes with time. In the simplest case the stationary regime is time independent. For instance, a free pendulum with friction will come to its lowest point with time and will remain there forever. If a periodic force is applied to a pendulum with friction, it will after some time sync with the periodic force. This is an example of a time-dependent stationary regime.

The stationary regimes considered above are absolutely stable. The system will enter these regimes with time starting from any initial state. There are regimes stable only with respect to small perturbations, such as a system being at a local (but not global) minimum of its potential energy. After a small deviation from this state the system will return to it. But after a larger deviation the system might go to another energy minimum.

There are unstable stationary solutions of DEs. For instance, a particle left at the point of its maximal potential energy with a zero speed will remain there forever. But after an arbitrarily small deviation the particle will leave this point.

Pendulum

As an example consider a pendulun with a viscous friction

ϕ''[t]+(αϕ)'[t]+(Subscript[ω, 0]^2)Sin[ϕ[t]]==0 .

In the vicinity of the potential energy minimum, ϕ=0, this DE can be linearized in ϕ

ϕ''[t]+(αϕ)'[t]+(Subscript[ω, 0]^2)ϕ[t]==0 .

Its solution with the initial condition ϕ[0]=(Subscript[ϕ, 0]) and, say, ϕ'[0]=0 is given by

ϕt[t_]=ϕ[t]/.DSolve[{ϕ''[t]+α ϕ'[t]+(Subscript[ω, 0]^2)ϕ[t]==0,ϕ[0]==ϕ0,ϕ'[0]==0},ϕ,t][[1]]

The real parts of the exponentials are negative, thus the solution decreases with time and the pendulum returns to its stable stationary state

Plot[ϕt[t]/.{(Subscript[ω, 0])->1,α->0.1,ϕ0->0.01},{t,0,(15/α)/.α->0.1},PlotRange->All]

In this case, the stability is entirely a property of the potential energy that has a minimum at ϕ=0. To demonstrate this stability, one does not need to know the exact dynamics of the oscillator. For instance, one can consider such a big damping α that the second derivative in the equation of motion can be neglected. This yields the first-order linear relaxational equation of motion

ϕ'[t]==-Γϕ[t],

where Γ≡(Subscript[ω, 0]^2)/α. The relaxational DE is one of the most popular in physics, and a system it describes can be called linear relaxator. The relaxator evolves to zero with time as

ϕ[t]=(Subscript[ϕ, 0])(e^(-Γt)),

thus it is stable. Let us now consider the dynamics of the pendulum near the maximum of its potential energy at  ϕ=π/2. In this case the DE can be linearized in the small deviation δϕ=ϕ-π/2

(δϕ)''[t]+(αδϕ)'[t]-(Subscript[ω, 0]^2)δϕ[t]==0 .

Its solution with the initial condition δϕ[0]=(Subscript[ϕ, 0]) and, say, δϕ'[0]=0 is given by

Clear[α,ϕ0]
δϕt[t_]=δϕ[t]/.DSolve[{(δϕ)''[t]+α (δϕ)'[t]-(Subscript[ω, 0]^2)δϕ[t]==0,δϕ[0]==ϕ0,δϕ'[0]==0},δϕ,t][[1]]

The real parts of one of the exponentials is positive, thus the solution increases with time and the pendulum goes away from its unstable stationary state

Plot[δϕt[t]/.{(Subscript[ω, 0])->1,α->0.1,ϕ0->0.01},{t,0,(1/α)/.α->0.1},PlotRange->All]

Again, if we consider very large α, we obtain the first-order equation ϕ'[t]==Γϕ[t] that describes exponential divergence.

Thermal runaway model

As a more interesting example of stability vs instability in solutions of DEs consider the simplest model of the thermal runaway or the explosive instability, This phenomenon can occur, in particular, in explosive's depots in the absence on a sufficient cooling. Cases of explosive instability were frequent in the Soviet Union during the WW-II and physicists (some in detention) had been asked to provide a theory. Important contributions were made, in particular, by Semenov and Frank-Kamenetskii. The very first and most primitive Semenov's model uses the first-order DE for the temperature of the explosive

C((dT)/(dt))  ==Q[T],    Q[T]= (Subscript[Q, Reaction])[T]-(Subscript[Q, Cooling])[T].

Here C is the heat capacity, (Subscript[Q, Reaction])[T] describes the heat-release rate due to the chemical reaction (decay of the explosive), and  (Subscript[Q, Cooling])[T] is the cooling rate due to the contact with the environment. The rate of the chemical reaction strongly increases with temperature and can be well approximated by the Arrhenius dependence

(Subscript[Q, Reaction])[T]=(Subscript[Q, 0])Exp[-((ΔU)/((Subscript[k, B])T))] ,

where ΔU is the energy barrier of the metastable state of the explosive that prevents its immediate decay. (Subscript[Q, Cooling])[T] can be approximated by

(Subscript[Q, Cooling])[T]=κ(T-(Subscript[T, 0])) ,

where (Subscript[T, 0]) is the temperature of the environment and κ is an effective thermal conductivity between the explosive and the environment.  The stability-instability of Semenov's temperature equation is explained in the plot below.

At the temperatures not much above (Subscript[T, 0]) the reaction rate (Subscript[Q, Reaction])[T] grows with temperature usually slower than the cooling rate (Subscript[Q, Cooling])[T] (case A). Thus the temperature increases from (Subscript[T, 0]) until the two rates become equal at T=(Subscript[T, 1]). This the stable stationary state of the system. However, at higher T the reaction rate increases much faster and the two rate curves cross again at (Subscript[T, 2]). This is the unstable stationary state. For the starting temperatures (Subscript[T, 1])<T<(Subscript[T, 2])  one has Q[T]<0 and the system cools down to (Subscript[T, 1]). But for the starting temperatures  T>(Subscript[T, 2]) one has Q[T]>0 and the system heats up infinitely (thermal runaway). If cooling is too weak (case B) there are no stationary states, one has Q[T]>0 at all temperatures and thermal runaway occurs unconditionally. A recent example a thermal runaway is the Bhopal disaster of 2-3 December 1984 when an explosion at a chemical plant spread pesticides over a densely populated area.

It is useful to represent Q[T] in the form Q=-(∂_T)W[T], where W plays the role of an effective potential energy. The latter is given by

Clear[ΔU,kB,κ,T0]
W[T_]=-∫(Q0 Exp[-((ΔU)/(kB T))]-κ(T-T0))\[DifferentialD]T

ΔU=10;kB=1;Q0=100;T0=1;  κ=2;(* Parameters have to be carefully chosen *)
Plot[W[T],{T,0,5}]

Semenov's equation for the temperature in the form

C((dT)/(dt))  ==-((dW)/(dT))

is mathematically equivalent to that for an overdamped particle sliding down the potential energy W[T].

Semenov's model does not take into account that the amount of the explosive decreased as a result of burning. Thus this model yields an infinite increase of the temperature in the thermal runaway. A more realistic model uses two first-order differential equations, one for the temperature T and one for the mass of the fuel (explosive) m

C((dT)/(dt))  ==-(Subscript[E, m])((dm)/(dt)) -κ(T-(Subscript[T, 0]))
((dm)/(dt))=-Γ[T]m,          Γ[T]=(Subscript[Γ, 0])Exp[-((ΔU)/((Subscript[k, B])T))].

Here (Subscript[E, m]) is the energy released by burning of one unit mass of the fuel and Γ[T] is the burning rate of the fuel. Substituting the second equation into the first equation, one obtains Semenov's equation with (Subscript[Q, Reaction])[T]=(Subscript[E, m])m Γ[T]. To recover Semenov's theory, one has to fix m to its initial value, to make the equation closed. Below is the numerical solution

ΔU=3000;   T0=300;    Em=3000;   c=1;    
Γ0=20;   kB=1;   m0=1;      (* κ=0.189; *)
Γ[T_]:=Γ0 Exp[-((ΔU)/(T kB))]
tMax=500;
Eqm=m'[t]==-Γ[T[t]]m[t];
EqT[κ_]=c T'[t]+Em m'[t]==-κ(T[t]-T0);
sol[κ_]:=NDSolve[{Eqm,EqT[κ],m[0]==m0,T[0]==T0},{m,T},{t,0,tMax}];
mt[κ_,t_]:=m[t]/.First[sol[κ]]
Tt[κ_,t_]:=T[t]/.First[sol[κ]]

For κ=0.188 there is a thermal runaway with fast and complete burning but for κ=0.189 there is only a slow burning. Such dependence of the solution of a differential equation on a parameter is called bifurcation.

Plot[{mt[0.188,tt],mt[0.189,tt]},{tt,0,tMax},PlotRange->{0,1},PlotStyle->{{Thick,Red},{Thick,Blue}}]       (* cannot use t *)

In the case of thermal runaway temperature peaks and then returns to normal after all fuel is burned.

Plot[{Tt[0.188,tt],Tt[0.189,tt]},{tt,0,0.4tMax},PlotRange->All,PlotStyle->{{Thick,Red},{Thick,Blue}}]     (* cannot use t *)