Differential Equations
From ComputingForScientists
1. Differential Equations
1.1. Objectives
 To introduce the concept of a differential equation and its relationship to a difference equation.
 To introduce how the solution to differential equations can be estimated using the Forward Euler Method.
 To solve basic physics problems that involve differential equations.
1.2. References
In these notes, I assume that the reader has taken differential and integral calculus and that they have not had a course on differential equations. The following two lecture notes also cover this topic at the level covered here:
1.3. Introduction
A common ordinary differential equation that we want to solve is of the form
To solve such an equation means to find a function y(t) that, when diffentiated, gives f(t). In some cases, we can simply guess y(t). An example of an equation for which we can guess the solution is
Because the derivative of − e^{ − t} is e^{ − t}, a solution is
where the constant C has been added because its derivative is zero.
The function y(t) = − e^{ − t} + C is a solution to dy / dt = e^{ − t} in the sense that
To determine a value for the constant C, we would need additional information. Usually the additional information is in the form of a socalled initial condition: the value of y(t) at t = 0.
For certain functions f(t), one can use direct integration to find the exact solution for y(t) if guessing does not work by rewriting dy / dt = f(t) as dy = f(t)dt and integrating both sides
For example, if f(t) is a simple function such as f(t) = e^{ − t}, we need to evaluate
which is
where C_{1} and C_{2} are arbitrary integration constants that can be combined by defining C = C_{2} − C_{1} to give
To verify that this form of y is indeed a solution to dy / dt = e^{ − t}, take the derivative of y(t) = − e^{ − t} + C and see it is f(t) = e^{ − t}.
In some cases, f(t) is complicated and one cannot use guessing or the basic analytical integration techniques learned in an undergraduate integral calculus course to find an exact solution. For some special cases of f(t), methods covered in a differential equations course can be used to obtain an exact solution. However, it is often the case that there is no known method to obtain an exact solution and a numerical method is needed to get an approximate solution. As an example, consider
written in integral form, we have
To find y(t), the integral on the righthand side must be evaluated, but it is not one that has a simple representation that appears in tables of integrals. As a result, we need to resort to finding a numerical approximation for y(t). An obvious candidate method for finding a numerical approximation is numerical integration. However, in this section, we will use the Forward Euler method to solve equations of the form dy / dt = f(t) even though numerical integration could be used. In the #Appendix it is shown that for equations of the form dy / dt = f(t), the Forward Euler method is equivalent to the left rectangle integration method. The motivation for using Forward Euler for this type of problem here is that when we consider equations of the form dy / dt = f(t,y) or problems that involve more than one variable (e.g., dy / dt = f_{y}(t,x,y), dx / dt = f_{x}(t,x,y)), the numerical integration methods considered previously do not apply.
There are many methods for numerically estimating the y(t) that satisfies an equation of the form dy / dt = f(t) (or, more generally, an equation of the form dy / dt = f(t,y)). The easiest numerical method to implement and interpret is named Forward Euler. It is important to note that in practice, this method is rarely used as the approximations it produces are often poor. However, the method is simple to interpret and implement, and for these reasons, it is used to introduce numerical solutions of differential equations.
There is an entire field of research devoted to developing methods to obtain more accurate solutions to differential equations. Many mathematical software programs have builtin functions where one needs to only specify the equation for the derivative, the initial conditions, and then a solution is given. The approach taken here is for you to first understand in detail one simple numerical method used to solve differential equations (Forward Euler). Then you will be asked to use this method to solve a problem and then compare with the solution from a builtin function that someone else wrote (in MATLAB, these functions are ODE23
and ODE45
, for example; they implement the RungeKutta method).
The motivation for this approach is that for some problems you may encounter, the builtin methods may not work, and you will need to write your own numerical method. And when interpreting numerical solutions to differential equations, it is useful to understand how the solutions work, even if you did not implement the program that created the solutions.
In summary,
 Sometimes one can guess an analytic solution to an equation of the form dy / dt = f(t),
 Sometimes direct integration can be used to find an analytic solution to an equation of the form dy / dt = f(t),
 For function of the form dy / dt = f(t), numerical integration can be used to compute a numerical solution for y(t), and
 Here we will use the Forward Euler method to solve equations of the form dy / dt = f(t) and then use it to obtain numerical solutions for equations that cannot be solved using numercial integration.
