First-Order Ordinary Differential Equations

Background--simple iteration

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

ExampleFunction[i_,  alpha_]  :=  
      ExampleFunction[i - 1, alpha]  +  
        alpha*ExampleFunction[i - 1, alpha]

What would happen if the following functions were evaluated?

ExampleFunction[25, 0.25] 

    ExampleFunction[25, -0.125]

It is clear that the function needs some value at some time (an initial condition) from which it obtains all its other values:

ExampleFunction[0, 0.25]  =  π/4

ExampleFunction[18, 0.25]

This function is fairly inefficient, it calculates the the first value about 18 times before it returns the answer, the following function sets the initial value and remembers what "trajectory" it is on--i,e., the trajectory is determined from the start

GrowsPropToSize[i_ ,  alpha_ ,  F0_]  :=  
      If[i == 0, GrowsPropToSize[0, alpha, F0] = F0, 
        GrowsPropToSize[i, alpha, F0]  =  
          GrowsPropToSize[i - 1, alpha, F0]  +  
            alpha*GrowsPropToSize[i - 1, alpha, F0]  ]

For example:

Table[GrowsPropToSize[i,  .01, 1], {i, 18}]

Might as well write a function to calculate a vector representing trajectory, here is an implementation to compute a trajectory vector of length 300

Trajectory[alpha_,  F0_]  :=  
      Table[GrowsPropToSize[i, alpha, F0], {i, 300}]

ListPlot[Trajectory[0.01, 1], PlotJoined -> True]

<< Graphics`MultipleListPlot`

Trajectory[alpha_,  F0_]  :=  
      Table[GrowsPropToSize[i, alpha, F0], {i, 300}]

Plotting a bunch of curves, each belonging to a different initial value, demonstrates that a point in space "belongs" to a curve.

MultipleListPlot[Trajectory[ .01, - .5], Trajectory[ .01,  .5], 
      Trajectory[ .01, 1], Trajectory[ .01, 1.5], 
      PlotJoined -> True]

MultipleListPlot[Trajectory[- .01, - .5], 
      Trajectory[- .01,  .5], Trajectory[- .01, 1], 
      Trajectory[- .01, 1.5], PlotJoined -> True]

Background- Forward Differencing

