Solving Differential Equations

The Equation To Be Studied

Let us consider the Duffing equation:

\(TraditionalForm\`\(x\&¨\)(t)\  + \ 
          2  γ\ \(\(x\& . \)(t)\)\  + x(t)\  + \ \(x\^3\)(t)\  = \ 
        f\ \(cos(ω\ t)\)\)

exDuffing[g_, f_, w_] := x''[t] + 2  g x'[t] + x[t] + x[t]^3 - 
        f Cos[w t]

exDuffing[g, f, w]

and it's linearization

\(TraditionalForm\`\(x\&¨\)(t)\  + \ 
          2  γ\ \(\(x\& . \)(t)\)\  + x(t)\  = \ 
        f\ \(cos(ω\ t)\)\)

exDuffLin[g_, f_, w_] := x''[t] + 2  g x'[t] + x[t] - 
        f Cos[w t]

exDuffLin[g, f, w]

Notice that the variable parameters in these definitions are γ (g)—the viscous damping, f (f)—the amplitude and ω (w)—the frequency of the driving force. The choice of "x" as the dependent function of the argument "t" we leave as is.

Notice also that the functions "exDuffing[g,f,w]" and "exDuffLin[g,f,w]" are expressions, not equations; we'll make equations by setting "exDuffing[g,f,w]==0", and "exDuffLin[g,f,w]==0".

A Special Case, Analytically

Consider the special case with γ = f = ω = 0, which leaves

exDuffing[0, 0, 0]

still a non-linear 2nd order differential equation.

DSolve[exDuffing[0, 0, 0] == 0, x[t], t]

ans000 = 
      Simplify[DSolve[exDuffing[0, 0, 0] == 0, x[t], t]]

Mathematica always cites the solution as a substitution rule. Most often, this is useful, so as to use the solution somewhere else. Here, however, we need the solution itself. We also note that Mathematica could have saved some ink by redefining the first integration constant, C[1]:

ans000[[1]]

x[t] /. ans000[[1]]

xSol1[t_] := 
      Simplify[x[t] /. ans000[[1]]] /. {Sqrt[1 + 2  C[1]] -> 
            C[1]}

xSol1[t]

xSol2[t_] := 
      Simplify[x[t] /. ans000[[2]]] /. {Sqrt[1 + 2  C[1]] -> 
            C[1]}

xSol2[t]

This cannot be plotted as it depends on the two integration constants, C[1] and C[2]. Of course, some initial conditions could fix those. Let's aim for x[0]=0, and x'[0]=1.

Simplify[D[xSol1[t], t]]

Simplify[D[xSol1[t], t]] /. {t -> 0}

Solve[{xSol1[t] /. {t -> 0}, 
          Simplify[D[xSol1[t], t]] /. {t -> 0}} == {0, 1}, {C[1], 
        C[2]}]

So already in this, simplest case ( γ = f = ω = 0 ) the analytic solution is not very helpful in trying to see how this looks like.

Numerical Solving

Solving differential equations numerically is of course also possible, and it generally provides a solution starting from the initial conditions—thus explicitly building them into the solution.

The major drawback is that we must provide numerical values to all parameters, since numerical solutions by definition do not exist for symbolic parameters.

Let's choose γ=0.1 (g->0.1), f = 1 (f->1) and ω=4 (w->4):

aNum = 
      NDSolve[{exDuffing[ .1, 1, 4] == 0, x'[0] == 1, 
          x[0] == 0}, x, {t, 0, 10}]

As usual, Mathematica returns the answer as a substitution. The function is then specified as "x[t]/.aNum":

Plot[x[t] /. aNum, {t, 0, 10}]

This plot can be decorated by specifying options:

Plot[x[t] /. aNum, {t, 0, 10}, AxesLabel -> {"t", "x"}, 
      PlotLabel -> "x[t] for x[0]=0 and x'[0]=1"]

To obtain a plot for a different frequency, the process must be repeated:

aNum2 = 
      NDSolve[{exDuffing[ .1, 1, 3] == 0, x'[0] == 1, 
          x[0] == 0}, x, {t, 0, 10}]

Plot[x[t] /. aNum2, {t, 0, 10}, AxesLabel -> {"t", "x"}, 
      PlotLabel -> "x[t] for x[0]=0 and x'[0]=1"]

Plot[x[t] /. aNum2, {t, 0, 10}, AxesLabel -> {"t", "x"}, 
      PlotLabel -> "x[t] for x[0]=0 and x'[0]=1"]

pDuffing[g_, f_, w_, ti_, tf_, x0_, 
        xp0_] := 
Module[{nFun, nPlot}, 
        nFun = NDSolve[{exDuffing[g, f, w] == 0, x'[0] == 
                xp0, x[0] == x0}, 
            x, {t, ti, tf}]; 
nPlot = 
          Plot[x[t] /. nFun, {t, ti, tf}, 
            AxesLabel -> {"t", "x"}, 
            PlotLabel -> "x[t] for x[0]=x0 and x'[0]=xp0"]]

pDuffing[ .05, 1, 1, 0, 70, 0, 0]

pDuffing[ .05,  .5, 1, 0, 70, 0, 0]

pDuffing[ .05,  .1, 1, 0, 70, 0, 0]

Do[pDuffing[ .03, n, 1, 0, 70, 0, 0], {n, 0, 4}]

Do[pDuffing[ .03,  .01 n, 1, 0, 70, 0, 0], {n, 0, 100}]

Do[pDuffing[ .03,  .01 n, 1.5, 0, 70, 0, 0], {n, 0, 100}]

Do[pDuffing[ .03,  .01 n, 2, 0, 70, 0, 0], {n, 0, 100}]