Numerical Methods for Differential Equations
[Pages:25]1
Numerical Methods for Differential Equations
1
2 NUMERICAL METHODS FOR DIFFERENTIAL EQUATIONS
Introduction
Differential equations can describe nearly all systems undergoing change. They are ubiquitous is science and engineering as well as economics, social science, biology, business, health care, etc. Many mathematicians have studied the nature of these equations for hundreds of years and there are many well-developed solution techniques. Often, systems described by differential equations are so complex, or the systems that they describe are so large, that a purely analytical solution to the equations is not tractable. It is in these complex systems where computer simulations and numerical methods are useful.
The techniques for solving differential equations based on numerical approximations were developed before programmable computers existed. During World War II, it was common to find rooms of people (usually women) working on mechanical calculators to numerically solve systems of differential equations for military calculations. Before programmable computers, it was also common to exploit analogies to electrical systems to design analog computers to study mechanical, thermal, or chemical systems. As programmable computers have increased in speed and decreased in cost, increasingly complex systems of differential equations can be solved with simple programs written to run on a common PC. Currently, the computer on your desk can tackle problems that were inaccessible to the fastest supercomputers just 5 or 10 years ago.
This chapter will describe some basic methods and techniques for programming simulations of differential equations. First, we will review some basic concepts of numerical approximations and then introduce Euler's method, the simplest method. We will provide details on algorithm development using the Euler method as an example. Next we will discuss error approximation and discuss some better techniques. Finally we will use the algorithms that are built into the MATLAB programming environment.
The fundamental concepts in this chapter will be introduced along with practical implementation programs. In this chapter we will present the programs written in the MATLAB programming language. It should be stressed that the results are not particular to MATLAB; all the programs in this chapter could easily be implemented in any programming language, such as C, Java, or assembly. MATLAB is a convenient choice as it was designed for scientific computing (not general purpose software development) and has a variety of numerical operations and numerical graphical display capabilities built in. The use of MATLAB allows the student to focus more on the concepts and less on the programming.
1.1 FIRST ORDER SYSTEMS
A simple first order differential equation has general form
?? ???????
??
(1.1)
???? where
means the c?hange in y with respect to time and ?? ??
is any function of y?a? nd time.
Note that the
d? erivat?iv! e of the variable, , depends upon itself. There are many different notations for , common ones include
and .
One of the simplest differential equations is
?? ????#"
?$
(1.2)
We will concentrate on this equation to introduce the many of the concepts. The equation is convenient because the
easy analytical solution will allow us to check if our numerical scheme is accurate. This first order equation is also
relevant in that it governs the behavior of a heating and cooling, radioactive decay of materials, absorption of drugs
in the body, the charging of a capacitor, and population growth just to name a few.
To solve the equation analytically, we start by rearranging the equation as
?? ?%?&"
??
(1.3)
and integrate once with respect to time to obtain '
( ? ? ?&" ?0)21
(1.4)
where C is a constant of integration. We remove the natural log term by taking the exponential of the entire equation
3547698A@CB ? 3EDGFIH
(1.5)
FIRST ORDER SYSTEMS 3
which finally can be written as
? ? 1 3?FIH $
(1.6)
You can check that this answer satisfies the equation by substituting the solution back into the original equation.
Since we obtained the solution by integration, there will always be a constant of integration that remains to be
? ?? specified. This constant (C in our
For simplicity of this chapter, we
awbiollvpersoocleuetidonw)iitshstpheeciinfiietidalbycoanndiintiiotinalthcoatn?di? t? io? n
or
the ?
initial state of the , yielding C=1.
system.
1.1.1 Discrete derivative
You should recall that the derivative of a function is equivalent to the slope. If you plotted the position of a car
traveling along a long, straight, Midwestern highway as a function of time, the slope of that curve is the velocity -
the derivative of position. We can use this intuitive concept of slope to numerically compute the discrete derivative
of a known function.
On the computer we represent continuous functions as a collection discrete, sampled values. To estimate the
slope (derivative) at any point on the curve we can simply take the change in rise divided by the change in run at
? ? ??? ?? ? any of the closely spaced points, and ,
?? ?? ?
? ?
" "
? ?
$
(1.7)
We can demonstrate this concept of the numerical derivative with a simple MATLAB script.
Program 1.1: Exploring the discrete approximation to the derivative.
t = linspace(0,2,20);
%% define a time array
y = exp(-t);
%% evaluate the function y = e^(-t)
plot(t,y);
%% plot the function
hold on
dydt = diff(y)./diff(t);
%% take the discret derivative
plot(t(1:end-1),dydt,'r--'); %% plot the numerical derivative
plot(t,-y);
%% plot the analytical derivative
This program simply creates variable t. The program then
eavalilsutaotefs2th0eefquunacltliyonsp? a?ced3
F
times between 0 H at these sample
and 2 points
and and
stores these numbers in plots the funct?ion. U? sing
the the
MATLAB diff command, we can evaluate the difference between neighboring points in the arrays and , which
? ? ? ? ipsouinstesdintoacloismt pouftneuamnbeesrtsim? a?te of? t? he d$Ae$ $r? iv? a( tivea.s Th? e?d? i?f f? co? m? m a" nd? s? im p? ly? t ak" es? t? he Cd$ i$Af$ ?fe? r( e nc" e? o?f( ne" igh b.orTinhge
resulting list is one element shorter than the original function. Finally, in the script we plot the numerical and
analytical function for the derivative. The plot that results from the script is shown in Figure 1.1. We see derivative
is approximate, but appears to be generally correct. We will explore the error in this approximation in the exercises
below and more formally in a later section.
1.1.2 Euler's method
We can use the numerical derivative from the previous section to derive a simple method for approximating the
solution to differential equations. When we know the the governing differential equation and the start time then we
know the derivative (slope) of the solution at the initial condition. The initial slope is simply the right hand side
#? %$ ! "? ioasnf d?E? pqrue?ad? tiioc? tnth1" e.1f..uOtIunurerF.ifigFrusotrrenthu1em.2cewariseceasolhfmotwehtehtohfuden,fckutninoocnwtio?nn?aas?n?Ed? u? tlheer" 'se? $ x,mt? rea?tp?hoo? lda,tiwo ni?llbuass,eetdhtheoinssloitnhpieetiaainltistthlioaelpiecnoittnoiadelixtcitoornanp.doiTtliaohtnee
extrapolation is valid for times not to far in the future (i.e.
), but the estimate eventually breaks down.
Another way to think about the extrapolation concept, is imagine you are in a car traveling on a small country
highway. You see sign stating the next gas station is 10 miles away, you look at your speedometer and it says you
are traveling 60 miles per hour. By extrapolation, you might predict that you will be at the gas station in 10 minutes.
The extrapolation assumes you will continue at your current speed (derivative of position) until you reach the next
gas station. If there is no traffic and your cruise control is working this assumption will be accurate. However,
your prediction will not be accurate if there are many stop lights along the way, or you get stuck behind a large,
4 NUMERICAL METHODS FOR DIFFERENTIAL EQUATIONS
1
0.8
0.6 y=e-t
0.4
0.2
y
0
-0.2
-0.4 dy/dt
-0.6
-0.8
-1
0
0.5
1
1.5
2
time
?????? ? Fig. 1.1 Graphical output from running program 1.1 in MATLAB. The plot shows the function
, the derivative of that
function taken numerically and analytically.
1
0.9
0.8
0.7
0.6
y
0.5
0.4
y=e-t
0.3
0.2
0.1 extrapolation of derivative
0
0
????? ? ? Fig. 1.2 Extrapolation of the function
0.5
1
1.5
time
based on the initial condition,
2
. For very short times, the estimate is
quite good, but clearly breaks down as we go forward in time.
slow-moving truck. Extrapolation is a good guess for where the system might be in a the near future, but eventually
the predictions will break down except in the simplest systems.
aanssdSuiminn? ipctetiiaoloncuoirsnpdgrioetoidodinc,,tfiiwogneusrkefnaooruwtinwthtoheetvrhaeeluwfeeuotaufrretehaeanrfdeuntnhcoetitnoanecxactnruadrpatothele,atsweloefopwrewialatlrtdthaeakgeinaisitnmia. alUltlismisnteeg,p?osu? wr ehqi.lueTahttiheoenvae?lxu? tera?ap?toa?laltai" toe?nr
time, where
, can be predicted
?? ??#"
the notation
by
H
extrapolations as
? ? ??
means the derivative
) of
?
?? ??
H
??
evaluated
" at
?!E $
time
equals
zero.
For
our
specific
(1.8) equation, the
extrapolation formula becomes
?$ ? ? " ? ? ? " ? $
(1.9)
This expression is equivalent to the discrete difference approximation in the last section, we can rewrite Equation
1.9 as
?%
?? ????
? ?
" "
?? ?
?&"
?$
(1.10)
?%&
' ( Once the value of the
call the time interval
function at over which
is we
known we can re-evaluate extrapolate, the time step
th?e
d? er? iva" tiv? e
and move forward to . Equation 1.9 is used
. We as an
typically iteration
1 [t ,y ]
00
0.8
FIRST ORDER SYSTEMS 5
y
0.6
[t ,y ]
11
0.4
Euler
y=e-t
[t ,y ]
0.2
22
[t ,y ]
33
[t ,y ]
0
44
0
0.5
1
1.5
2
time
Fig. 1.3 Graphical output from running program 1 in MATLAB. The points connected by the dashed line are the results of
?? the numerical solution and the solid line is the exact solution. The time step size is
?
? . This large time step size results in
large error between the numerical and analytical solution, but is chosen to exaggerate the results. Better agreement between the
numerical and analytical solution can be obtained by decreasing the time step size.
equation to simply march forward in small increments, always solving for the value of y at the next time step given
? #? ' Wst? he?eeeTktthhhneeaowntwrehtrhnseeuer-ielneetvfxwooatfrlerumathlapaationestildoatmhtnaiee.ottnshT(l1ohoo,dip0fse.ft2po,h5rrweo)o.hciuneiTcridthhmiuairisloesdspinlesrooolcpwceoeeqm,?sus??am?trieo?opn??nel? y?aut?scsia".nl"glWe$?ad? h,tEiaiglmneuedlteteshurseu'stseseepmrttorhseoiattzrhtheiosenlodopf.Fpoieignu?ttor? (ee01x.5.t$3r,? 0asp.i5esoe)lsmahatfoestweltrahnretghinneeeFtfixhirtgestutibrmtaeisme1ics.e3tte.srpteWentpode.
seems correct. As we make the time step size smaller and smaller the numerical solution comes closer to the true analytical solution.
A simple example of MATLAB script that will implement Euler's method is shown below. This program also plots the exact, known solution as a comparison.
Program 1.2: Euler's method for the first order equation.
clear; y = 1; dt = 0.5; time = 0; t final = 2; Nsteps = round(t final/dt); plot(time,y,'*'); hold on;
%% clear exisiting workspace %% initial condition %% set the time step interval %% set the start time=0 %% end time of the simulation %% number of time steps to take, integer %% plot initial conditions %% accumulate contents of the figure
for i = 1:Nsteps y = y - dt*y; time = time + dt plot(time,y,'*');
end t = linspace(0,t final,100); y = exp(-t); plot(t,y,'r') xlabel('time'); ylabel('y');
%% number of time steps to take %% Equation 1.9 %% Increment time %% Plot the current point
%% plot analytical solution
%% add plot labels
1.1.3 Evaluating error using Taylor series
When solving equations such as 1.2 we typically have information about the initial state of the system and would like to understand how the system evolves. In the last section we used the intuitive approach of discussing extrapolation. We simply used information about the slope, to propagate the solution forward. In this section we will place the extrapolation notion in a more formal mathematical framework and discuss the error of these approximations using the Taylor series.
6 NUMERICAL METHODS FOR DIFFERENTIAL EQUATIONS
? Consider an arbitrary function and assume that we have all the information about the function at the origin (
?
)
and we would like to construct an appro? ximation around the origin. Let's assume that we can create a polynomial
? ? approximation to the original function, , i.e. ??
& ) ?0) ? ? )2???? ) $ $A$
(1.11)
where we will need a method to solve fo? r the unknown coefficients a, b, c, d, etc. The simplest approximation for the would be to use the method of the last section to match the derivative of
the true and approximated functions, precisely we mean,
? ? ? ?? ?
)
?? ??
H
?
(1.12)
? where the
derivative
?? notation at the point
?
??#"
?
H
.
? means take the derivative of the function with respect to
and then evaluate that
We can improve the polynomial approximation by matching the second derivative of the real function and the
approximate function at the origin, i.e.
?
? ?? ? & ??
?? ? ?
??
?& H
H
?
? ?
? ? ??
?
? ?? ??
&
??
?
? &
?
H
H
$
(1.13) (1.14) (1.15)
If
we
continued
to
match ? ? ? ?
higher ? ?? ?
derivatives
0)
?
?? ??
of
H
the
)
tr? u& e
a?n&? d?&
approximated
H )
??
?
?? ?? ?
functions
H ) $ $ $
we ?6 (??
would
6? ?? 6
H
obtain
the
expression (1.16)
which is known as the Taylor series. Taylor series are covered in most Calculus text where you can find more detail,
examples, and generalizations.
?
To test this series we will return to our model function
?
3 F H . Substituting this function in Equation 1.16
yields the approximation as
? #? ? ? ? ?
& " ? )
?
"
?? )
?
$A$ $ ?"
6
?6 (??
(1.17)
?
The simplicity of this expression is due to the fact that all derivatives of evaluated at the origin are 1. In Figure
1.4 we plot? the first few terms of the series in comparis? on to the true function. We see that the approximation works well when is small and deviates for large values of . We also find that more terms included in the Taylor series
results in better agreement between the true and approximate functions.
'ini? tNiaolwcowndeirteiotunrntotothtehetimcoen?te?xt into the future. ?
o' f? ?'
the ,a
?
initial value problem. In the initial value
short ? ??
time
) '
into the
?
?? ??
H
future.
) '
T? & he?r&e? ?&forH e
we
)
problem evaluate
$C$ $
we the
want to move Taylor series
from the at a time
(1.18)
which can be rearranged as
?
?'
?
'
"?
?
?
?
?? ??
H
)
'
?
&
??
? &
H
)
$C$C$
(1.19)
which looks like Euler's method presented in the last section with the exception of the extra term on the right hand
side. This extra term? is the error of using Euler's method. Technically, the error is a series including terms t?o infinity
& ' of higher powers of (denoted by
' ? ' ! ? ' ? and therefore, the first error term
the repeating is dominant.
dots). We? When
re? tain$
only the then
fi? rst?
und$ er
the ,
?
a? ss? ump$ t? ion,tahnadt
so
is small on. We
' only need to retain the first neglected term in the Taylor series? to understand the error, when is small. Equation
' 1.18 shows that the error in Euler's method will scale with . By the scaling, we mean that if we halve the time
' step then we halve the error. In later sections we will introduce m? ethods that exhibit better scaling, if we halve the
time step then the error might decrease as a higher power of . Knowing the scaling allows one to check that
FIRST ORDER SYSTEMS 7
1
0.9
0.8
0.7
0.6 1+t+t2/2
0.5
exp(t)
0.4
0.3
1-t+t2/2-t3/6
0.2
1+t
0.1
0
0
0.2
0.4
0.6
0.8
1
t
y(t)
??? ? Fig. 1.4 Taylor series approximation of the function . We have included the first three terms in the Taylor expansion. It is
clear that the approximation is valid for larger and larger x as more terms are retained.
10-1
10-2
~ t 10-3
error
10-4
error
10-5 10-4
10-3
10-2
t
10-1
Fig. 1.5 Error of the Euler method as we change the time step size. For comparison we show the linear plot of . Since the slope of the error curve matches the slope of the linear curve we know that the error scales as .
there result is behaving as expected. The knowledge of the scaling is also used in more complex methods that use
adaptive control of the time step size to ensure that the solution is within acceptable error bounds.
One way to confirm the scaling of the numerical methods is to plot the error on a log-log plot. In Figure 1.5 we
plot the error in ap?plying Euler's method to our model equation as a function of the time step size. We also plot
' the line, error ' plot follows a
?
. power
We find law, i.e
that both error
?
li? n6 e.s
have the same slope. The s? lope of the line
A straight line on a a provides the power.
log-log We find
plot means that the from this plot that
' our Euler method error is scaling linearly with as the slopes of the two displayed curves match. This graphical
result agrees with the prediction of the error using Taylor series.
To provide a little insight into how such graphical results can be easily generated we present the program that
created Figure 1.5. The program is a simple modification of Program 1.2. In the program we simply loop the Euler
solver for several time step sizes and store the values of the error for plotting. The error is defined as the difference
8 NUMERICAL METHODS FOR DIFFERENTIAL EQUATIONS
??b? e? twe,ebnutthtehetrtuime eantdhantuwmeeirnicteaglrsaoteluutinotnilsisatirtrheeleevnadnto. f the integration period. In this example we integrate until
Program 1.3: Program to check error scaling in Euler method. The result is shown in Figure 1.5.
clear;
%% clear exisiting workspace
dt = [0.0001 0.0005 0.001 0.005 0.01 0.05 0.1 ];
for j = 1:length(dt)
y = 1;
%% initial condition
time = 0;
%% set the time=0
t final = 1.;
%% final time
Nsteps = round(t final/dt(j));
%% number of steps to take
for i = 1:Nsteps
y = y - dt(j)*y;
%% extrapolate one time step
time = time + dt(j);
%% increment time
end
X(j,2) = exp(-t final) - y;
%% compute the error and store
X(j,1) = dt(j);
%% store the time step
end
loglog(X(:,1),X(:,2))
%% display on log-log plot
1.1.4 Programming and implementation
In this section we provide a few different ways to create the program that implements Euler's method. Program 1.2 is a perfectly acceptable way to solve Euler's equation. Since we will be interested in applying this same basic method to different sets of equations there is some benefit in breaking the program into functions. While the benefit with something as simple as Euler's method applied to first order systems is not substantial, this programming effort will become worthwhile as we proceed in the next section and deal with systems of equations. We proceed with a demonstration of several programs that are equivalent to Program 1.2. It is recommended to try these programs as you proceed through the chapter, making sure to understand each program before proceeding to the next. Consult the section on programming functions in MATLAB in order to understand the syntax of the programs presented.
One possible program is provided in Program 1.4. This program is makes use of different functions to keep the Euler algorithm separate so that it only needs to be programmed once, and then left alone for all other equations. Program 1.4 has the advantage in that the to solve a different system you do not need to modify lines of code that deal with the Euler method.
The program has the program broken into two functions. The functions can then be called from the command line. For example, the command euler1storder(1,0.01,1) means the function will solve the first order equation with an initial condition of 1, a time step size of 0.01, until a time of 1. The first function, derivs, contains the governing differential equation we would like to solve. The input to this function is the current value of y and time and the output of the function is the derivative, or the right-hand-side of the differential equation. The function euler1storder contains all the code directly related to Euler's method. The benefit of this arrangement is a different equation is easy to solve without the risk of corrupting the Euler method itself. You would never need to modify the code inside the function euler1storder.
Program 1.4: Use of functions to generalize the Euler program. All the code should be in one m-file named modeleqns.m.
%% derivative functions function dy = derivs(time,y) dy = -y; %% governing equation
%% Euler Solver function euler1storder(y,dt,t final) clf; Nsteps = round(t final/dt); %% number of steps to take time = 0;
for i =1:Nsteps dy = derivs(time,y); y = y + dy*dt; time = time+dt;
%% compute the derivatives %% extrapolate one time step %% increment time
................
................
In order to avoid copyright disputes, this page is only a partial summary.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related download
- unit 2 differentiation lesson 1 the derivative
- 5 numerical differentiation
- 3 1b derivative of a function
- chapter 4 fourier series and integrals
- a reduced frequency approach for calculating dynamic
- difference quotient csusm
- of derivative jj ii
- numerical methods for differential equations
- derivative of 1 by root x
- dn1 3 gradients tang ents and derivatives
Related searches
- differential equations sample problems
- differential equations problems and solutions
- differential equations practice problems
- differential equations review sheet
- differential equations formula sheet pdf
- differential equations review pdf
- differential equations cheat sheet pdf
- differential equations pdf free download
- linear differential equations problems
- introduction to differential equations pdf
- linear differential equations definition
- solving ordinary differential equations pdf