1.4. dy/dt = f(t) using Forward Euler
As described in the #Introduction, equations of the form dy / dt = f(t) can be solved numerically using numerical integration. However, here we are going to develop an alternative solution method that can be used to solve equations that cannot be solved using numerical integration.
The key assumption that is made in the Forward Euler method for solving differential equations is that we can approximate the continuous differential dy / dt with a discrete differential Δy / Δt.
That is, we assume that we can approximate the derivative (slope) of a curve y(t) at a given point of interest from the slope of line made by connecting a straight line from that point (at t) to another point nearby (at t + Δt). In introductory calculus textbooks, derivatives are introduced by taking the limit of the ratio Δy / Δt as Δt becomes infinitesimally small. Here we are going to use a small (but not infinitesimally small) value for Δt as an approximation. (See the problem #Approximation_of_Derivatives for a problem that asks you to compapre the exact and approximate values of the slope of a function.)
The distinction between a differential equation and a difference equation obtained from approximating a differential equation is that the differential equation involves dt, which is an infinitesimally small increment of time, and a difference equation approximation to a differential equation involves a small, but noninfinitesimal, value of Δt. As expected, we will find that the approximation improves as the value used for Δt decreases (but the amount of computation required to get an approximation will increase).
With the Forward Euler approximation, we can write the differential equation
as a difference equation
or, using the definition ,
and finally
(1)
To use the Forward Euler method, we next need to choose a value of Δt that is small enough to make the assumption "accurate enough". The procedure for choosing an appropriate Δt, and a clarification of the meaning of "accurate enough" is given in the next subsection.
The previous equation can now be used to estimate y(t_{2}) by using t_{1} as the argument to f(t):
where
Now, given an initial value of y(t_{1}), we can use the above equation to provide an estimate of y(t_{2}).
To find and estimate of y(t_{3}), we rewrite the previous equation but add one to each subscript:
In this equation, we use the last computed value of y(t_{2}) to estimate y(t_{3}). This process can be repeated for an arbitrarily long time range to produce an estimate of y(t) at times t_{1}, t_{2}, ....
1.5. dy/dt = f(t) using Forward Euler Example
1.5.1. Hand Calculation
By hand (paper and pencil), estimate y(t) that satisfies dy / dt = at over the time range of t = 0 to 0.3 assuming that y(t = 0) = 1.1 and a = 0.9. To estimate y(t), use Δt = 0.1.
The calculation is
(Given)
(Given)
(By definition)
Note that if we were told to use Δt = 0.05 that we would have needed to do approximately twice as many calculations to obtain an estimate at y(t = 0.3).
1.5.2. Program Calculation
clear a = 0.9; Dt = 0.1; % Δt variable % First array element corresponds to t = 0. t(1) = 0; y(1) = 1.1; % Initial value of y. % Second array element corresponds to t = Dt = 0.1. t(2) = t(1) + Dt; y(2) = y(1) + Dt*a*t(1); % Gives 1.1 % Third array element corresponds to t = 2*Dt = 0.2. t(3) = t(2) + Dt; y(3) = y(2) + Dt*a*t(2); % Gives 1.109 % Fourth array element corresponds to t = 3*Dt = 0.3. t(4) = t(3) + Dt; y(4) = y(3) + Dt*a*t(3); % Gives 1.27 % Display t and y arrays t y
When executed, this program displays
t = 0 0.1000 0.2000 0.3000 y = 1.1000 1.1000 1.1090 1.1270
Once a few steps are written out as a program in longhand form, it should become clear how to use a for
loop generalization.
clear % Forward Euler solution to dy/dt = at from t=0 through t = 0.3 in steps of Δt = 0.1. a = 0.9; Dt = 0.1; % Δt variable % Forward Euler method t(1) = 0; y(1) = 1.1; % Note that last value of y calculated is y(4) due to the i+1 in equation for y below. for i = 1:3 % Stop at i = 3 so that last y value is y(4) t(i+1) = t(i) + Dt; y(i+1) = y(i) + Dt*a*t(i); end t y
1.5.3. Exact Solution
The equation dy / dt = at is an example of an equation that can be solved analytically using direct integration. First, rewrite as an integral
then do the integration
To determine the constant C, use the initial condition y(0) = 1.1:
which gives C = 1.1, so the analytic solution is
Although the exact solution is continuous, in order to plot it, we must choose values of t to plot. In general for exact solutions, one should use just enough values of t so that the plotted line appears continuous when plotted. The code for this is
% Exact solution array calculation. "ex" means exact. a = 0.9; tex = linspace(0,0.3,100); % Create 100 points between 0 and 0.3 inclusive. yex = a*tex.^2/2 + 1.1; % Compute exact solution at those points.
1.5.4. Comparison of Numerical and Exact
% Forward Euler solution to dy/dt = at from t=0 through t = 0.3 % in steps of Δt = 0.1. clear a = 0.9; Dt = 0.1; % Δt % Forward Euler method. % "fe" in variable names indicate that it is a Forward Euler solution tfe(1) = 0; yfe(1) = 1.1; for i = 1:3 % Note that last value of y calculated is y(4) tfe(i+1) = tfe(i) + Dt; yfe(i+1) = yfe(i) + Dt*a*tfe(i); end % Exact solution. % "ex" in variable names indicate that it is an exact solution a = 0.9; tex = linspace(0,0.3,1000); % Create 1000 points between 0 and 0.3 inclusive. yex = a*tex.^2/2 + 1.1; % Compute exact solution at those points. figure(1);clf % Show approximate (Forward Euler) solution as dots. % Increase dot size to 30 so more visible. plot(tfe,yfe,'k.','MarkerSize',30); hold on;grid on; % Plot exact solution as connected lines so that it appears continuous. plot(tex,yex,'b'); xlabel('$t$'); ylabel('$y$ ','Rotation',0,'HorizontalAlignment','Right'); title('$dy/dt = 0.9t; y(t=0) = 1.1; \Delta t = 0.1$'); legend('Euler Method Soln.','Exact Soln.','Location','NorthWest');
1.5.5. Discussion
Consider a graphical representation of the calculations performed for the Forward Euler solution. In computing y(t_{2}) given y(t_{1}), we multiplied the slope at t_{1} () by Δt to give an increment in y.
Because the slope was zero at t_{1} = 0, the first step from t_{1} to t_{2} resulted in no increase in y. A key assumption in the Forward Euler method is that over the time interval from t to t + Δt, the slope dy / dt does not change from what it was at the beginning of the interval at time t (which was zero for t = 0).
In the following image, the exact solution of y(t) to the equation dy / dt = at is shown as a blue line and the Forward Euler method solution is shown as black dots. The dashed red line is the slope at t = 0 that is used in the Forward Euler solution to step from t_{1} = 0.0 to t_{2} = 0.1. The dashed green line is the slope at t = 0.1 used in the Forward Euler solution to step from t_{2} = 0.1 to t_{3} = 0.2.
1.5.6. Convergence Test
When a numerical algorithm is used, one should compute and visualize how the solution changes when a parameter involved in the accuracy of the solution is decreased. For the Forward Euler method, this parameter is Δt. If an exact solution is known, the relative difference between the exact and approximate solution should be considered. There are two possible ways of quantifying the error: (1) at the final time step and (2) at each time step.
In this example, we are going to compute the error at the final time step, at t=0.3
. Before doing this, we need to generalize the program used. In the code excerpt below, the number of steps taken was 3 because we wanted the final timestep calculated to correspond to t=0.3
. If we simply change Dt
to 0.05, the final timestep, corresponding to the y(4)
will correspond to t=0.015
. Therefore, to ensure the last value of y
corresponds to t=0.3
, we need to make the number of iterations in the loop depend on Dt
.
clear a = 0.9; Dt = 0.1; % Δt tfe(1) = 0; yfe(1) = 1.1; for i = 1:3 % Note that last value of y calculated is y(4) tfe(i+1) = tfe(i) + Dt; yfe(i+1) = yfe(i) + Dt*a*tfe(i); end
There are many ways to do this. One way is to calculate tfe
before the for
loop and use its length to determine the number of iterations. Note that the number of iterations in the following code is length(tfe)1
. This is done so that the length of yfe
equals the length of tfe
.
A second change that we make is to compute the exact solution at the same time values that we have computed the Forward Euler solutions. In the following program, the lines used to compute tex
and yex
previously have been commented out, and in the line following them tfe
has been used instead of tex
.
clear a = 0.9; Dt = 0.1; % Δt tfe(1) = 0; yfe(1) = 1.1; tfe = [0:Dt:0.3]; for i = 1:length(tfe)1 % 1 here because of i+1 in yfe below. yfe(i+1) = yfe(i) + Dt*a*tfe(i); end % Exact solution. % "ex" in variable names indicate that it is an exact solution a = 0.9; %tex = linspace(0,0.3,1000); % Create 1000 points between 0 and 0.3 inclusive. %yex = a*tex.^2/2 + 1.1; % Compute exact solution at those points. yex = a*tfe.^2/2 + 1.1; % Compute exact solution at those points. % Compute the absolute difference between final points (which correspond to t = 0.3). Eabs = yfe(end)yex(end) % Compute the relative difference between final points (which correspond to t = 0.3). Erel = abs(yfe(end)yex(end))/yex(end)
The above program will not compute the difference between the Forward Euler and exact solutions at t =0.3 for any choice of Dt
. When testing a typical program, one would run the program for increasing smaller values of Dt
to ensure that the relative error becomes smaller as Dt
is decreased. The following program can be used to plot the dependence of the relative error on Dt
by wrapping the above program in a for
loop that iterates over values of Dt
.
clear; Dtvals = 0.1./10.^[0:5]; % Dt values to try % or %Dtvals = [0.1,0.01,0.001,0.0001,0.00001]; a = 0.9; tfe(1) = 0; yfe(1) = 1.1; for k = 1:length(Dtvals) Dt = Dtvals(k); tfe = [0:Dt:0.3]; for i = 1:length(tfe)1 % 1 here because of i+1 in yfe below. yfe(i+1) = yfe(i) + Dt*a*tfe(i); end yex = a*tfe.^2/2 + 1.1; % Relative difference between final points Erel(k) = abs(yfe(end)yex(end))/yex(end); % Display error in console window fprintf('Dt = %.1e, Error = %.1e\n',Dt,Erel(k)); end figure(1);clf; loglog(Dtvals,Erel,'k.','MarkerSize',30); grid on; title('$dy/dt = 0.9t; y(t=0) = 1.1$'); xlabel('$\Delta t$'); legend('Error = $ \frac{y_{euler}(0.3)y_{exact}(0.3)}{y_{exact}(0.3)} $',... 'Location','NorthWest') ylabel('Error');
1.6. Choosing Δt
When numerically solving any ordinary differential equation, a choice must be made for the timestep parameter Δt so that the numerical solution is at a level of accuracy that is desired. In general, the smaller the Δt, the more accurate the solution. However, smaller values of Δt require more computation so one should choose this parameter with consideration of both accuracy and computation time.
As a general rule, after identifying an initial value of Δt try Δt / 10 and Δt / 100 and compare the solutions.
In the example given for dy / dt = 0.9t, we used Δt = 0.1 initially. This value was intentionally chosen to be large so that the errors were easy to see. If the problem statement is to solve this equation from t = 0.0 through t = 0.3, we would make the following observations:
 The rate of change increases with time;
 The maximum rate of change of dy / dt = at will occur at t = 0.3; and
 If we want the maximum increments in y for the Forward Euler solution to be on the order of, say, 0.01, then our initial value of Δt should be .
There is no single recipe for determining the initial value of Δt to try. For example, for periodic solutions, one should think about the smallest timescales of variation expected from the differential equation and then use 1/10th or 1/100th of that time scale. For example, if a periodic solution is expected with a period of T, choose the initial value of Δt = T / 10.
Note that in a numerical methods course you would learn more rigid rules for choosing Δt and more advanced methods for understanding how the solution accuracy depends on Δt. For each numerical method used to solve a differential equation (e.g., Forward Euler, Heun's method, etc.), there are mathematical results that can be derived that tell you approximately how much error is expected for a given Δt. Many advanced differential equation solvers use adaptive time stepping in which the Δt value is decreased when the estimated value of y is rapidly changing.
In summary, an initial value of Δt should be selected based on reasoning about the problem. The value of Δt used in the program should be smaller than this with the understanding that smaller Δt will require a longer computation time. In problems where an exact solution is not available, one should use several algorithms to compute an approximate solution and use physical reasoning to determine if the numerical solution makes sense.
1.7. Variations on Forward Euler
Another method is the socalled Backward Euler method for equations of the form dy / dt = f(t) uses the slope at the end of the time interval instead of the slope at the start of the time interval.
The Forward Euler algorithm is
The Backward Euler algorithm is
Instead of using only the slope at the beginning or end of the time interval, one could approximate the slope using averages of the slopes at and
1.8. dy/dt = f(y,t) using Forward Euler
The Forward Euler method can also be used to solve the more general problem of finding the y(t) that satisfies
by computing new values of y given past values of y:
Note that the Backward Euler algorithm
has a complication for this type of differential equation: the righthand side requires knowledge of y(t_{2}) which is unknown. In the #Appendix, Heun's method is described, which addresses this complication.
1.9. dy/dt = f(y,t) using ODE23 /ODE45
1.9.1. 1D ODEs
Now that you have had experience with solving an equation of the form dy / dt = f(y,t) using the straightforwardtoimplement Forward Euler method, it is time to use more sophisticated methods that at more difficult to implement, but provide a better estimation of the solution.
In all of the examples given thus far, one needs two pieces of information to estimate y(t): (1) The time rate of change of y(t), and (2) An initial value for y(t).
The functions ODE23
and ODE45
have syntax
[t,y] = ode23(dydt,tspan,yo)
The argument dydt
is a the name of a function that takes inputs of t
and y
and returns dy / dt. This function will be called many times by ODE23
as it develops the solution.
tspan
is a twoelement array with the initial and final times for which the solution is to be computed. Unlike the Forward Euler method described above, the values for the estimate of y(t) provided by ODE23
and ODE45
will not have uniformly spaced time values. If one desires the solution on a uniform time grid, use an array of values for tspan
, e.g., tspan = [0:0.01:10]
.
yo
is simply the value of y(0).
As an example, consider solving
dy / dt = − y
with the initial condition y(0) = 1 over the time range t = 010.
dydt = @(t,y) [y]; % inline function tspan = [0,10]; % Provide solution over this time range yo = 1; % Initial condition [t,y] = ode23(dydt,tspan,yo); % Estimate t and y given dydt and yo. plot(t,y,'.','MarkerSize',10) xlabel('t'); ylabel('y'); title('ODE23 solution to dy/dt = y, y(0) = 1')
In the above example, the output time values for the solution are not uniform. To verify this, enter diff(t)
on the command line. Often it is important to have values of y
at a fixed set of time values. As an example consider comparing the solutions of Forward Euler, ODE23
, and ODE45
. Anytime one uses a program, one should verify that you understand how it works, and that the program works as expected, by studying the program's solution to a problem for which an exact solution is known.
In the following code, solutions to the problem dy / dt = − y; y(0) = 1 plotted individually and then compared at time values of t = 0, 0.1, 0.2, ..., 10.0. To do this, all solutions must have solutions for the same set of time steps.
clear; yo = 1; dt = 0.1; % Will write code so that all solutions provide y values for these t values. t = [0:dt:10]; % Exact solution to dy/dt = y; y(0) = 1. yex = exp(t); % Forward Euler (fe) yfe(1) = yo; for i = 1:length(t)1 yfe(i+1) = yfe(i)  yfe(i)*dt; end % yfe and yex will be row vectors, but ODE23/45 % return column vectors. Transpose them so that % they have the same shape. yfe = yfe'; yex = yex'; dydt = @(t,y) [y]; % inline function [t23,y23] = ode23(dydt,t,yo); [t45,y45] = ode45(dydt,t,yo); % Plot all solutions on linear scale figure(1);clf; plot(t,yex,''); hold on;grid on; plot(t,yfe,'k.','MarkerSize',20) plot(t,y23,'r.','MarkerSize',15) plot(t,y45,'g.','MarkerSize',10) xlabel('t'); ylabel('y','Rotation',0); legend('Exact','Forward Euler','ODE23','ODE45') title('Solutions to dy/dt = y, y(0) = 1') % Solution is exponential, so better graphical % representation is with a semilog yaxis. (Solution % will be linear on this type of axis.) figure(2);clf; semilogy(t,yex,''); hold on;grid on; semilogy(t,yfe,'k.','MarkerSize',20) semilogy(t,y23,'r.','MarkerSize',15) semilogy(t,y45,'g.','MarkerSize',10) xlabel('t'); ylabel('y','Rotation',0); legend('Exact','Forward Euler','ODE23','ODE45') title('Solutions to dy/dt = y, y(0) = 1') % Difference between solutions is difficult to % see on other plots, so plot relative error with respect to exact. figure(3);clf; semilogy(t,abs((yfeyex)./yex),'k.','MarkerSize',15) hold on;grid on; semilogy(t,abs((y23yex)./yex),'r.','MarkerSize',15) semilogy(t,abs((y45yex)./yex),'g.','MarkerSize',15) xlabel('t'); ylabel('y','Rotation',0); legend('Forward Euler  Exact',... 'ODE23  Exact',... 'ODE45  Exact',... 'Location','SouthEast') title('Solutions to dy/dt = y, y(0) = 1')
Note that the parameter dt
is not used in ODE23
and ODE45
. The step size is determined automatically by the algorithms associated with these functions. These functions have optional parameters that indirectly control the step size; these optional parameters may be modified using the ODESET
function to create a structure that is passed to ODE23
and ODE45
as the last arguement. For example,
opts = odeset('RelTol',1e8,'AbsTol',1e8); [t45,y45] = ode45(dydt,t,yo,opts);
The inline function dydt = @(t,y) [y];
had two input arguments: t
and y
. This should always be the case, even if one of the variables is not used. For example, if the ODE to be solved is
dy / dt = 1 − t
will have
dydt = @(t,y) [1t];
and if
dy / dt = 1 − yt
then
dydt = @(t,y) [1y*t];
1.9.2. 2D ODEs
The function ODE23
and ODE45
also work for coupled sets of ordinary differential equations. In this case, the function dydt
must return derivatives for each variable and the initial condition yo
is no longer a scalar  it is an array with as many elements as there are variables. As an example, consider the two differential equations (corresponding to the van der Pol equation):
The first step is to change the names of the variables so that it is easier to interpret them as array elements. Let x = Y_{1} and y = Y_{2}.
In code, the variable Y
will be an array where the first element corresponds to Y_{1} and the second element corresponds to Y_{2}. The inline equation now takes an input of Y
at a given time t
and returns dY_{1} / dt as the first element and dY_{2} / dt as the second element of an array. In the following command, the inline function returns an array with two elements; the first element is dY_{1} / dt and the second is dY_{2} / dt.
dYdt = @(t,Y) [Y(2); (1Y(1)^2)*Y(2)  Y(1)];
In the above, we have used the notation Y
instead of y
to indicate that Y
is an array and not a scalar.
We need two initial conditions, one for x (math>Y_1</math>) and the other for y (math>Y_2</math>)). This is specified using
Yo = [0; 2] % Initial conditions for x and y
With these changes, we can now use ODE23
or ODE45
:
[t,Y] = ode23(dYdt,tspan,Yo);
The output Y
has two columns. The first column corresponds to x (Y_{1}) and the second to y (Y_{2}). The full program is given below:
dYdt = @(t,Y) [Y(2); (1Y(1)^2)*Y(2)  Y(1)]; Yo = [0; 2] % Initial conditions for x and y [t,Y] = ode23(dYdt,tspan,Yo); x = Y(:,1); y = Y(:,2); plot(t,x); hold on; plot(t,y); legend('x','y'); title('ODE23 solution to VDP eqn with x(0) = 0 and y(0) = 2');
1.9.3. Using external function files
In the above examples, the time derivative function was returned by an inline function
dydt = @(t,Y) [y]; % For dy/dt = y dYdt = @(t,Y) [Y(2); (1Y(1)^2)*Y(2)  Y(1)]; % For van der Pol equations.
For simple problems, this approach is preferred. However, it has the limitation that one can not easily pass parameters to it. For example, if we wanted to have mu
as a variable, we might try (this won't work):
mu = 0.9; dYdt = @(t,Y) [Y(2); mu*(1Y(1)^2)*Y(2)  Y(1)]; Yo = [0; 2] % Initial conditions for x and y [t,Y] = ode23(dYdt,tspan,Yo);
To work around this limitation, you can define the function dYdt
in an external file and then pass parameters as shown. For example, if the following is saved in a file named vdp.m
function dYdt = vdp(t,Y,mu) dYdt = [Y(2); mu*(1Y(1)^2)*Y(2)  Y(1)]; % or % dYdt(1) = Y(2); % dYdt(2) = mu*(1Y(1)^2)*Y(2)  Y(1);
then the syntax is now
mu = 1.0; Yo = [0; 2] % Initial conditions for x and y [t,Y] = ode23(@vdp,tspan,Yo,mu);
Note that to pass more than one parameter, make mu
an array in the above code and adjust the function file vdp.m
accordingly.
2. Problems
2.1. Exact Solutions
An example of an ordinary differential equation that can be solved by exactly by direct integration is encountered in kinematics is
Assuming that v(0) = v_{o}, show that the result is
Consider the definition of velocity as the rateofchange of position:
.
Assuming that x(0) = x_{o}, use direct integration to show
2.2. Exact Solutions
Repeat the previous problem assuming that v(t_{f}) = v_{f} and x(t_{f}) = x_{f} (that is, we know the final condition and not the initial condition).
2.3. Exact Solutions
Suppose that
.
Explain the problem with this solution:
giving, after defining C = C_{2} − C_{1},
2.4. Exact Solutions
2.4.1.
Show that the function
satisfies the differential equation
and the initial condition y(0) = 0.
2.4.2.
Guess an expression for y(t) that satisfies the differential equation
and the conditions
 y(0) = 1,
 y(0) = 0, and
 y(10) = 1.
(You should have three expressions, one for each condition. An expression for y(t) can be said to be a solution to a differential equation if it satisfies the differential equation when plugged in and satisfies the initial condition.)
2.4.3.
Show that the equation
and the definition of velocity
can be used to write the equation
2.4.4.
Show that the function
satisfies the differential equation
and the initial conditions x(0) = 0, and dx / dt  _{t = 0} = v(0) = 1.
Using the definition of velocity:
2.4.5.
Find an expression x(t) that satifies the differential equation
and the initial conditions x(0) = 0, and dx / dt  _{t = 0} = 1.
2.4.6.
Show that y(t) = ae^{t / 2} + be^{t / 2} satisfies the differential equation
Find a and b so that the initial conditions y(0) = 2 and dy / dt  _{t = 0} = 1 / 3 are also satisfied.
2.5. Approximation of Derivatives I
A key assumption that is made when solving differential equations is that we can approximate the slope dy / dt of a function y(t) at a given value of t using Δy / Δt.
Consider the function y(t) = t^{2}. Analytically, the exact value of the slope at any value of t is dy / dt = 2t.
Use y(t) = t^{2} and Δt = 0.1 and the definition
to estimate the slope of the curve y(t) = t^{2} at t = 1,2 and 3.
Enter the values in the table below.
Table for Δt = 0.1
t dy/dt Δy / Δt
1
2
3
Then create a table using Δt = 0.01
t dy/dt Δy / Δt
1
2
3
On graph paper, draw y(t), and at t = 1
 a line showing the exact slope,
 a line with the estimated slope computed using Δt = 0.1 and
 a line with the estimated slope computed using Δt = 0.01.
2.6. Approximation of Derivatives II
A key assumption that is made when solving differential equations is that we can approximate the slope dy / dt of a function y(t) at a given value of t using Δy / Δt.
Consider the function y(t) = e^{ − t}. Analytically, the exact value of the slope at any value of t is dy / dt = − e^{ − t}.
Use y(t) = e^{ − t} and Δt = 0.1 and the definition
to estimate the slope of the curve y(t) = e^{ − t} at t = 1,2 and 3.
Enter the values in the table below.
Table for Δt = 0.1
t dy/dt Δy / Δt
1
2
3
Then create a table using Δt = 0.01
t dy/dt Δy / Δt
1
2
3
On graph paper, draw y(t), and at t = 1
 a line showing the exact slope,
 a line with the estimated slope computed using Δt = 0.1 and
 a line with the estimated slope computed using Δt = 0.01.
2.7. Hand Calculation I
By hand, estimate y(t) that satisfies dy / dt = − 0.8t^{2} over the time range of t = 0 to 0.3 assuming that y(t = 0) = 1.1. Compute the estimate of y(t) first using Δt = 0.1 and then using Δt = 0.05.
2.8. Hand Calculation II
By hand, estimate y(t) that satisfies dy / dt = e^{ − t} over the time range of t = 0 to 0.3 assuming that y(t = 0) = 1.1. Compute the estimate of y(t) first using Δt = 0.1 and then using Δt = 0.05.
2.9. Forward Euler vs. Numerical Integration
Can the numerical integration methods covered previously be used to solve all types of differential equations (equations for which the lefthand side is dy / dt)? If yes, explain why we bother with the Forward Euler method. If no, given an example of a differential equation that cannot be solved using numerical integration but can be solved using the Forward Euler method.
2.10. Forward Euler vs. Left Rectangle Method
It was noted that for equations of the form dy / dt = f(t), the Forward Euler method is equivalent to the left rectangle method.
Let f(t) = t^{2} and y(0) = 0.
 Use the left rectangle method from Numerical_Integration with a rectangle width of 0.1 to find y(t = 1).
 Use the Forward Euler method with Δt = 0.1 to find y(t = 1).
2.11. 1D ODE and Forward Euler
Solve the ODE dy / dt = e^{t} numerically using forward Euler subject to the initial condition y(t = 0) = 1 from t = 0 to t = 1.
Compare your solution to the exact solution at t = 1.
Store your values for the numerical solution in an array named yn
and store the values for the time steps in an array named t
. Plot your solution using plot(t,yn)
.
Create an array named code ye
that contains the exact solutions evaluated at the values in the array t
. Compare the numerical and analytical solution using
plot(t,yn,'b','LineWidth',2); hold on; plot(t,ye,'r','LineWidth',2);
Answer 

clear yn(1) = 1; Dt = 0.1; t = [0:Dt:1]; for i = 1:length(t)1 yn(i+1) = yn(i) + Dt*exp(t(i)); end % Exact solution ye = exp(t); plot(t,yn,'b','LineWidth',2); hold on; plot(t,ye,'r','LineWidth',2); legend('Numerical Solution','Exact Solution','location','northwest'); title('dy/dt = e^t'); xlabel('t'); ylabel('y'); Alternative method for numerical solution where yn(1) = 1; Dt = 0.1; t(1) = 0; for i = 1:10 yn(i+1) = yn(i) + Dt*exp(t(i)); t(i+1) = t(i) + Dt; end 
2.12. 1D ODE and Forward Euler
Solve dy / dt = t^{2} using forward Euler with y(t_{1}) = 0, t_{1} = 0, and Δt = 0.1.
Create a plot comparing the solution to the exact solution over the range t = 0 − 1. The plot should include the a title that indicates the difference between the numerical solution and the exact solution at t = 1. If the variable E
contains the error, this can be done using the command
title(sprintf('(y(t=1)y_e(t=1))/y_e(t=1) = %.4g',E))
Create a plot showing the error (y(t = 1) − y_{e}(t = 1)) / y_{e}(t = 1) for Δt = 0.1, 0.01, 0.001, 0.0001, and 0.00001.
Answer 

2.13. 1D ODE and ODE23
Solve the differential equation dy / dt + y = ty with initial condition of y(0) = 0 over the time range of t = 0 through t = 1 using ODE23
.
2.14. 1D ODE and ODE23
Solve the differential equation dy / dt + y = ty with initial condition of y(0) = 2 over the time range of t = 0 through t = 2 using ODE23
.
2.15. 1D ODE and ODE23/ODE45
Compute the solutions of dy / dt + y = f(t) over the time range t = [0,1] with f(t) = 1 and y(0) = 0 using
 Forward Euler

ode23

ode45
then
 Plot the solutions on the same axes along with the analytic solution, y_{a}.
 Compute (y(1) − y_{a}(1)) / y_{a}(1), where y_{a} is the analytic solution.
 Provide an explanation for why the Forward Euler method overestimates the amplitude for all values of t.
2.16. 1D ODE and ODE23
Repeat the previous problem for f(t) = t / 5.
2.17. 2D ODE and Forward Euler
Use the Forward Euler method to solve
dv / dt = − x + sin(2t)
dx / dt = v
with x(0) = 0, and dx / dt  _{t = 0} = 1 with a time step of Dt=0.1
and at time values of t = 0, 0.1, 0.2
. Do this on a piece of paper using only a calculator and show your work. Show your solutions in a handwritten table with columns of t
, x
, and v
. Turn in your solution on a piece of paper.
Answer 

v(i+1) = v(i) + Dt*(x(i) + sin(2*t)); x(i+1) = x(i) + Dt*v(i); x(1) = 0; v(1) = 1; v(2) = v(1) + 0.1*(x(1) + sin(2*0.1)); % 1.0199 x(2) = x(1) + 0.1*v(1); % 0.1000 v(3) = v(2) + 0.1*(x(2) + sin(2*0.2)); % 1.0488 x(3) = x(2) + 0.2*v(2); % 0.3040 
2.18. 2D ODE and ODE23
The secondorder equation of motion for a harmonic oscillator is
Let
then the secondorder equation can be written as two firstorder equations
or
This can be written in vector notation
2.18.1. I
If x(0) = 0 and dx / dt  _{t = 0} = 1, what is the analytic solution to the homogeneous equation d^{2}x / dt^{2} + x = 0?
Plot the exact solution over exactly three cycles and properly label the axes.
Answer to I and II: [3]; images for this problem can be found in [4].
2.18.2. II
Use the Forward Euler method to solve
dv / dt = − x
dx / dt = v
with x(0) = 0, and dx / dt  _{t = 0} = 1. Plot the approximate solution on the same axes as the exact solution and in your code provide a justification for your choice of Dt
.
Your plot should look similar to the following (the formatting of the labels on the xaxis will look different, and you do not need to have the same title).
2.18.3. III
Use the Forward Euler method to solve
dv / dt = − x + sin(2t)
dx / dt = v
with x(0) = 0, and dx / dt  _{t = 0} = 1
Answer to III and IV [5]; images for this problem can be found in [6].
2.18.4. IV
Look for examples of using MATLAB's ODE solver ode23
and then use ode23
to solve II and II. Create a plot the ode23
solution with the exact solution.
Your plot for II should look similar to the following (the formatting of the labels on the xaxis will look different, and you do not need to have the same title).
2.19. Physics Application: Friction I
A box is released with initial velocity of 10 m/s and slides along a horizontal floor with a coefficient of kinetic friction of μ which is time dependent and varies as μ = 0.1(1 − t / τ_{o}), where t is the time in seconds since the box was released and τ_{o} = 1 s.
Use the Forward Euler method to estimate the speed of the box after 0.1 seconds have elapsed. Use a step size of 0.05 seconds.
Estimate how long it takes the box to stop moving. Compare this with the exact time it will take the box to stop moving.
2.20. Physics Application: Friction II
A box is released with initial velocity of 10 m/s and slides along a horizontal floor with a coefficient of kinetic friction of μ which is time dependent and varies as μ = 0.1(1 − tx / (Lτ), where t is the time in seconds since the box was released, x is the distance the object has traveled in meters, τ = 1 s, and L = 2 m.
Use the Forward Euler method to estimate the speed of the box after 0.1 seconds have elapsed. Use a step size of 0.05 seconds.
Estimate how long it takes the box to stop.
2.21. 1D ODE and Forward Euler
The drag force on a 1kg object is , where s is the speed of the object and b = 1 N/(m/s). The object is dropped from a height of 10 meters with an initial speed of 0.0 m/s. Estimate the speed of the object 0.02 seconds after it has been released using Forward Euler with a step size of 0.01 seconds.
Answer 

% Equation of motion is % m*ds/dt = m*g + b*abs(s) % m = 1 and b = 1 s(1) = 0; Dt = 0.01; s(2) = s(1) + Dt*(g+abs(s(1))); s(3) = s(2) + Dt*(g+abs(s(2))); 
2.22. Physics Application: Drag
The drag force on an object is , where s is the speed of the object and b = 1 N/(m/s)^{2}. The object is dropped from a height of 10 meters with an initial speed of 0.0 m/s. Estimate the speed of the object 0.02 seconds after it has been released using Forward Euler with a step size of 0.01 seconds.
2.23. Physics Application: Drag
Compute the solution of dv / dt = kv^{2} − g with initial condition of v(0) = 0 using ODE23
and Forward Euler. Create a plot showing the Forward Euler solution and the ODE23
solution. The plot should have appropriate legend and axes labels.
2.24. Physics Application: Newton's Law of Cooling
Newton's law of cooling states that the time rate of change of the temperature T of an object is proportional to the difference between the temperature of the object and the ambient temperature, T_{a}. Assume a proportionality constant of 1/(10^{7} seconds) when answering the following questions.
 A rock in space has a temperature of 300 K. Assume the ambient temperature is 2 K. How long will it take for the rock to cool to 200 K, 3 K, and 2 K? You are encouraged to check your answer using an analytic solution.
 Suppose the ambient temperature happens to increase at a rate of 1 K/year. How long will it take for the rock to cool to 200 K, 3 K, and 2 K?
Answer these questions using Forward Euler and/or ODE23
. When executed, the following information should be displayed with the ?
replaced with calculated values. Use the fprintf
function to create the display; see Writing_Text_Files.
Constant ambient temperature using ??  Time to cool to 200 K: ? years Time to cool to 3 K: ? years Time to cool to 2 K: ? years Increasing ambient temperature using ??  Time to cool to 200 K: ? years Time to cool to 3 K: ? years Time to cool to 2 K: ? years
2.25. Physics Application: The Heat Equation
The equation that describes how the temperature of an object, T
changes with time (Newton's Law of Cooling) is
 Use a spreadsheet to determine how the temperature of an object that starts out at 100 degrees C cools when the ambient temperature,
Ta
is 0. Assume that q = 0.1 and Δt = 0.1. At approximately what time does the temperature equal onehalf of its initial temperature?  Use a spreadsheet to determine if or how the time it takes for the temperature to fall by onehalf depends on the initial temperature of the object. Describe in words you how you used the spreadsheet to answer this question and your conclusion.
2.26. Physics Application: Free Fall
2.26.1.
In a constant gravitational field, the equations of motion for a particle that does not experience drag is
Determine the exact (analytic) amount of time that the particle will take to hit Earth's surface using y(t = 0) = 2R_{E} above Earth's surface and an initial velocity v(t = 0) = 0. Store the value in a variable named tcga
.
Use Forward Euler to solve for both y and v using y(t = 0) = 2R_{E} above Earth's surface, an initial velocity v(t = 0) = 0, and a step size of Dt = 10
seconds and determine the last time value for which y(t) > 0. Store this time value in a variable named tcg
. Use R_{E} = 6,371 km.
2.26.2.
In reality, the gravitational force that an object experiences in free fall is not constant but rather increases as the object moves towards Earth. In this case, the equation of motion for no drag is
where y is the distance of the object relative to the surface of Earth.
Suppose an object is dropped from a height of y(t = 0) = 2R_{E} above Earth's surface with an initial velocity v(t = 0) = 0.
Determine the exact (analytic) amount of time that the particle will take to hit Earth's surface.
Using a step size of Dt = 10
seconds, determine the last time value for which y(t) > 0. Store the value in a variable named tvg
.
Use the following statements to print your answer to the screen
fprintf('Time to reach surface with constant grav. force (analytic): %.2e [s]\n',tcga) fprintf('Time to reach surface with constant grav. force (euler) : %.2e [s]\n',tcg) fprintf('Time to reach surface with variable grav. force (euler) : %.2e [s]\n',tvg)
On the same axes (use hold on
), and with correct labels, plot the position versus time for
 constant gravitational force using the Forward Euler method
 variable gravitational force using the Forward Euler method
2.27. Physics Application: Nonideal spring
A spring has a restoring force F = kx, where k = k_{o}(1 + x^{2} / L^{2}), k_{o} = 1N / m. and L = 10 meters.
A horizontally traveling 1 kg particle with a speed of 2 m/s strikes this spring.
 Use a graphical method to estimate the maximum amount of spring compression.
 Solve the equation of motion for the object while it is on the spring and use the solution to find the amount of spring compression.
 If the particle was traveling in the vertical direction downward and hit the spring, will the amount of compression differ? If so, compute the amount of compression using a graphical method.
Answer 

The energy equation for an ideal spring is
where PE = < math > k / 2x^{2} in general and PE_{o} = 0 for this problem because spring is initially uncompressed. At maximum compression, v_{f} = 0, so KE_{f} = 0. This leaves
When k is constant, the analytical solution is
If k is not constant, the potential energy of a spring is not simply k / 2x^{2}. Note that the energy equation used previously comes from
where W is the work done by any forces on the object as it moves along its path.
Using gives
Doing the integral and assuming x_{o} = 0 gives
So the energy equation
can be written as
This is a 4thorder polynomial. One could use the function
and plot y vs. x_{f}. The x_{f} value where the curve crosses zero is the solution. To solve this, first use the bruteforce method using a constant k and compare with the analytical solution. Then ask what is expected to happen physically when k increases with x. The spring gets stiffer, so we expect the nonconstant k case to result in a smaller x_{f}. The following creates the plot and the zoom tool was used to find . (The zoom tool is available when the figure window is undocked.) This is smaller than that computed for k = constant, which was x_{f} = 2.0. Additional code is given for finding the solution using other methods. L = 10; ko = 1; vo = 2; m = 1; xf = [0:0.0001:3]; y1 = ko*xf.^2  m*vo^2; y2 = ko*(1 + xf.^2/(2*L^2)).*xf.^2  m*vo^2; figure(1);clf plot(xf,y1,'r','LineWidth',2) hold on;grid on; plot(xf,y2,'b','LineWidth',2) xlabel('x_f') ylabel('\Delta KE  W') legend('k=k_o','k=k_o(1+x^2/L^2)') % Additional code below. % To solve this using <code>roots</code>, simplify equation to % <math>y=x_f^4/200 + x_f^2  4</math>. r = roots([1/(2*L^2),0,1,0,4]); fprintf('Estimate using roots: %.5f\n',real(r(end))); % Another method of obtaining <math>x_f</math> is to search for % the location where the <code>y</code> array transitions from below to above zero: % Find index where crossing. I = find(y2(1:end1) < 0 & y2(2:end) >= 0); %y2(I) % Should be below zero %y2(I+1) % Should be above or equal to zero %xf(I) % Value when y2 below zero %xf(I+1) % Value when y2 above or equal to zero % Estimate as average value before and after crossing zero xf_estimate = (xf(I) + xf(I+1))/2; fprintf('Estimate using energy equation: %.5f\n',xf_estimate); clear; dydt2 = @(t,y) [y(2);(1+y(1)^2/10^2)*y(1)]; opts = odeset('RelTol',1e9,'AbsTol',1e9); [t,y] = ode45(dydt2,[0,3],[0,2],opts); clf; %plot(y(:,1),y(:,2)); %grid on; %xlabel('x'); %ylabel('v'); I = find(y(1:end1,2) > 0 & y(2:end,2) <= 0); %y(I,2) % Should be below zero %y(I+1,2) % Should be above or equal to zero %y(I,1) % Value when y2 below zero %y(I+1,1) % Value when y2 above or equal to zero % Estimate as average value before and after crossing zero xf_estimate = (y(I,1) + y(I+1,1))/2; fprintf('Estimate using ODE45: %.5f\n',xf_estimate); 
2.28. Physics Application: Kinematics
A projectile is launched vertically from the surface of Earth with a velocity of 10 m/s.
 Compute and plot the analytical solution for its position and velocity as a function of time.
 Compute and plot a numerical solution for its position and velocity as a function of time.
Answer 

clear; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Analytic (note that these are not rate of change equations) % y = 10*t  (1/2)*9.8*t^2; % v = 10  9.8*t % Create arrays for y and v for analytic solution ta = [0:0.1:1.0]; for i = 1:length(ta) ya(i) = 10*ta(i)(1/2)*9.8*ta(i)*ta(i); va(i) = 10  9.8*ta(i); end % Note  there is a way of doing the above without a loop: % ya = 10*t  (1/2)*9.8*t.^2; % va = 10  9.8*t; % The "." before the ^ means to take the power of each element in the % array. figure(1) clf; plot(ta,va,'b','LineWidth',2); xlabel('t (s)'); ylabel('v (m/s)'); title('Analytic solution'); grid on; figure(2) clf; plot(ta,ya); xlabel('t (s)'); ylabel('y (m)'); title('Analytic solution'); grid on; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Numerical % v(t = 0) = 10;% Initial velocity % y(t = 0) = 0; % Initial height (surface of earth) % dv/dt = 9.8; % RateofChange equation for velocity % dy/dt = 109.8*t % RateofChange equation for position % Discrete version of above equations: % v(i+1) = v(i) + dt*(9.8); % y(i+1) = y(i) + dt*(109.8*t(i)) vn(1) = 10; % Initial velocity yn(1) = 0; % Initial height tn = [0:0.01:1.0]; dt = tn(2)  tn(1); for i = 1:length(tn)1 vn(i+1) = vn(i) + dt*(9.8); yn(i+1) = yn(i) + dt*(109.8*tn(i)); end figure(3); plot(tn,vn,'r','LineWidth',2); xlabel('t (s)'); ylabel('v (m/s)'); title('Numerical solution'); grid on; figure(4); plot(tn,yn,'r','LineWidth',2); xlabel('t (s)'); ylabel('y (m/s)'); title('Numerical solution'); grid on; 
2.29. Physics Application: Kinematics
The time rate of change of the position of an object is equal to − gt, where g = 9.8 m/s^{2}:
Assume that the initial position of the object is 10 m.
 On a piece of paper, estimate the position of the object at t = 0.2 s assuming Δt = 0.1 seconds using the Forward Euler method.
 Write a program that uses the Forward Euler method to estimate y(t) for t = 0 − 10 s assuming Δt = 0.1 seconds.
 Plot position computed using the Forward Euler method versus time for t = 0 − 10 s.
 On the same plot, plot position versus time for the analytic solution for t = 0 − 10 s.
For future reference, experiment with smaller and larger values of Δt = 0.1 and determine how the difference between the solutions depends on Δt = 0.1.
Your solutions should be almost identical for small values of t. For larger values of t, the solutions will diverge.
Make sure to label your axes and include a legend.
See [7] for how to plot two lines on the same figure.
Answer 

% Paper estimate % y(1) = 0; % t(1) = 0; % y(2) = y(1) + dt*(g*t(1)) % y(2) = 10 + 0.1*(9.8*t(1)) => y(2) = 10 % t(2) = 0.1 % y(3) = y(2) + 0.1*(9.8*t(2)) % y(3) = 10 + 0.1*(9.8*0.1) => 9.9020 clear % Forward Euler Solution g = 9.8; ye(1) = 10; t(1) = 0; dt = 0.1; t = [0:dt:10dt]; for i = 1:length(t) ye(i+1) = ye(i)  g*t(i)*dt; t(i+1) = t(i) + dt; end % Show Forward Euler at steps computed by hand. ye(2) ye(3) % Analytic solution % v = vo + at ya = 10  (g/2)*t.^2; figure(1);clf % Use thicker line for Euler so that Analytic line does not completely % cover. plot(t,ye,'k','LineWidth',1); hold on;grid on; plot(t,ya,'b','LineWidth',1); plot(t,yeya,'g','LineWidth',1); legend('Forward Euler','Analytic','Euler  Analytic','Location','SouthWest'); ylabel('y (m/s)'); xlabel('t (s)'); % Display difference between Euler and Analytic ye(end)ya(end) % dt = 0.1 gives 4.94 % dt = 0.01 gives 0.4905 
2.30. Physics Application: Kinematics
The time rate of change of the velocity of an object is constant and equal to − g, where g = 9.8 m/s^{2}:
Assume that the initial velocity of the object is 0 m/s.
 On a piece of paper, estimate the velocity of the object at t = 0.2 s assuming Δt = 0.1 seconds using the Forward Euler method.
 Write a program that uses the Forward Euler method to estimate v(t) for t = 0 − 10 s assuming Δt = 0.1 seconds. Verify that the first few velocity values computed match your paper estimate.
 Plot velocity computed using the Forward Euler method versus time for t = 0 − 10 s.
 On the same plot, plot velocity versus time for the analytic solution for t = 0 − 10 s.
Make sure to label your axes and include a legend.
See [8] for how to plot two lines on the same figure.
(In this special case, Forward Euler will give the same solution as the analytic solution for any choice of Δt = 0.1.)
Save your answer in a file named HW5b.m
.
Answer 

% Paper estimate % v(1) = 0; % v(2) = v(1) + dt*(g) % v(2) = 0 + 0.1*(9.8) => v(2) = .98 % v(3) = v(2) + 0.1*(9.8) % v(3) = 0.98 + 0.1*(9.8) => 1.96 % Note that for this particular problem, Forward Euler matches % Analytic exactly. clear % Forward Euler Solution g = 9.8; ve(1) = 0; t(1) = 0; dt = 0.1; for i = 1:100 ve(i+1) = ve(i)  g*dt; t(i+1) = t(i) + dt; end % Show Forward Euler at steps computed by hand. ve(2) ve(3) % Analytic solution % v = vo + at va = g.*t; figure(1);clf % Use thicker line for Euler so that Analytic line does not completely % cover. plot(t,ve,'k','LineWidth',3); hold on;grid on; plot(t,va,'g','LineWidth',1); legend('Forward Euler','Analytic'); ylabel('v (m/s)'); xlabel('t (s)'); 
2.31. Physics Application: Kinematics
A small box slides down an L = 11 meter plane inclined 30 degrees from horizontal.
 Draw the freebody diagram to write down an equation for dv / dt, the rateofchange of velocity along the direction of the plane. Clearly label your axes and the value for which position is zero. Verify that the equation makes sense for a horizontal plane and a plane that is nearly vertical.
Use the Forward Euler method to estimate the time that it takes to reach the bottom of the plane for the following cases:
 The coefficient of kinetic friction is
0.1
. Save the time in a variable namedt1
.  The coefficient of kinetic friction increases with time according to 0.1(1 + t / T), where T = 10 seconds. Save the time in a variable named
t2
. Indicate in a comment in your code if this value is expected to be greater, equal, or less thant1
and provide a physical justification for your answer.  The coefficient of kinetic friction decreases with time according to 0.1(1 − t / T), where T = 10 seconds. Save the time in a variable named
t3
. Indicate in a comment in your code if this value is expected to be greater, equal, or less thant1
and provide a physical justification for your answer.  The coefficient of kinetic friction increases with the distance, s, the object has slid along the plane according to . Save the time in a variable named
t4
. Indicate in a comment in your code if this value is expected to be greater, equal, or less thant1
and provide a physical justification for your answer.
When your program is executed, the four values of time corresponding to the four cases should be displayed (in order). To do this, enter the following four lines at the end of your program:
fprintf('t1 = %.2f\n',t1); fprintf('t2 = %.2f\n',t2); fprintf('t3 = %.2f\n',t3); fprintf('t4 = %.2f\n',t4);
Answer 

Basic solution given below. A better solution would use a while loop to run while Note that if clear v(1) = 0; x(1) = 0; dt = 0.01; mu = 0.1; g = 9.8; t(1) = 0; for i = 1:1000 v(i+1) = v(i) + dt*g*(sind(30)mu*cosd(30)); x(i+1) = x(i) + dt*v(i); t(i+1) = t(i) + dt; end I = find(x<=11); t1 = t(I(end)); T = 10; for i = 1:1000 v(i+1) = v(i) + dt*g*(sind(30)mu*(1+t(i)/T)*cosd(30)); x(i+1) = x(i) + dt*v(i); t(i+1) = t(i) + dt; end I = find(x<=11); t2 = t(I(end)); T = 10; for i = 1:1000 v(i+1) = v(i) + dt*g*(sind(30)mu*(1t(i)/T)*cosd(30)); x(i+1) = x(i) + dt*v(i); t(i+1) = t(i) + dt; end I = find(x<=11); t3 = t(I(end)); T = 10; for i = 1:1000 v(i+1) = v(i) + dt*g*(sind(30)mu*(1sqrt(x(i)/10))*cosd(30)); x(i+1) = x(i) + dt*v(i); t(i+1) = t(i) + dt; end I = find(x<=11); t4 = t(I(end)); fprintf('t1 = %.2f\n',t1); fprintf('t2 = %.2f\n',t2); fprintf('t3 = %.2f\n',t3); fprintf('t4 = %.2f\n',t4); 
2.32. Physics Application: Kinematics
A 1 kg mass is dropped from a height of 10 meters at which point it strikes a spring with spring constant of 10 N/m.
2.32.1. Analytical
In this part, use techniques covered in your firstsemester physics class.
 On a piece of paper, sketch the expected position versus time from the time the mass is dropped until after it hits the spring and returns back to its maximum height.
 Write an equation for the position and velocity versus time while the particle is moving downward and before it hits the spring.
 Write an equation for the position and velocity versus time while the particle is on the spring.
 Plot the position versus time from the time the mass is dropped until it hits the spring.
2.32.2. Numerical
In this part, use techniques covered in this course [9] to develop a numerical (and approximate) solution for the position and velocity of the particle.
 Write down an equation for the acceleration of the particle before it hits the spring. Use the Forward Euler method to compute the velocity of the particle as a function of time while it is falling and before it hits the spring.
 Write down an equation for the acceleration of the particle while it is on the spring. Use the Forward Euler method to compute the velocity of the particle as a function of time while it is on the spring.
 On the same plot, show the Analytic and Numerical position versus time curves for one cycle. Label your plots.
 On the same plot, show the Analytic and Numerical velocity versus time curves for one cycle. Label your plots.
Note that the analytic equation for the position of the mass while it is on the spring (y < 0) is
where , t is measured from the time the mass hits the spring, and v_{o} is the velocity when the mass hits the spring.
Answer: https://github.com/rweigel/computingforscientists/blob/master/ODEs/FreeFallSpring.m
2.33. Physics Application: Mass on Spring
In the configuration shown, the mass is in equilibrium at position x_{1o}. x_{1} and v_{1} are the position and velocity of the spring, respectively.
2.33.1. Part I
 Write an equation that relates dv_{1} / dt to one or more of the variables m_{1}, k, x_{1o}, x_{1}, and v_{1}.
 Write an equation that relates dx_{1} / dt to one or more of the variables m_{1}, k, x_{1}, and v_{1}.
 If m_{1} is given a swift kick at t = 0 so that v_{1}(0) = v_{o} and x_{1}(0) = x_{1o}, sketch the expected v_{1}(t), x_{1}(t). Label important features on the sketch in terms of the parameters m_{1}, k, v_{o}, and x_{1o}.
 Sketch the expected KE(t), PE(t), and E(t) = PE(t) + KE(t).
2.33.2. Part II
 Rewrite the equation for dv_{1} / dt in a form that can be solved using the Forward Euler method.
 Rewrite the equation for dx_{1} / dt in a form that can be solved using the Forward Euler method.
 What is a reasonable value to use for Δt in the Forward Euler method (in terms of one or more of the parameters m_{1}, k, v_{o}, and x_{1o})?
2.33.3. Part III
The mass of 10^{6} kg is given a swift kick to the right at a velocity v_{1o} = 1 m/s. The equilibrium position x_{1o} is 5000 m. Both springs have a spring constant of 1 N/m.
 Write a program that uses the Forward Euler method to compute x_{1}(t) and v_{1}(t). When your program is executed, two figure windows should appear showing x_{1}(t) and v_{1}(t), respectively. Your plots should have appropriate labels with units.
 In your program where you define the time step, include a statement for how or why you used the value.
2.33.4. Part IV
MATLAB has a buildin function named ODE23
that can be used to give a more accurate numerical solution.
Use the MATLAB documentation and online resources to figure out how to use ODE23
to compute the position and velocity of the mass considered.
When I run your program, I should see two plots
 A with two lines corresponding to position versus time using Forward Euler and
ODE23
 A with two lines corresponding to velocity versus time using Forward Euler and
ODE23
Your plots must have appropriate axis labels, titles, and legends.
2.33.5. Part V
Repeat all of the steps in Part IIV except use v_{1}(0) = v_{o} and v_{2}(0) = − v_{o}. Both masses are initially at their equilibrium position.
2.34. Physics Application: Pendulum
The image above shows a freebodydiagram for a point mass at the end of a string that has zero mass. This type of problem is typically covered in introductory physics textbooks.
By definition, the rateofchange of θ is ω (equivalent to the rateofchange of position being velocity). So we can write
The torque equation associated with the free body diagram is
2.34.1. Part I
 If the pendulum is released at at an angle of θ_{o} with ω_{o} = 0 at t = 0, sketch the expected curves for θ(t) and ω(t). Label important features on the sketch in terms of one or more of the parameters m, L, g, and θ_{o}.
 Sketch the expected KE(t), PE(t), and E(t) = PE(t) + KE(t).
Answer 

To draw a sketch, note that at t = 0, θ = θ_{o},ω = 0 and after a period T the pendulum will return to this state: t = T, θ = θ_{o},ω = 0 Another point of interest is when the pendulum is vertical: at t = T / 4, θ = 0,ω = − ω_{max} and when the pendulum is at its maximum height on the other side of the vertical: t = T / 2 and θ = − θ_{o},ω = 0 Finally, at t = 3T / 4, θ = 0,ω = ω_{max} These points should be plotted along with a smooth curve that connects them to obtain sketches of θ(t) and ω(t). At this point, we don't know ω_{max} or T. To determine ω_{max}, use energy conservation. From an introductory physics textbook, the potential energy of is mgy and the kinetic energy is mL^{2}ω^{2} / 2. The variable y is related to θ via
so that the (constant) total energy is E = mgL(1 − cos(θ)) + mL^{2}ω^{2} / 2 using t = 0,θ = θ_{o},ω = 0 gives E = mgL(1 − cos(θ_{o})). Using t = T / 4,θ = 0,ω = − ω_{max} gives
or
The above information can be used to plot KE(t), PE(t), and E(t). Note that there is no way of obtaining T at this point unless we invoke the approximation
In this case, giving . From an introductory physics textbook, for small θ, θ(t) varies sinusoidally with an angular frequency of . To determine a shift of the sin wave, we can use the initial conditions. We are given that initially θ = θ_{o} with ω_{o} = 0. The equation
satisfies the condition θ(0) = θ_{o}. That it satisfies ω_{o} = 0 can be shown by using the fact that ω = dθ / dt:
which, for t = 0 gives ω = 0. Also from an introductory physics textbook, the potential energy of is mgy and the kinetic energy is mL^{2}ω^{2} / 2. The variable y is related to θ via
so that y = L(1 − cos(θ)), which in the smallangle approximation is Lθ^{2} / 2. Thus the total energy in the smallangle approximation is is
or, using θ(t) = θ_{o}sin(ω_{o}t) and ω(t) = θ_{o}ω_{o}cos(ω_{o}t) is
or, using
or
So that the total energy is constant and
To sketch these, plot points of time for which ω_{o}t = 0, π / 2, π, 3π / 2, 2π, 5π / 2, etc., and then connect the points with a smooth curve. 
2.34.2. Part II
 Rewrite the equation for dθ / dt in a form that can be solved using the Forward Euler method.
 Rewrite the equation for dω / dt in a form that can be solved using the Forward Euler method.
 What is a reasonable initial value to try for Δt in the Forward Euler method (in terms of one or more of the parameters (m, L, g, and θ_{o})?
2.34.3. Part III
Suppose θ_{o} = 30^{o} and L = 10^{ − 3} m.
 Write a program that uses the Forward Euler method to compute θ(t) and ω(t). When your program is executed, two figure windows should appear showing θ(t) and ω(t), respectively. Your plots should have appropriate labels with units.
 In your program where you define the time step, include a statement for how or why you used the value.
2.34.4. Part IV
MATLAB has a buildin function named ODE23
that can be used to give a more accurate numerical solution.
Use the MATLAB documentation and online resources to figure out how to use ODE23
to compute the position and velocity of the mass considered.
When I run your program, I should see two plots
 A plot with two lines corresponding to θ versus time using Forward Euler and
ODE23
 A plot with two lines corresponding to ω versus time using Forward Euler and
ODE23
Your plots must have appropriate axis labels, titles, and legends.
2.34.5. Part V
 Plot ω versus θ using the
ODE23
solution using ω_{o} = 0 with values θ_{o} = 10^{o}, 20^{o}, ..., 170^{o} and 179^{o}. Each line is called a phase space trajectory. In your plot, each phase space trajectory should have a different color. Be prepared to answer the following questions Why are the line lengths different?
 The initial conditions selected always had ω = 0. Could you have selected θ_{o} = 0 and a nonzero value of ω_{o} and obtained trajectories that overlapped the ones shown? If yes, what are the values of ω_{o} that you would use?
 How is a phase space trajectory related to a streamline?
2.35. Direct Integration
Solve dy / dt = t^{2} using numerical integration (any method, e.g., rectangle left/right/middle, trapezoid, Simpson) with y(t_{1}) = 0, t_{1} = 0, and Δt = 0.1.
Create a plot of the solution to the exact solution over the range t = 0 − 1. The plot should include the a title that indicates the difference between the numerical solution and the exact solution at t = 1.
Create a plot showing the error (y(t = 1) − y_{e}(t = 1)) / y_{e}(t = 1) for Δt = 0.1, 0.01, 0.001, 0.0001, and 0.00001.
Answer 

2.36. Heun's Method
The equation for the velocity of a spherical mass near Earth's surface is mdv / dt = − g. The effect of air resistance can be modeled as a force F = − CρAv^{2}, where ρ is the density of air, A is the crosssectional area of the sphere, and C is a constant.
Write a program that solves dv / dt = − 9.8 + v^{2} / L using v(0) = 0 m/s, Δt = 0.1 s, and L = 10 m.
Create a sketch (on a piece of paper) comparing the expected solution to dv / dt = − 9.8 + v^{2} / L with the exact solution of dv / dt = − g over the range t = 0 − 10 s.
Create a plot comparing the numerical solution to dv / dt = − 9.8 + v^{2} / L with the exact solution of dv / dt = − g over the range t = 0 − 3 s.
Answer 

3. Appendix
3.1. dy/dt = f(t) using numerical integration
As as alternative to the Forward Euler method, we note that
can be solved by direct integration by writing
and evaluating the integral on the lefthand side with limits of y(t_{N}) and y(t_{1}) and approximating the integral on the righthand side using numerical integration
or
where the continuous integral on the righthand side has been approximated using the "left" rectangle method and and the answer will be the same as that found by the above method. A sketch of the rectangles associated with this method will show that the estimate of the area will be an underestimate and the first rectangle has an area of zero. In the above equation, to find y
at t_N
= t_2
, t_3
, etc., one needs to evaluate the righthand side once for each t_N
.
As an example, consider f(t) = at. A straightforward (but slow  don't use it) implementation of this equation for a=1
is
% y(N) = y(1) + sum of Dt*f(t) from 1 to N1 % For each N, sum Dt*f(t) for t = 1 to N1. y(1) = 0; Dt = 0.1; t = [0:Dt:1]; for i = 1:length(t)1 s = 0; for j = 1:i s = s + Dt*t(j) end y(i+1) = y(1) + s; end
A problem with this implementation is the number of computations that must be performed due to the nested for
loops. The number of steps that must be taken is much larger than that for the Forward Euler, so we seek an alternative method by writing out a few steps in longhand.
y(2) = y(1) + Dt*f(1) y(3) = y(1) + Dt*f(1) + Dt*f(2) y(4) = y(1) + Dt*f(1) + Dt*f(2) + Dt*f(3) ... y(N) = y(1) + Dt*f(1) + ... + Dt*f(N1)
This can be implemented using
y(1) = 0; Dt = 0.1; t = [0:Dt:1]; s = 0; for i = 1:length(t)1 s = s + Dt*t(i); y(i+1) = y(1) + s; end
To save an addition operation at each step, we can write
s = y(1); for i = 1:length(t)1 s = s + Dt*t(i); y(i+1) = s; end
To improve the accuracy, we can reduce or seek an alternative algorithm that will has a higher accuracy give the same . Given the limitation noted above, that the slope at the start of the integral was assumed to be the slope over the entire time span of we can use the average slope between and :
and this method will give the same result as using the numerical integration approach when the Trapezoid method is used instead of the rectangle method.
3.2. dy/dt = f(y,t) using Heun's method
The methods described thus far were applied to differential equations of the form
More generally, we will want to find a numerical solution to
We can use Euler's method, for which the slope over the interval to was approximated using the slope at :
but the method in which the slope over the interval is approximated using averages of the slopes at and
has a complication  the righthandside of the equation contains the unknown . One possibility for addressing this complication is to rewrite the above equation as
and use guesses for to compute the lefthandside until it equals the righthandside.
Another possibility is called Heun's method. In Heun's method, in the term
is approximated by using Euler's method
giving
or
for which the righthandside is completely determined. Note that if , Heun's method is equivalent to the numerical integration method using trapezoids.
Example
Solve with .
To implement this in MATLAB, one could define f
using an anonymous function to make the code look more like the mathematical notation:
y(1) = 0; % Initial condition. Dt = 0.1; % Time step t = [0:Dt:3]; % Time grid f = @(t,y) 1 + y*t; % Euler estimate of y(2) yeu(2) = y(1) + Dt*f(t(1),y(1)); % Heun estimate of y(2) y(2) = y(1) + (Dt/2)*( f(t(1),y(1)) + f(t(2),yeu(2)) );
The more traditional notation is (spaces added for readability)
y(1) = 0; % Initial condition. Dt = 0.1; % Time step t = [0:Dt:3]; % Time grid yeu = y(1) + Dt*(1 + y(1)*t(1)); % Euler estimate of y(2). k1 = 1 + y(1)*t(1); % f(t, y) k2 = 1 + yeu *t(2); % f(t+Dt, yeu) y(2) = y(1) + (Dt/2)*(k1 + k2); % Heun estimate of y(2).
or, if we don't want to have a time array, the first step is
Dt = 0.1 t = 0; y(1) = 0; yeu = y(1) + Dt*(1 + y(1)); k1 = 1 + y(1)*t; k2 = 1 + yeu*(t+Dt); y(2) = y(1) + (Dt/2)*(k1 + k2);