MATH 302-010
Ordinary Differential Equations
Prof. D. A. Edwards
Copyright 2024

Numerics and Homogeneous ODEs (Revised)

We begin with our standard command:

Quit[]

Euler’s Method

IMPORTANT: In this notebook, h is used for ∆t!

As previously discussed, when DSolve doesn’t return an analytic expression for the solution of an ODE, we have to use numerical methods:

eq1 = {y'[t]==Cos[50*Sin[t]],y[0]==0}
DSolve[eq1,y[t],t]

The most basic numerical method of all is Euler’s method.  We can tell Mathematica to use Euler’s method by adding the option Method->”ExplicitEuler”  after the normal arguments of NDSolve.  The stepsize is given using the option StartingStepSize->h value.

eulsolve = NDSolve[eq1,y,{t,0,2*Pi},Method->"ExplicitEuler",StartingStepSize->1/10];
Plot[Evaluate[y[t]/.eulsolve],{t,0,2*Pi}]

The jagged edges in the graph should give us pause; it seems as if something isn’t right.  Euler’s method is rather crude and has large errors, so let’s reduce the step size.  We again use a table, but to change variable values by amounts other than 1, we have to use the syntax {variable,low value,high value, step between values}.

eultab= Table[NDSolve[eq1,y,{t,0,2*Pi},Method->"ExplicitEuler",StartingStepSize->1/N],{N,10,100,10}];
Plot[Evaluate[y[t]/.eultab],{t,0,2*Pi}]

Note that decreasing the step size gets rid of the jagged edges.  In particular the value with h=1/100 is very close to the results using standard NDSolve.

eul100= NDSolve[eq1,y,{t,0,2*Pi},Method->"ExplicitEuler",StartingStepSize->1/100];
sol1 = NDSolve[eq1,y,{t,0,2*Pi}];
Plot[Evaluate[y[t]/.eul100],{t,0,2*Pi}]
Plot[Evaluate[y[t]/.sol1],{t,0,2*Pi}]

When using NDSolve, Mathematica automatically makes adjustments to handle areas when the function changes rapidly.

Exercises

1.  Consider the differential equation y'=(e^(10sin(30t))) with the initial condition y(0)=0.

a)  Use Euler’s method with h=1/10 and compare with the default solution from NDSolve by looking at plots of both.

Solution

ex1 = {y'[t]==Exp[10*Sin[30*t]],y[0]==0}
eul100= NDSolve[ex1,y,{t,0,2*Pi},Method->"ExplicitEuler",StartingStepSize->1/10];
sol1 = NDSolve[ex1,y,{t,0,2*Pi}];
Plot[Evaluate[y[t]/.eul100],{t,0,2*Pi}]
Plot[Evaluate[y[t]/.sol1],{t,0,2*Pi}]

The step size is obviously too large for Euler’s method, as we are missing the stair-step pattern.

b)  Vary h in Euler’s method from 1/10 to 1/100.  For what values of h does it begin to resemble the correct solution?

Solution

extab= Table[NDSolve[ex1,y,{t,0,2*Pi},Method->"ExplicitEuler",StartingStepSize->1/N],{N,10,100,10}];
Table[Plot[Evaluate[y[t]/.extab[[j]]],{t,0,2*Pi}],{j,1,10}]

The 4th graph and beyond look good, so we have that an upper bound for h should be h=1/40.

Second Order

When solving second-order ODEs, the easiest way to write the equation using the prime notation:

sec1=y''[t]+3*y'[t]+y[t]==0

When solving these equations, we usually have two unknown constants:

DSolve[sec1,y[t],t]

Mathematica can also handle the case where only one condition is given, in which case you have only one constant:

ic1=y[0]==1
DSolve[{sec1,ic1},y[t],t]

We also specify a derivative initial condition using the prime notation:

ic2=y'[0]==0
DSolve[{sec1,ic1,ic2},y[t],t]

The same ideas hold true for numerical solutions as well.

numsec = NDSolve[{sec1,ic1,ic2},y,{t,0,3}];
y[3]/.numsec
y'[3]/.numsec

Note that we can get both y and its derivative using NDSolve, which can be useful for certain purposes.  Plotting also works the same way:

Plot[y[t]/.numsec,{t,0,3}]
Plot[y'[t]/.numsec,{t,0,3}]

However, we must be more careful asking Mathematica to solve a nonlinear second-order differential equation, even numerically.  First, we try to solve analytically.  If it is computing for a long time, hit CMD+. to abort the calculation.

sec2=(y''[t])^2-Sin[y[t]]==0
DSolve[sec2,y[t],t]

Next we try numerically, which takes a short while.

numsec2=NDSolve[{sec2,ic1,ic2},y,{t,0,3}]

Note that there are two solutions (corresponding to the plus or minus square root of y'').  Plotting them, we have

Plot[y[t]/.numsec2[[1]],{t,0,3}]
Plot[y[t]/.numsec2[[2]],{t,0,3}]

Note that neither solution can go all the way to 3, because the second derivative becomes imaginary at the endpoints shown.

Exercise

1.  Solve the differential equation y''+3y'-7y=1, both in general and with the initial condition y(0)=3, y'(0)=4.

Solution

eq1=y''[t]+3*y'[t]-7*y[t]==1
ic1=y[0]==3
ic2 = y'[0]==7

DSolve[eq1,y[t],t]
DSolve[{eq1,ic1,ic2},y[t],t]

Text Formatting

In order to make your Mathematica notebooks more appealing, it may be necessary to format the text.  The most common formatting options may be accessed from the toolbar or the  Writing Assistant, which can be viewed by selecting Palettes...Writing Assistant.  There you will find a B button (for bold), an I button (for italics), and a U button (for underlined).  (The keyboard shortcuts for these commands are the same as commonly used for other applications.)  These styles may also be combined.

There are also buttons for aligning a paragraph on the left,

centering it,

or aligning it right.  (Note that it is easiest to have different paragraphs in different cells for this type of formatting.)

You may also change the font or size of the text from the Writing Assistant.