Differential Equations

From ComputingForScientists

Jump to: navigation, search

Contents

  1. Differential Equations
    1. Objectives
    2. References
    3. Introduction
    4. dy/dt = f(t) using Forward Euler
    5. dy/dt = f(t) using Forward Euler Example
      1. Hand Calculation
      2. Program Calculation
      3. Exact Solution
      4. Comparison of Numerical and Exact
      5. Discussion
      6. Convergence Test
    6. Choosing Δt
    7. Variations on Forward Euler
    8. dy/dt = f(y,t) using Forward Euler
    9. dy/dt = f(y,t) using ODE23 /ODE45
      1. 1-D ODEs
      2. 2-D ODEs
      3. Using external function files
  2. Problems
    1. Exact Solutions
    2. Exact Solutions
    3. Exact Solutions
    4. Exact Solutions
    5. Approximation of Derivatives I
    6. Approximation of Derivatives II
    7. Hand Calculation I
    8. Hand Calculation II
    9. Forward Euler vs. Numerical Integration
    10. Forward Euler vs. Left Rectangle Method
    11. 1-D ODE and Forward Euler
    12. 1-D ODE and Forward Euler
    13. 1-D ODE and ODE23
    14. 1-D ODE and ODE23
    15. 1-D ODE and ODE23/ODE45
    16. 1-D ODE and ODE23
    17. 2-D ODE and Forward Euler
    18. 2-D ODE and ODE23
      1. I
      2. II
      3. III
      4. IV
    19. Physics Application: Friction I
    20. Physics Application: Friction II
    21. 1-D ODE and Forward Euler
    22. Physics Application: Drag
    23. Physics Application: Drag
    24. Physics Application: Newton's Law of Cooling
    25. Physics Application: The Heat Equation
    26. Physics Application: Free Fall
    27. Physics Application: Non-ideal spring
    28. Physics Application: Kinematics
    29. Physics Application: Kinematics
    30. Physics Application: Kinematics
    31. Physics Application: Kinematics
    32. Physics Application: Kinematics
      1. Analytical
      2. Numerical
    33. Physics Application: Mass on Spring
      1. Part I
      2. Part II
      3. Part III
      4. Part IV
      5. Part V
    34. Physics Application: Pendulum
      1. Part I
      2. Part II
      3. Part III
      4. Part IV
      5. Part V
    35. Direct Integration
    36. Heun's Method
  3. Appendix
    1. dy/dt = f(t) using numerical integration
    2. dy/dt = f(y,t) using Heun's method

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

\frac{dy}{dt} = f(t)

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

\frac{dy}{dt} = e^{-t}

Because the derivative of et is et, a solution is

\frac{}{}y(t) = -e^{-t} + C

where the constant C has been added because its derivative is zero.

The function y(t) = − et + C is a solution to dy / dt = et in the sense that

\frac{d(-e^{-t} + C)}{dt} = e^{-t}

To determine a value for the constant C, we would need additional information. Usually the additional information is in the form of a so-called 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 re-writing dy / dt = f(t) as dy = f(t)dt and integrating both sides

\int dy = \int f(t) dt

For example, if f(t) is a simple function such as f(t) = et, we need to evaluate

\int dy = \int e^{-t} dt

which is

\frac{}{}y + C_1 = -e^{-t}+C_2

where C1 and C2 are arbitrary integration constants that can be combined by defining C = C2C1 to give

\frac{}{}y = -e^{-t}+C

To verify that this form of y is indeed a solution to dy / dt = et, take the derivative of y(t) = − et + C and see it is f(t) = et.

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

\frac{dy}{dt} =-\sin(t^2)

written in integral form, we have

\int dy = \int \sin(t^2)dt

