differential_equations_lab

Part 2 - Sensitivity to initial conditions and NDSolve

A differential equation is said to be insensitive to initial conditions (or "stable") if a small change in the initial condition causes a (relatively) small change in the behavior of the solution (over time). It is sensitive to initial conditions if a small change in the initial conditions causes a (relatively) significant change in the behavior of the solution. So, as a simplistic example, consider . This has a solution of $FormBox[RowBox[{y, =, RowBox[{k, , RowBox[{, ^, RowBox[{(, RowBox[{0.4, t}], )}]}]}]}], TraditionalForm]$ , and if we graph this for several closely spaced values of , we get:

RowBox[{{, RowBox[{1, ,, 1.01, ,, 1.02, ,, 1.03, ,, 1.04, ,, 1.05, ,, 1.06, ,, 1.07, ,, 1.08, ,, 1.09, ,, 1.1}], }}]

$RowBox[{Plot, [, RowBox[{RowBox[{Evaluate, [, RowBox[{k, , RowBox[{, ^, RowBox[{(, RowBox[{0.4, t}], )}]}]}], ]}], ,, {t, 0, 5}}], ]}]$

[Graphics:../HTMLFiles/index_37.gif]

As you can see, these solutions draw further and further apart over time. If you graph over a larger time period, this becomes even more evident:

$RowBox[{Plot, [, RowBox[{RowBox[{Evaluate, [, RowBox[{k, , RowBox[{, ^, RowBox[{(, RowBox[{0.4, t}], )}]}]}], ]}], ,, {t, 99, 100}}], ]}]$

[Graphics:../HTMLFiles/index_40.gif]

(Notice, these solutions that started out only covering a range of 0.1, now spread over a range of about $FormBox[RowBox[{RowBox[{0.25, *, 10^17}], =, RowBox[{2.5, *, 10^16}]}], TraditionalForm]$ .) The opposite case can be shown by the equation :

{-10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10}

$RowBox[{Plot, [, RowBox[{RowBox[{Evaluate, [, RowBox[{k, , RowBox[{, ^, RowBox[{(, RowBox[{RowBox[{-, 0.4}], t}], )}]}]}], ]}], ,, {t, 0, 50}, ,, PlotRangeAll}], ]}]$

[Graphics:../HTMLFiles/index_47.gif]

Here, all of the solutions converge towards 0, so this equation clearly is pretty insensitive to initial conditions (stable).

A useful fact (which I'm not even going to try to prove here) is that if you have a first order differential equation of the form , then the equation is stable if ∂f/∂y<0 and unstable if ∂f/∂y>0 .  (Things become a bit more complicated if is sometimes greater and sometimes less than 0.)  Thus, in the first equation, FormBox[RowBox[{∂f/∂y, =, 0.4}], TraditionalForm] , which means it is unstable (as we saw).

Now, this can cause a problem if you are using a numerical method to approximate the solution to a differential equation (i.e., NDSolve).  If the equation is stable, then life is good.  Your approximation should be pretty safe.  However, if the equation is unstable, then each small error in the approximation has the effect of kicking the solution over onto a slightly different solution curve; since the equation is sensitive to initial conditions, that means that when the error moves it over to a new curve, that curve diverges from the original solution curve more and more over time, potentially greatly increasing your total error.  By default, Mathematica tries to make the error for approximating differential equations for each step of the approximation.  If the solutions are stable, then the overall error should remain pretty small.  However, if the solutions are unstable, then those small errors can cause a much larger error, especially over longer time intervals.  So, in our example, if we use NDSolve to approximate the solution and evaluate that at , compared to the actual solution:

RowBox[{approxSolution, =, RowBox[{RowBox[{y[t], /., RowBox[{NDSolve, [, RowBox[{RowBox[{{, Ro ... {0.4, y[t]}]}], ,, y[0] 1}], }}], ,, y[t], ,, {t, 0, 1000}}], ]}]}], /., t1000}]}]

$RowBox[{{, 5.2215*10^173, }}]$

$RowBox[{exactSolution, =, RowBox[{RowBox[{, ^, RowBox[{(, RowBox[{0.4, t}], )}]}], /., t1000}]}]$

$RowBox[{{, 2.96591*10^168, }}]$

I think most people will agree that this is a pretty significant difference (though notice that it will actually be quite accurate for small values). You can correct for this by increasing the accuracy goals for NDSolve. The details of this are a bit complicated, but generally this is pretty easy to do:

RowBox[{approxSolution, =, RowBox[{RowBox[{y[t], /., RowBox[{NDSolve, [, RowBox[{RowBox[{{, Ro ... ;1}], }}], ,, y[t], ,, {t, 0, 1000}, ,, WorkingPrecision16}], ]}]}], /., t1000}]}]

′ NDSolve :: precw : T ... y [t] 0.4` y[t], y[0] 1}) is less than WorkingPrecision (16.`). More…

$RowBox[{{, 5.221498135681*10^173, }}]$

$RowBox[{{, 1.03654*10^165, }}]$

This still isn't great, but it IS better (by several orders of magnitude). For our purposes, WorkingPrecision should generally be about 10 greater than the number of digits you plan to trust (the default value for WorkingPrecision is about 16). The error message here is complaining that the 0.4 in the equation doesn't even come close to the requested WorkingPrecision. You can fix that by replacing all constants in the equation with exact values:

approxSolution = y[t]/.NDSolve[{y '[t] 4/10y[t], y[0] 1}, y[t], {t, 0, 1000}, WorkingPrecision25]/.t1000

$RowBox[{{, 5.221469700129443230662*10^173, }}]$

$RowBox[{{, 3.01831*10^161, }}]$

(Notice that this also makes your result noticeably more accurate as well.) Don't go too crazy with this, however. If you decide to shoot for 25 digits accuracy, for example, you run into the following problem:

approxSolution = y[t]/.NDSolve[{y '[t] 4/10y[t], y[0] 1}, y[t], {t, 0, 1000}, WorkingPrecision75]/.t1000

InterpolatingFunction :: dmval : Input value {1000} lies outside the range of data in the interpolating function. Extrapolation will be used. More…

$RowBox[{{, 1.082403894091420337511553862438689913109891323022925208794651*10^7, }}]$

The error message here is warning you that it isn't able to maintain the requested accuracy over the entire interval. This means the solution is not trustworthy (as you can see by the pretty wildly inaccurate solution it gives). There are ways to deal with this, but I don't plan on covering this (look it up in the Help system if you are curious).

Problem

Consider the differential equation .  Answer the following questions about it:

a)  Is this equation sensitive to initial conditions or not?  Why/why not?

b)  Use DSolve to plot solutions with initial conditions near (say, from 0.9 to 1.1, perhaps with a stepsize of 0.01).  Does this support or contradict your answer in part a?  (Hint:  Plot all solutions for this problem out to at least .)

c)  Use NDSolve to find an approximate solution with initial condition and plot that on a graph with the solution obtained with DSolve (same initial condition).  Explain why NDSolve gives a different answer (hint: "because it's an approximation" won't get you any credit here).  Is this difference "significant"?  Why/why not?

d)  Increase the WorkingPrecision for NDSolve gradually (remember, the default value starts at 16) until it agrees with DSolve.  Using that WorkingPrecision, increase the domain for t to be twice as long and plot it against DSolve again.  Find a new WorkingPrecision that will fix this one.

e)  Plot NDSolve with WorkingPrecision of 25 (.  What is going on with this?  (It should look qualitatively different than most of the solutions you plotted in part d.)

Created by Mathematica (June 4, 2004)