Understanding NDSolve

AP1601
Columbia University

Connecting the Dots

Mathematica can easily make lists of numbers.

sineList = Table[{x, Sin[x]}, {x, 0, 6 π, 6 π/100}];

p1=ListPlot[sineList,PlotStyle->PointSize[0.02`]]

While a list of numbers can be plotted and displayed, it is not a function until we tell Mathematica how to connect the dots.

?Interpolation

?InterpolatingFunction

s= Interpolation[sineList //N]

The "standard form" of displaying an InterpolatingFunction[...] suppresses the long list of numbers used to create the interpolation. We can see the hidden information by using FullForm[...]:

s // FullForm

We must use ListPlot[...] to display a list of numbers. For functions, we must use Plot[...].

p2=Plot[s[x],{x,0,6 π}]

Show[p1,p2]

If you look carefully, the interpolation is slightly curved.
By default, Interpolation[...] connects the dots with third order polynomials.

Show[p1,p2,PlotRange->{{4.`,5.5`},{-1.`,-0.6`}},PlotLabel->"Close Up"]

We can examine the error between the true function  and the interpolation…

Plot[s[x]-Sin[x],{x,0,3 π},PlotPoints->200,PlotLabel->"Interpolation Error"]

Connecting Fewer Dots

When there are fewer points, we also have less information to represent the function.

sineList2 = Table[{x, Sin[x]}, {x, 0, 6 π, 6 π/25}];

blue = RGBColor[0,0,1];

p3=ListPlot[sineList2,PlotStyle->{PointSize[0.02`],blue}]

s2= Interpolation[sineList2 //N]

p4=Plot[s2[x],{x,0,6 π},PlotStyle->{blue}]

Show[p3,p4]

Show[p3,p4,PlotRange->{{4.`,5.5`},{-1.`,-0.6`}},PlotLabel->"Close Up"]

Show[p1,p2,p3,p4,PlotRange->{{4.`,5.5`},{-1.`,-0.6`}},PlotLabel->"Close Up"]

Reducing the Interpolation Order

If we lower the order of the polynomial used to construct the interpolation, then "smoothness" of the interpolation is reduced.

s2= Interpolation[sineList2 //N,InterpolationOrder -> 2]

p4=Plot[s2[x],{x,0,6 π},PlotStyle->{blue}]

Show[p1,p2,p3,p4,PlotRange->{{4.`,5.5`},{-1.`,-0.6`}},PlotLabel->"Close Up"]

InterpolatingFunctions are "Pretty Good" Functions

Mathematica treats interpolating functions as real functions. What I mean by this phrase is that Mathematica can integrate and differentiate interpolatingfunctions in the same way that it differentiates and integrates other functions.

Taking the derivative, returns:

{D[Sin[x],x], D[s[x],x]}

Plot[Evaluate[{(∂_x)Sin[x],(∂_x)s[x]}],{x,0,6 π},PlotStyle->{{},{blue}}]

Show[%,PlotRange->{{2.5`,4.`},{-1.`,-0.6`}},PlotLabel->"Close Up"]

Plot[Evaluate[{(∂_({x,2}))Sin[x],(∂_({x,2}))s[x]}],{x,0,6 π},PlotStyle->{{},{blue}}]

Show[%,PlotRange->{{4.`,5.5`},{0.6`,1.`}},PlotLabel->"Close Up"]

Plot[Evaluate[{(∂_({x,3}))Sin[x],(∂_({x,3}))s[x]}],{x,0,6 π},PlotStyle->{{},{blue}}]

Show[%,PlotRange->{{2.5`,4.`},{0.6`,1.`}},PlotLabel->"Close Up"]

You can control NDSolve

NDSolve is a powerful command with many options that you can control. Most of these control the number of points used to represent your solution and the accuracy of your solution. More accurate solutions require more points. So do solutions having a longer range. Of course, your computer requires more memory and more time to represent more accurate solutions defined on longer time ranges.

Equations for an Ellipse

eqel = {D[x[t],t] == y[t], D[y[t],t] == - b x[t]}

initial = {x[0] == x0, y[0] == y0}

eqs = eqel ~ Join ~ initial

sol1 = NDSolve[eqs /. x0 -> 1 /. y0 -> 0 /. b-> 2, {x[t], y[t]}, {t, 0, 10}]

Plot[Evaluate[{x[t],y[t]}/.sol1],{t,0,10},PlotStyle->{{},{blue}}]

ParametricPlot[Evaluate[{x[t],y[t]}/.sol1],{t,0,10}]

Show[%,AspectRatio->Automatic]

x[t] /. sol1 // FullForm

Options to NDSolve

??NDSolve

sol1 = NDSolve[eqs /. x0 -> 1 /. y0 -> 0 /. b-> 2, {x[t], y[t]}, {t, 0, 1000}]

sol1 = NDSolve[eqs /. x0 -> 1 /. y0 -> 0 /. b-> 2, {x[t], y[t]}, {t, 0, 1000}, MaxSteps -> 15000]

Plot[Evaluate[{x[t],y[t]}/.sol1],{t,0,500},PlotStyle->{{},{blue}}]

Show[%,PlotRange->{{100,120},All}]

Plot[Evaluate[{x[t],y[t]}/.sol1],{t,100,120},PlotStyle->{{},{blue}}]

sol2 = NDSolve[eqs /. x0 -> 1 /. y0 -> 0 /. b-> 2, {x[t], y[t]}, {t, 0, 1000}, PrecisionGoal-> 3, AccuracyGoal->3]

Plot[Evaluate[(x[t]/.sol1)-(x[t]/.sol2)],{t,0,500},PlotLabel->"Difference between More & Less Accurate"]

ParametricPlot[Evaluate[{{x[t],y[t]}/.sol1,{x[t],y[t]}/.sol2}],{t,0,500},PlotStyle->{{},{blue}},AspectRatio->Automatic,Compiled->False]

You can calculate with solutions to NDSolve

Remember, for most calculations that you want to perform with your solutions to the equations of motion, you can treat the solutions from NDSolve just like any other "real" function.

Plot[Evaluate[{x[t]/.sol1,(∂_t)(x[t]/.sol1)}],{t,0,50},PlotStyle->{{},{blue}}]

Plot[Evaluate[{y[t]/.sol1,(∂_t)(x[t]/.sol1)}],{t,0,50},PlotStyle->{{},{blue}}]

radius[t_] = Sqrt[x[t]^2 + y[t]^2] /. sol1

Plot[Evaluate[{radius[t],(∂_t)radius[t]}],{t,0,50},PlotStyle->{{},{blue}}]

Summary

NDSolve[…] is a powerful built-in function in Mathematica for solving differential equations. The solutions returned by NDSolve[…] are interpolating functions. Interpolating functions are "real functions". They can be combined in complex expressions, integrated and differentiated. However, they are only approximations to the actual solution to the differential equations. These approximations generally have increasing error for long-time integrations.