To find y(t), the integral on the right-hand 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 = fy(t,x,y), dx / dt = fx(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 built-in 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 built-in function that someone else wrote (in MATLAB, these functions are ODE23 and ODE45, for example; they implement the Runge-Kutta method).

The motivation for this approach is that for some problems you may encounter, the built-in 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.

\frac{dy}{dt} \simeq \frac{\Delta y}{\Delta t} \equiv \frac{y(t+\Delta t)-y(t)}{\Delta 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 non-infinitesimal, 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

\frac{dy}{dt} = f(t)

as a difference equation

\frac{\Delta y}{\Delta t} = f(t)

or, using the definition \Delta y/\Delta t \equiv y(t+\Delta t)-y(t)/\Delta t,

\frac{}{}y(t+\Delta t)-y(t) = \Delta t f(t)

and finally

(1)   \frac{}{}y(t+\Delta t) = y(t) + \Delta t f(t)

To use the Forward Euler method, we next need to choose a value of Δt that is small enough to make the assumption dy/dt \simeq \Delta y/\Delta t "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(t2) by using t1 as the argument to f(t):

\frac{}{}y(t_2) = y(t_1) + \Delta t f(t_1)

where \frac{}{}t_2 = t_1 + \Delta t

Now, given an initial value of y(t1), we can use the above equation to provide an estimate of y(t2).

To find and estimate of y(t3), we re-write the previous equation but add one to each subscript:

\frac{}{}y(t_3) = y(t_2) + \Delta t f(t_2)

In this equation, we use the last computed value of y(t2) to estimate y(t3). This process can be repeated for an arbitrarily long time range to produce an estimate of y(t) at times t1, t2, ....

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

\frac{}{}t_1 = 0 (Given)

\frac{}{}y(t_1) = 1.1 (Given)

\frac{}{}t_2 = t_1 + \Delta t = 0 + 0.1 (By definition)

\frac{}{}y(t_2) = y(t_1) + \Delta t f(t_1) = 1.1 + 0.1\cdot (at_1) = 1.1 + 0.1\cdot (0.9\cdot 0.0) = 1.1

\frac{}{}t_3 = t_2 + \Delta t = 0.1 + 0.1 = 0.2

\frac{}{}y(t_3) = y(t_2) + \Delta t f(t_2) = 1.1 + 0.1\cdot (at_2) = 1.1 + 0.1\cdot (0.9\cdot 0.1) = 1.109

\frac{}{}t_4 = t_3 + \Delta t = 0.2 + 0.1 = 0.3

\frac{}{}y(t_4) = y(t_3) + \Delta t f(t_3) = 1.082 + 0.1\cdot (at_3) = 1.109 + 0.1\cdot (0.9\cdot 0.2) = 1.27

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 long-hand 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, re-write as an integral

\int dy = \int atdt

then do the integration

\frac{}{}y = at^2/2 + C

To determine the constant C, use the initial condition y(0) = 1.1:

\frac{}{}1.1 = a0^2/2 + C

which gives C = 1.1, so the analytic solution is

\frac{}{}y = at^2/2 + 1.1

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');
From raw.githubusercontent.com on May 18 2019 22:25:27.

1.5.5. Discussion

Consider a graphical representation of the calculations performed for the Forward Euler solution. In computing y(t2) given y(t1), we multiplied the slope at t1 (dy/dt|_{t=t_1}) by Δt to give an increment in y.

Because the slope was zero at t1 = 0, the first step from t1 to t2 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 t1 = 0.0 to t2 = 0.1. The dashed green line is the slope at t = 0.1 used in the Forward Euler solution to step from t2 = 0.1 to t3 = 0.2.

From raw.githubusercontent.com on May 18 2019 22:25:27.

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');
From raw.githubusercontent.com on May 18 2019 22:25:27.

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 at/0.01 
 = 0.01/0.9\cdot 0.3 \simeq 0.04.

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 so-called 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

\frac{}{}y(t+\Delta t) = y(t) + \Delta t f(t)

The Backward Euler algorithm is

\frac{}{}y(t+\Delta t) = y(t) + \Delta t f(t+\Delta t)

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 \frac{}{}t_1and \frac{}{}t_2

\frac{}{}y(t_2) = y(t_1) + \Delta t \frac{f(t_1) + f(t_2)}{2}

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

\frac{dy}{dt} = f(y(t),t)

by computing new values of y given past values of y:

\frac{}{}y(t_2) = y(t_1) + \Delta t f(t_1,y(t_1))

Note that the Backward Euler algorithm

\frac{}{}y(t_2) = y(t_1) + \Delta t f(t_1,y(t_2))

has a complication for this type of differential equation: the right-hand side requires knowledge of y(t2) 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. 1-D ODEs

Now that you have had experience with solving an equation of the form dy / dt = f(y,t) using the straightforward-to-implement 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 two-element 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 = 0-10.

dydt = @(t,y) [-y]; % in-line 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]; % in-line 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 y-axis. (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((yfe-yex)./yex),'k.','MarkerSize',15)
    hold on;grid on;
    semilogy(t,abs((y23-yex)./yex),'r.','MarkerSize',15)
    semilogy(t,abs((y45-yex)./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',1e-8,'AbsTol',1e-8);
[t45,y45] = ode45(dydt,t,yo,opts);

The in-line 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) [1-t];

and if

dy / dt = 1 − yt

then

dydt = @(t,y) [1-y*t];

1.9.2. 2-D 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):

\frac{dx}{dt} = -y

\frac{dy}{dt} = (1-x^2)y - x

The first step is to change the names of the variables so that it is easier to interpret them as array elements. Let x = Y1 and y = Y2.

\frac{dY_1}{dt} = -Y_2

\frac{dY_2}{dt} = (1-Y_1^2)Y_2 - Y_1

In code, the variable Y will be an array where the first element corresponds to Y1 and the second element corresponds to Y2. The in-line equation now takes an input of Y at a given time t and returns dY1 / dt as the first element and dY2 / 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 dY1 / dt and the second is dY2 / dt.

dYdt = @(t,Y) [-Y(2); (1-Y(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 (Y1) and the second to y (Y2). The full program is given below:

dYdt = @(t,Y) [-Y(2); (1-Y(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); (1-Y(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*(1-Y(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*(1-Y(1)^2)*Y(2) - Y(1)];
  % or
  % dYdt(1) = -Y(2);
   % dYdt(2) = mu*(1-Y(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

\frac{dv}{dt} = -g

Assuming that v(0) = vo, show that the result is

\frac{}{}v(t) = v_o - gt

Consider the definition of velocity as the rate-of-change of position:

\frac{dx}{dt} = v(t).

Assuming that x(0) = xo, use direct integration to show

\frac{}{}x(t) = x_o + v_ot - gt^2/2

2.2. Exact Solutions

Repeat the previous problem assuming that v(tf) = vf and x(tf) = xf (that is, we know the final condition and not the initial condition).

2.3. Exact Solutions

Suppose that

\frac{}{}dy/dt = e^{ty}.

Explain the problem with this solution:

\int dy = y + C_1

\int e^{ty}dt = \frac{e^{ty}}{y} + C_2

giving, after defining C = C2C1,

y = \frac{e^{ty}}{y}+C

2.4. Exact Solutions

2.4.1.

Show that the function

\frac{}{}y(t) = 1-e^{-t/\tau}

satisfies the differential equation

\frac{}{}dy/dt + y = 1

and the initial condition y(0) = 0.

2.4.2.

Guess an expression for y(t) that satisfies the differential equation

\frac{}{}dy/dt + y = 1

and the conditions

  1. y(0) = 1,
  2. y(0) = 0, and
  3. 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

\frac{}{}dv/dt + x = f(t)

and the definition of velocity

v = \frac{}{}dx/dt

can be used to write the equation

\frac{}{}d^2x/dt^2 + x = f(t)

2.4.4.

Show that the function

x(t) = \frac{5}{3}\mbox{sin}(t)-\frac{1}{3}\mbox{sin}(2t)

satisfies the differential equation

\frac{}{}dv/dt = -x + \mbox{sin}(2t)

and the initial conditions x(0) = 0, and dx / dt | t = 0 = v(0) = 1.

Using the definition of velocity:

\frac{}{}v = dx/dt

2.4.5.

Find an expression x(t) that satifies the differential equation

\frac{}{}d^2x/dt^2 + x = \mbox{sin(2t)}

and the initial conditions x(0) = 0, and dx / dt | t = 0 = 1.

2.4.6.

Show that y(t) = aet / 2 + bet / 2 satisfies the differential equation

\frac{}{}d^2y/dt^2 - dy/dt + y/4 = 0

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) = t2. Analytically, the exact value of the slope at any value of t is dy / dt = 2t.

Use y(t) = t2 and Δt = 0.1 and the definition

\frac{\Delta y}{\Delta t} = \frac{y(t+\Delta t)-y(t)}{\Delta t}

to estimate the slope of the curve y(t) = t2 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

  1. a line showing the exact slope,
  2. a line with the estimated slope computed using Δt = 0.1 and
  3. 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) = et. Analytically, the exact value of the slope at any value of t is dy / dt = − et.

Use y(t) = et and Δt = 0.1 and the definition

\frac{\Delta y}{\Delta t} = \frac{y(t+\Delta t)-y(t)}{\Delta t}

to estimate the slope of the curve y(t) = et 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

  1. a line showing the exact slope,
  2. a line with the estimated slope computed using Δt = 0.1 and
  3. 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.8t2 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 = et 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 left-hand 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) = t2 and y(0) = 0.

  1. Use the left rectangle method from Numerical_Integration with a rectangle width of 0.1 to find y(t = 1).
  2. Use the Forward Euler method with Δt = 0.1 to find y(t = 1).

2.11. 1-D ODE and Forward Euler

Solve the ODE dy / dt = et 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);

2.12. 1-D ODE and Forward Euler

Solve dy / dt = t2 using forward Euler with y(t1) = 0, t1 = 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) − ye(t = 1)) / ye(t = 1) for Δt = 0.1, 0.01, 0.001, 0.0001, and 0.00001.

2.13. 1-D 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. 1-D 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. 1-D 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

  1. Forward Euler
  2. ode23
  3. ode45

then

  • Plot the solutions on the same axes along with the analytic solution, ya.
  • Compute (y(1) − ya(1)) / ya(1), where ya is the analytic solution.
  • Provide an explanation for why the Forward Euler method overestimates the amplitude for all values of t.

Answer: Code | Images

2.16. 1-D ODE and ODE23

Repeat the previous problem for f(t) = t / 5.

2.17. 2-D 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 hand-written table with columns of t, x, and v. Turn in your solution on a piece of paper.

2.18. 2-D ODE and ODE23

The second-order equation of motion for a harmonic oscillator is

\frac{}{}d^2x/dt^2 + x = f(t)

Let

\frac{}{}v = dx/dt

then the second-order equation can be written as two first-order equations

\frac{}{}dv/dt + x = f(t)

\frac{}{}dx/dt = v

or

\frac{}{}dx/dt = v

\frac{}{}dv/dt = -x + f(t)


This can be written in vector notation

d\mathbf{Y}/dt = [dx/dt,dv/dt]^T = [v, -x + f(t)]^T

2.18.1. I

If x(0) = 0 and dx / dt | t = 0 = 1, what is the analytic solution to the homogeneous equation d2x / dt2 + 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 x-axis will look different, and you do not need to have the same title).

From raw.githubusercontent.com on May 18 2019 22:25:27.

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 ode23to 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 x-axis will look different, and you do not need to have the same title).

From raw.githubusercontent.com on May 18 2019 22:25:27.

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. 1-D ODE and Forward Euler

The drag force on a 1-kg object is b\cdot |s|, 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.

2.22. Physics Application: Drag

The drag force on an object is b\cdot s^2, 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 = kv2g 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.

Answer: Code | Figure

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, Ta. Assume a proportionality constant of 1/(107 seconds) when answering the following questions.

  1. 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.
  2. 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


\frac{\Delta T}{\Delta t} = -q(T-T_a)

  1. 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 one-half of its initial temperature?
  2. Use a spreadsheet to determine if or how the time it takes for the temperature to fall by one-half 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

\frac{}{}dy/dt = v

\frac{}{}dv/dt = -g

Determine the exact (analytic) amount of time that the particle will take to hit Earth's surface using y(t = 0) = 2RE 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) = 2RE 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 RE = 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

\frac{}{}dy/dt = v

\frac{}{}dv/dt = -gR_E^2/(y+R_E)^2

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) = 2RE 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

  1. constant gravitational force using the Forward Euler method
  2. variable gravitational force using the Forward Euler method

Answer: Code Figure

2.27. Physics Application: Non-ideal spring

A spring has a restoring force F = kx, where k = ko(1 + x2 / L2), ko = 1N / m. and L = 10 meters.

A horizontally traveling 1 kg particle with a speed of 2 m/s strikes this spring.

  1. Use a graphical method to estimate the maximum amount of spring compression.
  2. 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.
  3. 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.

2.28. Physics Application: Kinematics

A projectile is launched vertically from the surface of Earth with a velocity of 10 m/s.

  1. Compute and plot the analytical solution for its position and velocity as a function of time.
  2. Compute and plot a numerical solution for its position and velocity as a function of time.

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/s2:

\frac{dy}{dt} = -gt

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.


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/s2:

\frac{dv}{dt} = -g

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.

2.31. Physics Application: Kinematics

A small box slides down an L = 11 meter plane inclined 30 degrees from horizontal.

  • Draw the free-body diagram to write down an equation for dv / dt, the rate-of-change 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 named t1.
  • 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 than t1 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 than t1 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 0.1(1+\sqrt{s}/\sqrt{L}). 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 than t1 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);

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 first-semester physics class.

  1. 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.
  2. Write an equation for the position and velocity versus time while the particle is moving downward and before it hits the spring.
  3. Write an equation for the position and velocity versus time while the particle is on the spring.
  4. 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.

  1. 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.
  2. 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.
  3. On the same plot, show the Analytic and Numerical position versus time curves for one cycle. Label your plots.
  4. 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

\frac{}{}y(t) = -(v_o/\omega)\sin(\omega t)+(mg/k)\cos(\omega t) - mg/k

where \omega = \sqrt{k/m}, t is measured from the time the mass hits the spring, and vo 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

From github.com on January 20 2018 08:24:54.

In the configuration shown, the mass is in equilibrium at position x1o. x1 and v1 are the position and velocity of the spring, respectively.

2.33.1. Part I

  • Write an equation that relates dv1 / dt to one or more of the variables m1, k, x1o, x1, and v1.
  • Write an equation that relates dx1 / dt to one or more of the variables m1, k, x1, and v1.
  • If m1 is given a swift kick at t = 0 so that v1(0) = vo and x1(0) = x1o, sketch the expected v1(t), x1(t). Label important features on the sketch in terms of the parameters m1, k, vo, and x1o.
  • Sketch the expected KE(t), PE(t), and E(t) = PE(t) + KE(t).

2.33.2. Part II

  • Re-write the equation for dv1 / dt in a form that can be solved using the Forward Euler method.
  • Re-write the equation for dx1 / 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 m1, k, vo, and x1o)?

