First-Order Ordinary Differential Equations

Background--Simple iteration from a Fixed Boundary Condition

Suppose a function,  F[i], changes proportional to its current size, i.e.,  F[i+1] = F[i] + αF[i]

ExplFun[i_,α_] := ExplFun[i,α] =ExplFun[i-1,α] + α*ExplFun[i-1,α]

In the above, the symbol is assigned (ExplFun[i,alpha] = ...) as part of the function definition, so that intermediate values are ``remembered.''
The function needs some value at some time (an initial condition) from which it obtains all its other values:

ExplFun[0,0.25] = π/4

ExplFun[18,0.25]

Background--Simple iteration from a Specified Boundary Condition

ExplFun[0,0.25] = π/4 above serves as an initial value for the function. The initial value and α determine the value at any later time.  The initial value can be expressed as another parameter for the function. We use a recursive method to define the function and store computed values in memory.

ExplFun[0 , α_, InitVal_] := InitVal
ExplFun[Inc_Integer, α_, InitVal_] := ExplFun[Inc,α, InitVal] =ExplFun[Inc-1,α, InitVal] + α*ExplFun[Inc-1,α,InitVal]

One might as well write a function to calculate a vector representing `trajectory', here is an implementation to compute a trajectory vector of specified length.  We use AppendTo to grow an initially empty list.  We will do the same thing with another method (Nest) below.

Traj[Steps_Integer?Positive,α_, InitVal_] :=Traj[Steps,α,InitVal]= AppendTo[Traj[Steps-1,α,InitVal],ExplFun[Steps,α,InitVal]]
Traj[0,_,_] = {};

Traj[12,.01,.001]

We define a function, Evolve, producing  an interactive tool with input : Initial values and α.

Evolve[NumSteps_Integer/;(NumSteps>0 && NumSteps < 10000)] := 
Manipulate[Block[{$RecursionLimit=4*NumSteps},UpBound=
N[Last[Traj[NumSteps,α,InitVal]]];
yrange={Min[{0,InitVal,1.05UpBound}],Max[{0,InitVal,1.05UpBound}]};
ListPlot[Traj[step,α, InitVal],PlotRange->{{0,NumSteps},yrange},PlotStyle->{Red,PointSize[0.005]},BaseStyle->{Large},ImageSize->Full]],{{step,Quotient[NumSteps,10]},1,NumSteps,1,Appearance->"Labeled"},{{α,0.001},-1,1,Appearance->"Labeled"},{{InitVal,.01},-1,1,Appearance->"Labeled"}]

Evolve[300]

Visualizing Trajectories: Initial Conditions Define a Set of Curves that "Cover" the Plane.

Plotting a bunch of curves for the same positive α value, but each corresponding to a different initial value, demonstrates that a point in space "belongs" to a curve. Here we make a function to demonstrate for a range of inital values and a function of α

PlotTrajs[α_] := Block[{$RecursionLimit=10^4},ListPlot[Evaluate[Table[Traj[300,α,iv],{iv,-1,1,0.25}]],PlotRange->All,Joined->True,ImageSize-> Large,BaseStyle-> Large,PlotStyle->Thick]]

PlotTrajs[0.02]

Making a similar plot for a family of curves corresponding to a negative α value

PlotTrajs[-0.02]

Finite Differencing Methods (Method 1): Explicit Forward Differences

Approximation to (dy/dt)= f(y)
Approximate (dy/dt)= f(y) with  a finite difference ((y(t + Δt) - y(t))/((t + Δt)) - t))= f(y(t))
If all the time steps are the same size, Δt, then use the short-hand (Subscript[y, n]) = y(n × Δt)

Approximate f(y) with f((Subscript[y, i-1])); then solving the finite difference equation above,
(Subscript[y, i])= (Subscript[y, i-1 ])+ Δt f((Subscript[y, i-1])): the right-hand side depends only on (Subscript[y, i-1]), and therefore we can apply a recursion method to find any (Subscript[y, i])

PushMethod1 takes an arbitrary f for the right-hand-side, the current point and the size of the time step Δt

PushMethod1[f_, {ti_,yi_}, Δt_]:={ti + Δt, yi +Δt f[yi]}

Here is an example right-hand-side function:

FuncEx[y_] := -Sin[y]

Applying it once "pushes" the current point (the second argument) forward one increment in time.

PushMethod1[FuncEx,{0,1},.01]

Applying it twice "pushes" the current point (the second argument) forward two increments in time. We have "nested" one function call inside the other.

PushMethod1[FuncEx,PushMethod1[FuncEx,{0,1},.01],.01]

Nest is a useful function for these kinds of operatons: Nest[g,x,3] is g[g[g[x]]]. However, in our case, PushMethod1 takes three arguments, we need to use a Pure Function methods to specify that the PushMethod1 operates only on the point represented by the second argument.  This should provide a useful example of what Nest does: # plays the role of the functions argument and the & indicates that the function is to interpreted as a pure function.

Nest[PushMethod1[FuncEx,#,.01]&,{0,1},2]

We can use NestList to build up the trajectory,very similar to our method above that used AppendTo, but this is more efficient.

NestList[PushMethod1[FuncEx,#,.01]&,{0,1},2]

NestListWhile to build up the trajectory while a "test" evaluates to True.  We have to use a pure function for the test, which operates on the current point (or the last point entered into the list).

NestWhileList[PushMethod1[FuncEx,#,.01]&,{0,1},(First[#]< 0.03)&]

Visualizing Trajectories from the  Explicit Forward Differences Method

A Function PlotM1, taking arguments for α and initial condition is used with NestListWhile and ListPlot to produce graphics with a red line and green points.

PlotM1[Δt_,InitCond_]:=PlotM1[Δt,InitCond]=ListPlot[NestWhileList[PushMethod1[FuncEx,#,Δt]&,{0,InitCond},(First[#]<20)&],Joined->True,PlotStyle-> {Thickness[0.01],Red},PlotRange->All,Mesh->Full,MeshStyle-> {Green,PointSize[0.01],Opacity[0.5]},ImageSize-> Full]

PlotM1[0.1, 1]

PlotM1[1.5,1]

Finite Differencing Methods (Method 2): Implicit Forward Differences

Approximate f(y) with f((Subscript[y, i])); then solving the finite difference equation above,
(Subscript[y, i])= (Subscript[y, i-1 ])+ Δt [f((Subscript[y, i])) ].  So, (Subscript[y, i])= (Subscript[y, i-1]) + Δt (f((Subscript[y, i])) + f'((Subscript[y, i-1]))dy)
 (Subscript[y, i])= (Subscript[y, i-1]) + Δt (f((Subscript[y, i])) + f'((Subscript[y, i-1]))((Subscript[y, i]) - (Subscript[y, i-1])) )
(Subscript[y, i]) = ((Subscript[y, i-1]) - Δt [f((Subscript[y, i-1]))  - f'((Subscript[y, i-1]))(Subscript[y, i-1])])/(1 -  Δt  f'((Subscript[y, i-1])) )

PushMethod2[f_,df_, {ti_,yi_}, Δt_]:={ti + Δt, (yi +(Δt (f[yi] - df[yi] yi)))/(1 - Δt df[yi])}

dFuncEx[y_] :=Evaluate[D[FuncEx[y],y]]

NestList[PushMethod2[FuncEx,dFuncEx,#,0.1]&,{0,1},3]

We define a function, PlotM2, which takes arguments Δt and InitialCondition and then uses ListPlot with Blue lines and Gray points.

PlotM2[Δt_,InitCond_]:=PlotM2[Δt,InitCond]=ListPlot[NestWhileList[PushMethod2[FuncEx,dFuncEx,#,Δt]&,{0,InitCond},(First[#]< 20)&],Joined->True,PlotStyle-> {Thick,Blue},PlotRange->All,Mesh->Full,MeshStyle-> {Black,PointSize[0.01],Opacity[0.5]},ImageSize-> Full]

Visual Comparison of Trajectories from the  Explicit and Implicit Methods

funcextplot = Show[PlotM1[0.1,1],PlotM2[0.1,1],PlotRange->{{0,10},{0,1}}]

Show[PlotM1[1.5,1.0],PlotM2[1.5,1.0],PlotRange->{{0,10},{-0.5,1}}]

Method 2 will fail if the step size is increased to 2

Intuitive Understanding of First-Order ODEs

Newton's law of cooling (dT/dt)= -k(T - (Subscript[T, o])) can be written in the non-dimensional form  ((dΘ)/(dτ))= 1-Θ  (i.e, defining Θ as T/(Subscript[T, o]) and τ =  k t)

In the general case, ((dΘ)/(dτ)) will depend on both Θ and t, i.e., ((dΘ)/(dτ)) = ((dΘ)/(dτ))(Θ,t). This the equation of a surface in three dimensions, as shown in the following plot (in this specific case there is no t dependence):

ZeroPlane[xmin_, xmax_,ymin_,ymax_]:= Graphics3D[{Gray,Opacity[0.25],Cuboid[{xmin,ymin,-.001},{xmax,ymax,.001}]}]

Show[Plot3D[1-Θ, {tau,0,2}, {Θ, 0,2},AxesLabel->{"τ","Θ","((dΘ)/(dτ))"},DisplayFunction->Identity],ZeroPlane[0,2,0,2]]

Geometry of Flows of First-Order ODEs: Vector Fields (Simple Example: Newton's Law of Cooling)

Plot the vector (dτ,dΘ) = dτ(1, ((dΘ)/(dτ))) at each point:
We want to make a plot of  ((dΘ)/(dτ))in the Θ-τ plane. We can do so by plotting vectors of the form {dτ, dΘ} = dτ{1, ((dΘ)/(dτ))} which will be proportional to the vector {1, 1-Θ}. This is done as follows:

GraphicsRow[
{
VectorPlot[{1,1-Θ},{tau,0,4},{Θ,-2,4},Axes->True,AxesLabel->{"τ","Θ"},ImageSize-> Large],
StreamPlot[{1,1-Θ},{tau,0,4},{Θ,-2,4},Axes->True,AxesLabel->{"τ","Θ"},ImageSize-> Large]
}
]

Show[
ContourPlot[Norm[{1,1-Θ}],{tau,0,4},{Θ,-2,4},Axes->True,FrameLabel->{"t","normalized temperature"},ImageSize-> Full,BaseStyle-> Large,ColorFunction-> "Temperature"],StreamPlot[{1,1-Θ},{tau,0,4},{Θ,-2,4},Axes->True,FrameLabel->{"t","y"},ImageSize-> Full]
]

Note in the equations and the plot that((dΘ)/(dτ))is independent of τ.

Geometry of Flows of First-Order ODEs: Vector Fields (More Complex Example)

Slightly more complicated example: (dy/dt)= y sin((yt/(1 + t +y))),
(dt,dy) = dt(1,ysin(yt/(1 + t +y))))

As in the simpler example above we first plot the derivative's value as a surface using Plot3D, then we represent the behavior as a vector field

Show[Plot3D[y  Sin[ ((y t)/(t+y+1))],{t,0,10},{y,0,10},AxesLabel->{"t","y","(dy/dt)"},PlotPoints->40,BaseStyle-> Large,ImageSize-> Large],ZeroPlane[0,10,0,10]]

GraphicsRow[
{
VectorPlot[{1,y Sin[((y t)/(t+y+1))]},{t,0,10},{y,0,10},Axes->True,AxesLabel->{"t","y"},ImageSize-> Full],StreamPlot[{1,y Sin[((y t)/(t+y+1))]},{t,0,10},{y,0,10},Axes->True,AxesLabel->{"t","y"},ImageSize-> Full],
}
]

Show[
ContourPlot[Norm[{1,y Sin[((y t)/(t+y+1))]}],{t,0,10},{y,0,10},Axes->True,FrameLabel->{"t","y"},ImageSize-> Full,BaseStyle-> Large,ColorFunction-> "Temperature"],
StreamPlot[{1,y Sin[((y t)/(t+y+1))]},{t,0,10},{y,0,10}]
]

To summarize this part, we have seen how first-order differential equations have associated trajectories in the plane of the dependent variables, analogous to what we simulated above by incrementing a certain functional relationship. Numerical differential equation solvers use similar techniques.

Using Mathematica's DSolve function on Example used for Finite Differencing Examples

Mathematica's DSolve function finds closed-form expressions for solutions to differential equations, when possible. There is a related function called NDSolve that finds numerical solutions that can be used when DSolve gives up. Here are some examples of using DSolve:

funcextplot

dsol = DSolve[{y'[t]==FuncEx[y[t]]},y[t],t]

Note that the solution is given as a rule, just like for the function Solve. Because no initial condition was specified, the solution involves an unknown constant, C[1].   In the next  case an initial condition was specified for the differential equation, so there is no undetermined constant in the solution.

dsol = DSolve[{y'[t]==FuncEx[y[t]],y[0]==1},y[t],t]

The next statement extracts y (x) for plotting ...

ExactPlot =Plot[ y[t]/.dsol,{t,0,10},PlotStyle->{Thickness[0.025],Darker[Cyan]},PlotRange->All]

Show[ExactPlot,funcextplot,ImageSize-> Large]

Comparison of Exact Solution to Finite Differencing Methods

Show[PlotM1[0.1,1],PlotM2[0.1,1],ExactPlot,PlotRange->{{0,10},{0,1}}]

Show[PlotM1[1.5,1],PlotM2[1.5,1],ExactPlot,PlotRange->{{0,10},{-0.25,1}}]

DSolve on a First-Order ODE: More Complex Example

Here is another example:

dsol = DSolve[ y'[x]  +x *y[x] ==0,y[x],x]

We can insert boundary conditions other than y[0]

dsol = DSolve[ {y'[x]  + Sin[x] *y[x] ==0,y'[5]==1},y[x],x]

GraphicsRow[{p=Plot[y[x]/.dsol,{x,0,10},PlotStyle->Thick,ImageSize-> Medium],Show[p,PlotRange-> {{0,6}, {0,6}},AspectRatio->1,ImageSize-> Medium]}]