Create a function to return a list of values from forward differencing with a function (These examples were modified from those found in ``Computer Science with Mathematica'' by Roman E. Maeder, Cambridge University Press, (2000)

ForwardDifferenceListV1[AFunction_ ,  InitialValue_ ,  delta_,  
        ListLength_Integer]  := 
Module[{TheResultList  =  {InitialValue}, TheValue = InitialValue}, 
Do[
          TheValue  =  TheValue  +  delta AFunction[TheValue];  
          AppendTo[TheResultList, TheValue], {ListLength}]; 
        TheResultList]

exampleFunction[x_]  :=  0.1 x

result  =  
        ForwardDifferenceListV1[exampleFunction,  1,  0.01,  
          500];

result // Short

ListPlot[result, PlotJoined -> True]

Write another version of this forward difference function that returns a list of values
and the "x" value for subequent use in ListPlot, this one will take (Subscript[x, o]) and y((Subscript[x, o])) as an argument in a list

Clear[ForwardDifferenceListV2]

ForwardDifferenceListV2[AFunction_ , x0_ ,  fx0_ ,  delta_,  
        Xlast_]  := 
Module[{TheResultList  =  {{x0, 
                fx0}}, TheValue = fx0, 
          CurrentX   =  x0}, 
While[
          CurrentX  <  Xlast, 
CurrentX  = 
            CurrentX  +  delta; 
 
          TheValue =  TheValue  +  delta AFunction[TheValue];  
          AppendTo[
            TheResultList, {CurrentX, 
              TheValue}]]; 
TheResultList]

result  =  
        ForwardDifferenceListV2[exampleFunction, 0, 1,  0.01, 
          4];

result // Short

ListPlot[result]

Background-Repeated Operations on a Function

Because each iteration is the same, the iteration can be considered a functional operation, for the the case considered above,
{(Subscript[x, i + 1]), (Subscript[y, i + 1])} = {(Subscript[x, i])+ δ, (Subscript[y, i])+ δ f((Subscript[y, i]))}.  Therefore a "Incrementing Operator" can be obtained that updates the values:

StepOnce[{x_, y_},  AFunction_,  delta_]  :=  {x +  delta,  
        y  +  delta AFunction[y]}

Then, the trajectory should be obtained from:
    {{xo,yo},StepOnce[{xo,yo}], StepOnce[StepOnce[{xo,yo}]], ⃛⃛⃛}
    This is what NestList[function, initialvalue,depth] does:

OurStepOnce[{x_ ,  y_}]  :=  
      StepOnce[{x, y}, exampleFunction,  0.01]

result  =  NestList[OurStepOnce, {0, 1}, 400];

ListPlot[result]

Using a Mathematica trick of a ``pure function'' one can eliminate the intervening function (OurStepOnce) definition:

ListPlot[
      NestList[StepOnce[#, exampleFunction, 0.01] &, {0, 1}, 400]]

DisplayLater  =  
        DisplayFunction -> Identity;
DisplayNow  =  
        DisplayFunction -> $DisplayFunction;

lp[i_]  := 
        ListPlot[NestList[StepOnce[#, exampleFunction, 0.01] &, {0, i}, 400], 
          DisplayLater];

Show[Table[lp[i], {i, -4, 4,  .5}], DisplayNow]

Visualizing Geometry 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-Θ

\(Plot3D[
      1 - Θ, \ {tau, \ \(-1\), 
        1}, \ {Θ, \ \(-2\), 3}, 
      AxesLabel -> {"τ", "Θ", \*""d(Θ/d)τ""}]\)

Plot the vector (dτ,dΘ) = dτ(1, d(Θ/d)τ) at each point:

<< Graphics`PlotField`

PlotVectorField[{1, 1 - Θ}, {tau,  0, 
        4},  {Θ, -2, 4}, Axes -> True, 
      AxesLabel -> {"τ", "Θ"}]

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

\(Plot3D[y\ \ Sin[\ \(y\ t\)\/\(t + y + 1\)], {t, 0, 10}, {y, 0, 10}, 
      AxesLabel -> {"t", "y", \*""\!\(dy\/dt\)""}, 
      PlotPoints -> 40]\)

PlotVectorField[{1, y  Sin[ ((y t)/(t + y + 1))]}, {t,  0, 
        12},  {y, 0, 10}, Axes -> True, 
      AxesLabel -> {"t", "y"}]

Using Mathematica's DSolve function

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

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

exactplot  = Plot[y[x] /. Flatten[dsol], {x, 0, 10}]

Compare this to the forward differencing scheme:
Recall how
exampleFunction[x_]  :=  0.1 x
StepOnce[{x_, y_},  AFunction_,  delta_]  :=  {x +  delta,  
          y  +  delta AFunction[y]}
OurStepOnce[{x_ ,  y_}]  :=  
        StepOnce[{x, y}, exampleFunction,  0.01]
result  =  NestList[OurStepOnce, {0, 1}, 400]; 
ListPlot[result]
were used to create a plot for a function that grew at a rate that was 0.1 times its current size, try this for the ODE that was solved for above, let the new function be y':

NewExampleFunction[x_ ,  y_]  :=  -Sin[x] y

NewStepOnce[{x_,  y_},  AFunctionXY_ ,  
        delta_]  :=  {x  +  delta,  
        y  +  delta AFunctionXY[x, y]}

NewOurStepOnce[{x_,  y_}]  :=  
      NewStepOnce[{x, y}, NewExampleFunction, 0.01]

result  =  NestList[NewOurStepOnce, {0, 1}, 1000];

result // Short

forwarddifferenceplot  = 
      ListPlot[result, PlotStyle -> {Hue[1]}]

Show[forwarddifferenceplot, exactplot]