2.33.3. Part III

The mass of 106 kg is given a swift kick to the right at a velocity v1o = 1 m/s. The equilibrium position x1o is 5000 m. Both springs have a spring constant of 1 N/m.

  • Write a program that uses the Forward Euler method to compute x1(t) and v1(t). When your program is executed, two figure windows should appear showing x1(t) and v1(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 build-in 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

  1. A with two lines corresponding to position versus time using Forward Euler and ODE23
  2. A with two lines corresponding to velocity versus time using Forward Euler and ODE23

Your plots must have appropriate axis labels, titles, and legends.

Answers: Sketches | Code

2.33.5. Part V

Repeat all of the steps in Part I-IV except use v1(0) = vo and v2(0) = − vo. Both masses are initially at their equilibrium position.

From github.com on January 20 2018 08:24:54.

2.34. Physics Application: Pendulum

From github.com on January 19 2018 22:18:46.

The image above shows a free-body-diagram 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 rate-of-change of θ is ω (equivalent to the rate-of-change of position being velocity). So we can write

\frac{d\theta}{dt} = \omega

The torque equation associated with the free body diagram is

mL^2\frac{d\omega}{dt} = -mgL\sin(\theta)

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).

2.34.2. Part II

  • Re-write the equation for dθ / dt in a form that can be solved using the Forward Euler method.
  • Re-write 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 = 30o 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 build-in 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

  1. A plot with two lines corresponding to θ versus time using Forward Euler and ODE23
  2. 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

  1. Plot ω versus θ using the ODE23 solution using ωo = 0 with values θo = 10o, 20o, ..., 170o and 179o. 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
    1. Why are the line lengths different?
    2. The initial conditions selected always had ω = 0. Could you have selected θo = 0 and a non-zero value of ωo and obtained trajectories that overlapped the ones shown? If yes, what are the values of ωo that you would use?
    3. How is a phase space trajectory related to a streamline?

2.35. Direct Integration

Solve dy / dt = t2 using numerical integration (any method, e.g., rectangle left/right/middle, trapezoid, Simpson) with y(t1) = 0, t1 = 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) − ye(t = 1)) / ye(t = 1) for Δt = 0.1, 0.01, 0.001, 0.0001, and 0.00001.

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ρAv2, where ρ is the density of air, A is the cross-sectional area of the sphere, and C is a constant.

Write a program that solves dv / dt = − 9.8 + v2 / 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 + v2 / 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 + v2 / L with the exact solution of dv / dt = − g over the range t = 0 − 3 s.

3. Appendix

3.1. dy/dt = f(t) using numerical integration

As as alternative to the Forward Euler method, we note that

\frac{}{}dy = f(t)dt

can be solved by direct integration by writing

\int dy = \int f(t) dt

and evaluating the integral on the left-hand side with limits of y(tN) and y(t1) and approximating the integral on the right-hand side using numerical integration

\frac{}{}y(t_N)-y(t_1) = \int_{t_1}^{t_N} f(t) dt \simeq \sum_{t=t_1}^{t=t_{N-1}} f(t) \Delta t

or

\frac{}{}y(t_N) \simeq y(t_1) + \sum_{t=t_1}^{t=t_{N-1}} f(t) \Delta t

where the continuous integral on the right-hand side has been approximated using the "left" rectangle method and \frac{}{}\Delta t = (t_N-t_1)/N 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 right-hand side once for each t_N.

As an example, consider f(t) = at. A straight-forward (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 N-1
% For each N, sum Dt*f(t) for t = 1 to N-1.
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 long-hand.

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(N-1)

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 \frac{}{}\Delta t or seek an alternative algorithm that will has a higher accuracy give the same \frac{}{}\Delta t. 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 \frac{}{}\Delta t we can use the average slope between \frac{}{}t_1and \frac{}{}t_2:

\frac{}{}y(t_2) = y(t_1) + \Delta t a (t_1 + t_2)/2

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

\frac{}{}dy/dt = f(t)

More generally, we will want to find a numerical solution to

\frac{}{}dy/dt = f\bigl(t,y(t)\bigr)

We can use Euler's method, for which the slope over the interval \frac{}{}t_1 to \frac{}{}t_2 was approximated using the slope at \frac{}{}t_1:

\frac{}{}y(t_2) = y(t_1) + \Delta t f(t_1,y(t_1))

but the method in which the slope over the interval is approximated using averages of the slopes at \frac{}{}t_1and \frac{}{}t_2

\frac{}{}y(t_2) = y(t_1) + \Delta t \frac{f\bigl(t_1,y(t_1)\bigr) + f\bigl(t_2,y(t_2)\bigr)}{2}

has a complication - the right-hand-side of the equation contains the unknown \frac{}{}y(t_2). One possibility for addressing this complication is to rewrite the above equation as

\frac{}{}y(t_2) - f\bigl(t_2,y(t_2)\bigr)/2 = y(t_1) + \Delta t f\bigl(t_1,y(t_1)\bigr)/2

and use guesses for \frac{}{}y(t_2) to compute the left-hand-side until it equals the right-hand-side.

Another possibility is called Heun's method. In Heun's method, \frac{}{}y(t_2) in the term \frac{}{}\frac{}{}f\bigl(t_2,y(t_2)\bigr)

is approximated by using Euler's method

\frac{}{}y_E(t_2) = y(t_1) + \Delta t f\bigl(t_1,y(t_1)\bigr)

giving

\frac{}{}y(t_2) = y(t_1) + \frac{\Delta t}{2} f\bigl(t_1,y(t_1)\bigr) + f\bigl(t_2,y_E(t_2)\bigr)

or

\frac{}{}y(t_2) = y(t_1) + \frac{\Delta t}{2} \biggl[f\bigl(t_1,y(t_1)\bigr) + f\bigl(t_2,y(t_1) + \Delta t f(t_1,y(t_1))\bigr)\biggr]

for which the right-hand-side is completely determined. Note that if \frac{}{}f(t,y(t)) = f(t), Heun's method is equivalent to the numerical integration method using trapezoids.

Example

Solve \frac{}{}dy/dt = 1 + yt with \frac{}{}y(0)=0.

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);
Personal tools