Numerical Methods

Index of Unit VI Solution of Differential Equations

Introduction........................................................................................................................................1
6.1 One-step methods.................................................................................................................2
6.1.1 Euler's Method
[Link] Error Analysis for the Euler Method ............................................................... 5
[Link] Improved Euler Method
6.1.2 Method for Higher Order Taylor Series ............................................................... 7
6.1.3 Midpoint Method (or Improved Polygon Method)............................................................ 8
6.1.4 Runge–Kutta Method ...................................................................................................... 9
[Link] Second-order Runge-Kutta methods..................................................................... 11
6.2 Métodos de Pasos Múltiples......................................................................................................16
6.2.1 Heun's Method of Non-Auto Start ..................................................................................... 17
6.2.2 Higher Order Multi-Step Methods ............................................................................... 18
6.2.3 Milne Method.................................................................................................................. 19

6.2.4 Fourth order Adams method

6.3 Systems of Ordinary Differential Equations

6.3.1 Euler's Method................................................................................................................... 20

6.4 Applications
Appendices
Appendix A "Investigated Methods"............................................................................................ 32

Appendix B "Examples".................................................................................................................. 33

Conclusion

Bibliografía........................................................................................................................................ 47
 Unit VI 'Solution of Differential Equations'

Introduction

Numerical methods make us capable of understanding schemes.

numerical in order to solve mathematical, engineering and
scientists on a computer, reduce basic numerical schemes,
write programs and solve them on a computer and use it correctly
existing software for such methods not only enhances our ability
for the use of computers but also broadens mathematical expertise
and the understanding of basic scientific principles.

But this time we will apply numerical methods for equations.

differentials and their methods for problem-solving through the
one-step methods; as well as through Step methods
Multiple.

Portillo Contreras Misael ISC Group 1 Page 1

 Unit VI 'Solution of Differential Equations'

6.1 One-Step Methods

All one-step methods can be expressed in this general form,

that will only differ in the way the slope is estimated. As in
the problem of the falling skydiver, the simplest procedure is to use the
differential equation to estimate the derived slope atxiat the beginning of the
interval. In other words, the slope at the beginning of the interval is taken
as an approximation of the average slope over the entire interval.
This procedure, called the Euler method.

6.1.1 Euler Method

The first derivative provides a direct estimate of the slope at

xI

f(xI,yi) Wheref(xi,yi) it is the differential equation evaluated atx I y yi

it can be substituted in the equationyi 1 yi hwhere it will result for us
next equation:
y i 1 yi f  xi,yI  h
This formula is known as the Euler method (or Euler-Cauchy method)
midpoint). A new value of y is predicted through the slope.
(equal to the first derivative at the original value of x what will have to
linearly extrapolate over the step size h

Portillo Contreras Misael ISC Group 1 Page 2

 Unit VI 'Solution of Differential Equations'

Problem statement. (Example)

Use the Euler method to numerically integrate the equation

dy
2
 x 3 12x2 20x 8.5
dx

Sincex 0 untilx 4with a step size of 0.5. The initial condition at

x 0 it y 1 Remember that the exact solution is given by the equation

y 0 .5x4 4x 3 10x2 8.5x 1

Solution. The equation can be [Link] 1 yi f  xi,yi  h to implement

the Euler method:
y 0.5 y 0  f  0.10. 5
Where y 0 1 and the estimated slope atx 0 is
f  0,1  0  3 12
2  0 2 20  0  8.5 8.5

y 0.5 1.0 8.5  0.5  5.25

The real solution inx 0.5it

y 0 .5  0.5 4 4  0.5 3  10  0.5 2 8.5  0.5 1 3.21875

Portillo Contreras Misael ISC Group 1 Page 3

 Unit VI 'Solution of Differential Equations'

Portillo Contreras Misael Computer Science Group 1 Page 4

 Unit VI 'Solution of Differential Equations'

Thus, the error is:

Et =true - approximate = 3.21875 - 5.25 = -2.03125 or

Expresada como error relativo porcentual, et .63.1%For the second step,
y 1 y 0.5  f  0.5, 5.250.5

5.25   2  0.5 3 12  0.5 2  20 0.5
   8.5 0.5 5.875

[Link] Error analysis for the Euler method

The numerical solution of ODEs involves two types of error (remember the
chapters 3 and 4)
1. Errors of truncation, or deratization, caused by the nature of the
techniques used to approximate the values of y.
2. Rounding errors, which are the result of the limited number of digits
significant that a computer can retain.

Truncation errors consist of two parts. The first is an error of

Local truncation resulting from an application of the method in question on a
simple step. The second is a propagated truncation error that results from the
approximations produced during the previous steps. The sum of the two is the
total, global truncation error.

One can obtain certain knowledge about the magnitude and properties of the
truncation error when deriving the Euler method directly from the expansion
from the Taylor series. To do this, note that the differential equation is subject to
integration will be in the form
General:y' f(x,y)

dy
Wherey' ,y x y y they are the independent and dependent variables,
dx
respectively. If the solution (that is, the function that describes the
behavior of y) has continuous derivatives, can be represented by a
expansion of the Taylor series with respect to the initial values  xi,yI  like in
yi ´´ 2 yi(n) n
y i 1 yi yi ´h h  h  Rn
2! n!

Portillo Contreras Misael ISC Group 1 Page 5

 Unit VI 'Solution of Differential Equations'

y (n1)( ) n1
Where: h xi1x yR remaining term, defined as Rn h
i n
(n 1)!
Where is somewhere in the interval ofxia xi 1

[Link] Improved Euler Method

This method is based on the same idea as the simple Euler method, but it makes a
refinement in the approach, taking an average of the slopes
of the found tangent lines.

Thus, in the graph we see that the average slope m corresponds to the
pending the bisector line of the tangent line to the curve at the point of
the initial condition and the tangent line to the curve at the point ,  x1 ,y1 
wherey1It is the approximation obtained with the first Euler's formula.
Finally, this bisector is moved parallel to the point of
the initial condition, and the value of this line is considered at the
pointx x1like the improved Euler approximation.

Portillo Contreras Misael ISC Group 1 Page 6

 Unit VI 'Solution of Differential Equations'

The approximation at each step is then determined by the formula:

yi 1 yi
f  xi,yi   f  x  1,y  1
i i
0
h Being: yn1 y h  nf  xn,yn 
2

6.1.2 Method for the higher-order Taylor series

One way to reduce the error of Euler's method could be the inclusion of
higher order terms in the expansion of the Taylor series for its solution.
For example, with the inclusion of the second-order term as follows
equation:
f'(xi,yI)
y i 1 y i f  xI,yi  h h2
2!

A local truncation error

f''(xI,yi)
Ea h3
6
Although the incorporation of higher-order terms is sufficiently
simple to implement in the polynomials, its inclusion is not so trivial when
The ODE is more complicated. In particular, the ODEs that are a function of both the
dependent variable as well as the independent, require differentiation by the
chain rule. For example, the first derivative of f(x,y) it

f x,y  f x,y  dy
f ´ xi,yi  
x y dx

The second derivative is:

f f dy f f dy
f 
x y dx x y dx dy
f´  xi,yi  
x y dx

Portillo Contreras Misael ISC Group 1 Page 7

 Unit VI 'Solution of Differential Equations'

Higher order derivatives become much more complicated.

As a consequence, as described in the following sections, has been developed.
alternative one-step methods, These schemes are comparable in performance
with the procedures of the high-order Taylor series, but they require only
from the calculation of the first derivatives.

6.1.3 Midpoint method (or improved polygon method)

Another simple modification of the Euler method. Known as the point method.
average (of the improved polygon or the modified version of Euler), this technique uses the
Euler's method to predict a value of y at the midpoint of the interval (see
the following figure

Portillo Contreras Misael ISC Group 1 Page 8

 Unit VI 'Solution of Differential Equations'

h
yI 1
yI f(xI,yi)
2 2
Then, this predicted value is used to calculate a slope at the midpoint:

yI 1 f x ,yi 1 i 1
2 2 2

which is taken to represent a valid approximation of the slope

average for the entire interval. This slope is used later for extrapolation.
linearly from x i until xi1

yi 1 yi f x ,y h i 1 i 1
2 2

As in the previous section, this procedure can also connect with the
Newton-Cotes integration formulas

6.1.4 Runge-Kutta Method

The objective of the Runge-Kutta numerical methods is the analysis and solution of
the initial value problems of ordinary differential equations (ODE), these
they are an extension of the Euler method to solve the (ODEs), but with a
order of higher accuracy than this.

Runge-Kutta (RK) methods achieve the accuracy of the procedure of a

Portillo Contreras Misael ISC Group 1 Page 9

 Unit VI 'Solution of Differential Equations'

Taylor series without requiring the calculation of higher derivatives. There are many
variations, but all can be denoted in the generalized form of the equation

yI1 yI  xi,yi, h  h

Where  xi,yi, h  is known as the increment function, which can
to be interpreted as a representative slope over the interval. The function
increment is generally written as:

a1k1 a2 k 2  an k n

Where the a's are constants and the k's are:

Note that the k they are recurrence relations. That is,k1appears in the
equation fork 2, which appears in the equation for k 3, etc. As each k is a
functional evaluation, this recurrence causes the methods be efficient
for calculations on computer.
It is possible to conceive several types of Runge-Kutta methods by using different
numbers of terms in the increment function as specified by n.
Note that the first-order Runge-Kutta (RK) method with n 1yes, in fact,
the Euler method.

Once n is chosen, the a, p, and q are evaluated by equating equation 25.28 to the
terms in the Taylor series expansion. Thus, at least for the versions of

Portillo Contreras Misael ISC Group 1 Page 10

 Unit VI 'Solution of Differential Equations'

lower order, the number of terms n often represents the order of the
approximation. For example, the second-order RK methods use the function
increment with two terms Those second-order methods will be
exact if the solution of the differential equation is quadratic. Moreover, since the
terms withh 3and older ones are eliminated during the derivation, the error of
local truncation is and the global is In subsequent sections
we will develop third and fourth order RK methods For those cases,
global truncation errors are , respectively. y

[Link] Second-order Runge-Kutta methods.

The second order version of the previous equation is:

yI1 y i  a1k1 a2 k 2  h

Where:

k1 f  xi yI
k2 f  xI p1h yi q11k1h 

When using the equation, we must determine the values for the constants a1, a2, p1.
y p11. To do this, we recall that the second-order Taylor series for i1in y
terms of y t y f(x,y)
i i it is written as:

f'(x)i,yi) 2
yI1 yi f  xi,yi  h  h ecu. 1
2!

Where it must be determined by differences using the chain rules

Portillo Contreras Misael ISC Group 1 Page 11

 Unit VI 'Solution of Differential Equations'

f x,y  f x,y  dy
f ´ xi,yi  
x y dxecu. 2

If we substitute equation eq. 2 into equation eq. 1, we have

f x,y  f x,y  dy h 2
yi 1 yI f  xI,yi  h 
x y dx 2!

The basic strategy that will be highlighted in the Runge-Kutta methods is the use of
of algebraic manipulations to solve the values of which
it causes the equations

yi1 y i  a1k1 a2 k 2  h
and the previous ones are equivalent. To do this, we first use a Taylor series to
expand the equation.k 2 f x i p1h yiq11k1hThe  Taylor series for a
A function of two variables is defined as:

g g
e x r y s  g x,y   r s
x y

If this method is appliedpara ex

pandto go the
equation yi1 y i a 1k1  a2 k 2  hhas

f f
f  xI p1h yi q11k1h  f  xi,yi   p1h q11k1h O h 
2
x y

This result can be replaced un

j toon
c la
equationk1 f  xi,yI  y yi1 ya i1k1  a2k2 h
to give

Portillo Contreras Misael ISC Group 1 Page 12

 Unit VI 'Solution of Differential Equations'

f f
yi1 yi a1hf  xi,yi   a2 hf  xi,yI  a2p1h 2  a2 q11h 2f  xi,yI   
 O h3
x y

Oh, when grouping terms,

f f 2
yi1 yi a1f  xiyI   a 2f  xi,yi  h  a 2p 1
 a2 q11f  xI,yi  h O h 
3
x y

Now, if we compare common terms in the previous equations

we determined that to make the two equations equivalent, the following must be fulfilled
the following:

a1 a 2 1
1
a1p 2
2
1
a1 q11
2

The previous three simultaneous equations contain the four constants.

unknown. As there is one more unknown than the number of equations, there is no solution.
a unique set of constants that satisfy the equations. However, when
Suppose a value for one of the constants, we can determine the other three. In
As a consequence, there exists a family of second-order methods rather than just one.
version.

Since we have three equations with four unknowns, we must assume the value of
one of these unknowns to determine the other three. Suppose we specify
a value for a2. Then the equations can be solved simultaneously.
25.31 to 25.33 to obtain:

Portillo Contreras Misael ISC Group 1 Page 13

 Unit VI 'Solution of Differential Equations'

Since we can choose a finite number of values for a2, there is a number
interminable of second-order RK methods. Each version could give
exactly the same results if the solution of the ODE were quadratic, linear
or a constant. However, different results are obtained when the
the solution is more complicated. Below we present three of the most
commonly used and preferred:

Heun's method with only corrector (a2 = ½). If we assume that a2 is 1/2, the
Equations (25.34) and (25.35) could be solved for a1 = ½ and p1 = qI 1 = 1. These
parameters, when substituted into the equation (25.30), yield

Where

Note that k1 is the slope at the beginning of the interval and k2 is the one at the end.
consequence, this second-order Runge-Kutta method is indeed the technique
of Heun without iteration.

The midpoint method (a2 = 1). If we assume that a2 is 1,

then , and the equation is now

Portillo Contreras Misael ISC Group 1 Page 14

 Unit VI 'Solution of Differential Equations'

Where

This is the midpoint method.

Ralston method . Ralston and Rabinowitz determined that when
select a minimum limit on truncation error is obtained
for second-order RK algorithms. For this
version

Where

Portillo Contreras Misael ISC Group 1 Page 15

 Unit VI 'Solution of Differential Equations'

6.2 Multiple Step Methods

The one-step methods described in the previous sections use information

at a single point xi to predict a value of the dependent variable yi+1 in a
future point xi+1. Alternative procedures, called multi-step methods, are
they are based on the knowledge that once the calculation begins, it has
valuable information from the previous points is at our disposal. The
the curvature of the lines connecting those previous values provides information
regarding the trajectory of the solution. The multi-step methods that
we will explore how to utilize this information to solve the ODEs. Before
describe the higher-order versions, we will present a simple method of
second order that serves to demonstrate the general characteristics of the
multi-step procedures.

Graphic illustration of the fundamental difference between the methods for

solve ODE a) in one step and b)
Multiple steps.

Misael Portillo Contreras ISC Group 1 Page 16

 Unit VI 'Solution of Differential Equations'

6.2.1 Heun's method of Non-Self-starting

Let us remember that Heun's procedure uses Euler's method as a

predictor:
y 0yi1 yi'f(xi,yi)h
And the trapezoidal rule as a corrector:

f(xi yi)  f(xi1 yI01 )

yi1 yi hec.1
2

Thus, the predictor and the corrector have local truncation errors
of, y respectively. This suggests that the predictor is the weak link
in the method, as it has the biggest error. This weakness is significant due to
that the efficiency of the iterative corrector step depends on the accuracy of the
initial prediction. Consequently, one way to improve Heun's method is
through the development of a predictor that has a local error of This is
can be fulfilled by using Euler's method and the slope at, and information
extra from the previous point like in:

yi01 y  f(x
i1 y i
2hi ec.2

Observe that equation ec. 2 reaches at the expense of using a size of

largest step, 2h. Furthermore, note that equation eq. 1 is not self-starting,
since it involves a previous value of the dependent variable yi-1. That value might not
to be available in a common initial value problem. Because of this, the
Equations are called Heun's method of non-self-initiation. The estimated derivative
from the equation is now located at the midpoint rather than at the beginning of the interval
about which the prediction is made. As will be demonstrated later, this location
centered improves the predictor's error to However, before proceeding to
a formal deduction of the Heun method of non-self-starting, we will summarize the
method and we will express it using a slightly modified nomenclature:

Predictor:
Portillo Contreras Misael ISC Group 1 Page 17
 Unit VI 'Solution of Differential Equations'

Corrector:

Where the superscripts were added to denote that the correction is applied
iteratively from j 1 to obtain refined solutions. Observe
m m
what y i&y i 1 these are the final results of the iterations of the corrector in
los pasos de tiempo anteriores. Las iteraciones son terminadas en cualquier paso
of time based on the stop criterion:

yI1´  y I1
Ea 100% ec. 3
i1
y i
I1

WhenEais less than an error tolerance It is preset, it ends the

iterations. At this pointj m

6.2.2 Higher-order Multi-step Methods

Now that we have formally developed the integration formulas of

Newton-Cotes and Adams, we can use them to deduce multipoint methods of
higher order. As happened with the non-autonomous Heun method, the
integration formulas are applied in series as predictor-corrector methods.
Moreover, if the open and closed formulas have local truncation errors
In the same order, it is possible to incorporate listing-type modifiers. To improve
the accuracy and allow control of step size. Provides equations
generals for those modifiers. In the following section, we present two of the
most common higher-order multistep procedures: the Milne method and
the fourth order Adams method.

Misael Portillo Contreras ISC Group 1 Page 18

 Unit VI "Solution of Differential Equations"

Milne's Method.

The Milne method is the most common of the multiphase methods based on the
Newton-Cotes integration formulas. Use the three-point Newton-Cotes formula.
points like a predictor:

and the closed formula of Newton-Cotes of three points (Simpson's 1/3 rule) as
a corrector:

Fourth-order Adams method:

A popular multi-step method based on Adams integration formulas

use the fourth-order Adams-Bashforth formula (see table 26.1) as the
predictor:

and the fourth-order Adams-Moulton formula as the corrector:

Portillo Contreras Misael ISC Group 1 Page 19

 Unit VI 'Solution of Differential Equations'

The predictor and corrector modifiers for the fourth-order Adams method
They can be developed from the formulas and the error coefficients.

6.3 Systems of Ordinary Differential Equations

Many practical problems in engineering and science require the solution of a

system of simultaneous ordinary differential equations more than just one
equation. Such systems can generally be represented as:

dy1
f1(x y1 y 2 yn
dx
dy2
f 2(x y1 y 2 yn
dx

dyn
f n(x y1 y 2 yn
dx
The solution of such a system requires that the n conditions are known.
initials in the initial value of x.

Portillo Contreras Misael ISC Group 1 Page 20

 Unit VI "Solution of Differential Equations"

Note.- The methods used for the solutions of these systems of equations are the
used in the previous sections, therefore we will proceed to adjust the step size
directly, of course after having solved the system using one of the methods
seen before
Step size control.

Now that we have developed ways to estimate local truncation error, we

can be used to adjust the step size. In general, the strategy is to increase
increase the step size if the error is too small and decrease it if it is too large.
Press and Cols. (1992) have suggested the following criterion to comply with the above:

Where h-actual and h-new = current and new step size, ∆actual = accuracy
actual calculated, ∆new = desired accuracy, and a = constant exponent that is
equal to 0.2 when the step size is increased and 0.25 decreases the step size
step.

The key parameter in equation 25.47 is obviously ∆new since it is its

vehicle to specify the desired accuracy. One way to do this would be
relate ∆ new to a relative level of error. Although this works well only
when positive values occur, it can cause problems for solutions that
they pass through zero. For example, you could be simulating an oscillatory function.
which repeatedly passes through zero, but is limited by maximum values
absolutes. In such a case, I might need these maximum values to appear in the
desired accuracy.

A more general way to handle those cases is to determine the new ∆ as:

Portillo Contreras Misael ISC Group 1 Page 22

 Unit VI 'Solution of Differential Equations'

Where E = global tolerance level. Your choice of y-scale will determine then.
how the error has escalated. For example, if y-scale = y, the accuracy will be
managed in terms of fractional relative error. If you now deal with a case
where relative errors are desired constant to a pre-established maximum limit, there exists
a y-scale equal to that limit. A trick suggested by Press et al. To obtain
the constant relative errors except those that cross very close to zero, is:

Portillo Contreras Misael ISC Group 1 Page 23

 Unit VI 'Solution of Differential Equations'

6.4 Applications

Euler's method and modified Euler's method a circuit contains a

impedance, a resistance, and a capacitance, the equation that governs this problem
"LRC" when the system is not subject to any potential is of type:

An L reactance of 0.4H, R=300Ω, and a circuit characteristic will be taken into account.
capacitance of .001 F. At the initial time (t=0), the current is 3A and its derivative
(that is, the electric charge) of .5A/s. °C Solution First, this must be transformed
problem in a set of first-order equations. Q will be taken equal to the
derivative of current intensity.

If the traditional Euler method is used, the equations must be solved.

using the formulas:

Portillo Contreras Misael ISC Group 1 Page 24

 Unit VI 'Solution of Differential Equations'

The results table obtained with a step of .0005 is:

If now the modified Euler method is used, the formulas are:

Portillo Contreras Misael ISC Group 1 Page 25

 Unit VI 'Solution of Differential Equations'

It is important to emphasize that the problem becomes very unstable if higher values are used.
highs for L.

Butcher's method: implicit second order

Let the following PVI be:

Y|= .3y+et =f(t , y)
Y(0) = 1
Solve this problem using the 2nd order Runge-Kutta method.
built from the following Butcher matrix:

Solution:
It is worth noting that the above scheme is implicit as it is a dense matrix A.
Applying the generic second-order Runge-Kutta formulas to the array of
Butcher anterior falls:

Portillo Contreras Misael ISC Group 1 Page 26

 Unit VI 'Solution of Differential Equations'

Substituting in the function f by the expression from the example, the following remains
algoritmo de cálculo:

Note that it is now necessary to solve a system of equations in K1 and K2 for

every step of time.

Calculations begin with i=0, t=0, y0=1, that is, the initial value is assumed to be
temporal step value h=0.1. The sequence of the consequent calculations is
summary in the table below.

Portillo Contreras Misael ISC Group 1 Page 27

 Unit VI 'Solution of Differential Equations'

Rigid equation system and stability (SisRigid)

Let the following coupled system of first-order differential equations be:

Coya writing in matrix form leads to:

Solution:
To find an analytical solution to the problem, it is necessary to diagonalize the matrix A.
or decoupling the system of equations through a similar transformation. For

Portillo Contreras Misael ISC Group 1 Page 28

 Unit VI 'Solution of Differential Equations'

This requires calculating the eigenvalues and eigenvectors of the matrix A.

The eigenvalues are obtained by setting the determinant of |A-λI| equal to zero, which
results in the following quadratic equation:

And the nuance of the corresponding eigenvectors is:

Therefore, through the following change of variables:

The previous system is transformed from a decoupled one:

And the analytical solution is now immediate:

Portillo Contreras Misael ISC Group 1 Page 29

 Unit VI 'Solution of Differential Equations'

That by substituting the original variables and the initial conditionals, it remains
finally

If y1 is graphed against t, two time scales, t1 and t11, are observed where the
first term of the solution, e-11.1t, remains practically constant at the
long of t1 and then slowly decays in a time interval 0<t11<40.

It is said that the ODE system is stiff when there is an eigenvalue in the matrix.
of the system that contributes almost nothing to the solution, especially regarding the domain of
integration. A stiffness index is defined as follows:

Portillo Contreras Misael ISC Group 1 Page 30

 Unit VI "Solution of Differential Equations"

Where |Re|λ|| represents the modulus of the eigenvalue part. Indeed, let it be the
generic method given by the following difference equation, applied to f(t,
y)=ky=y|:

Where µ(hλ) represents the characteristic polynomial of the difference equation,

It is said that the method is stable if it meets the condition, |µ(hλ)|≤1
for all the values of hλ.

Portillo Contreras Misael ISC Group 1 Page 31

 Unit VI 'Solution of Differential Equations'

Appendix A 'Investigated Methods'

Misael Portillo Contreras ISC Group 1 Page 32

 Unit VI 'Solution of Differential Equations'

Appendix B "Examples"

One-Step Methods

Euler's Method

1. Use the Euler method to numerically integrate the equation.

dy
2
 x 3 12x2 20x 8.5
dx

Sincex 0 untilx 4with a step size of 0.5. The initial condition at

x 0 is y 1 Remember that the exact solution is given by the equation.

y 0 .5x4 4x 3 10x2 8.5x 1

Solution. The equation can be used yi 1 yi f  xi,yi  h to implement

the Euler method:
y 0.5 y 0  f  0,10. 5
Where y 0 1 and the estimated slope atx 0 es
f  0.1  0  3 12
2  0 2 20  0  8.5 8.5

y 0.5 1.0 8.5  0.5  5.25

The real solution inx 0.5is

y 0 .5  0.5 4 4  0.5 3  10  0.5 2 8.5  0.5 1 3.21875

Portillo Contreras Misael ISC Group 1 Page 33

 Unit VI 'Solution of Differential Equations'

Heun's method

2. Use the Heun method to numerically integrate the following equation.

Differential

dy
4e0.8x 0.5y
dc
From x=0 to x=4, with a step size h=1. With the initial condition that
When x=0 then y=2. Obtain the exact solution by integrating analytically.
And compare the results with those obtained by the Heun method.

Analytical solution of the differential equation:

y' 4e0.8x 0.5y
Initial conditions:
x 0,y 2
The differential equation in y' 4e0.8x 0.5y it is in the form
y'day g(x)
The proven solution for the differential equation is:
y this a x
In y this a x, c it is arbitrary, the solution represents an infinity of
solutions according to the value c.

Comparing y this a xiswith y' 4e0.8x 0.5y it is found that in a 0.5

The equation is rearranged. y' 4e0.8x 0.5y agreed oh then it multiplies the
result
e 0.5x of the solution y this to simplify the first member:
ax
By the term x
y' 4e0.8x  0.5y

Solution. Before solving the problem numerically, we can use

calculation to determine the following analytical solution

Portillo Contreras Misael ISC Group 1 Page 34

 Unit VI 'Solution of Differential Equations'

Runge Kutta method

Use the fourth-order RK method to solve the ODE

Solution. First, we need to solve for all the slopes at the beginning of
interval
K i1 f(0,4,6) 0 .5(4) 2
K i2 f(0,4,6) 4 0.3(6) 0.1(4) 1.8

WhereK iit j is the i-th value of the j-th dependent variable. Afterwards,
we must calculate the first values of deyxyy2 at the midpoint:
h 0.5
y1 k1 1 4 (2) 3.5
2 2
h 0.5
y 2 k1 2 6 (1.8) 6.45
2 2
These will be used to calculate the first set of point slopes.
half

K 2.1 f(0,25,3.5,6.45) 1 .75

K 2.2 f(0,25,3.5,6.45) 1.715
These are used to determine the second set of point predictions.
half
h
y1 k 2.1 4   1.75  0.5 3.5625
2 2
h
y 2 k 2.2 6  1.715 0.5 3.5625
2 2
which will be used to calculate the second set of point slopes
mid
k3.2 f  0.25,3.5625,6.42875 1
 .78125
k3.2 f  0.25,3.5625,6.42875 1.715125

Portillo Contreras Misael ISC Group 1 Page 35

 Unit VI "Solution of Differential Equations"

These will be used fortterminate the predictions at the end of the interval
y1 k 3.1h 4   1.781250  .5  3.109375
y 2 k 3.2 h 6   1.7151250  .5  6.857563
Which will be used to calculate the slopes at the end of the interval,

k4.1 f  0.5,3.109375,857563 1
 .554688
k4.2 f  0.5,3.109375,857563 1.631749

1
y1 (0.5) 4    2  2   1.75 1.78125   1.554688  0.5 3.115234
6
1
y20.5 6   1.8 2  1.75 1.78125 1.631794
  0.5 6.857670
6
By proceeding in a similar manner for the remaining steps, one obtains

Portillo Contreras Misael ISC Group 1 Page 36

 Unit VI "Solution of Differential Equations"

Second Order Runge Kutta Method

Comparison of several second-order RK schemes.

Statement: Use the midpoint method and Ralston's method for integration
numerically the equation:

From until using a step size of 0.5. The initial condition

in Compare the results with the values obtained from another.
second order RK algorithm: the Heun method without iteration corrector.

Solution: The first step in the midpoint method is the use of the equation to
calculate:

However, since the EDO is a function only of x, such a result lacks

relevance of the second step to calculate:

Note that such an estimation of the slope is much closer to the value
Average for the interval (4.4375) where the slope at the beginning of the interval (8.5)
that could have been used by the Euler procedure. The slope at the
The midpoint can then be substituted in equation 25.37 to predict.

The calculation is repeated, and the results are summarized in figure 25.14 and table 25.3.

Misael Portillo Contreras ISC Group 1 Page 37

 Unit VI "Solution of Differential Equations"

Through the Ralston method, k1 for the first interval is also equal to 8.5
y:

The average slope is calculated by

La cual se usara para predecir

Note that all second-order RK methods are superior to the method of

Euler

Portillo Contreras Misael ISC Group 1 Page 38

 Unit VI 'Solution of Differential Equations'

Multi-Step Methods

Heun's Autostart Method

Use the Heun method of non-autostart to perform the same calculations as

in example 25.5 using Heun's method. That is, to integrate
of using a step size of 1.
As in example 25.5, the initial condition in However,
As we are dealing with a multi-step method here, we require the information.
additional to that .

Solution: The predictor is used to linearly extrapolate from

The corrector is then used to calculate the value:

Which represents a relative percentage error of -5.73%. This error is somewhat more
smaller than the value of -8.18% incurred in the auto-start Heun.
Now, the predictor equation can be applied iteratively to improve.
the solution:

Portillo Contreras Misael ISC Grupo 1 Page 39

 Unit VI 'Solution of Differential Equations'

What does an Et of -1.92% represent? An error estimate can be determined.

approximately using the equation eq. 3:

The equation can be applied iteratively until Ea is below

a pre-specified value of Es. As was the case with Heun's method, the
Iterations converge to a value of 6.360865. However, as the value of
initial predictor is more accurate, the multi-step method converges somewhat
faster.

For the second step, the predictor is:

What is superior to the prediction of 12.08260 that was calculated with the method of
Heun original. The first corrector gives 15.76693 and subsequent iterations.
they converge on the same result as obtained with Heun's method
Autostart: 15.30224. As with the previous step, the rate of convergence of
The corrector has been improved due to better initial prediction.
Deduction and analysis of the error of the predictor-corrector formulas. Already
we use graphic concepts to deduce the non-autostart Heun. Now
we will show how the same equations can be mathematically deduced.
This deduction is particularly interesting because it links the ideas of adjustment of
curve, of numerical integration and of the ordinary differential equations. The exercise is also useful because
provide a simple procedure for developing multi-step methods
higher order and estimate its errors.

The deduction is based on solving the general ODE:

Portillo Contreras Misael ISC Group 1 Page 40

 Unit VI 'Solution of Differential Equations'

This equation can be solved by multiplying both sides by integrating between

the limits :

The left side can be integrated and evaluated using the fundamental theorem:

ec. 4

The equation represents a solution to the ODE if the integral can be evaluated. It is
to say, provides a means to calculate a new value of the variable
dependent based on a previous value of and the differential equation.
The formulas of numerical integration provide a way to do this.
evaluation. For example, the trapezoidal rule can be used to evaluate the integral,
like in:

ec. 5

Where it is the step size. When substituting equation ec.5 into the
equation ec.4 is given:

Portillo Contreras Misael ISC Group 1 Page 41

 Unit VI 'Solution of Differential Equations'

What is the correction equation for the Heun method? This is based on the
Trapezoidal rule, the truncation error can be taken directly from the table.
2:

ec.6

A similar procedure can be used to deduce the predictor. For this

in this case, the limits of integration go from :

What can be integrated and rearranged to obtain:

ec.7

Now, more than using the closed formula from table 2, the first formula in
open integration of Newton-Cotes can be used to evaluate the integral as
en:

Portillo Contreras Misael ISC Group 1 Page 42

 Unit VI "Solution of Differential Equations"

ec. 8

Which is called the midpoint method. Substituting equation ec. 8 into the
equation ec.7 is obtained:

ec.9

What is the predictor for the Heun method for no self-starting. Similar to the corrector, the error
local truncation can be taken directly:

ec.10

Where the subscript p designates that this is the error of the predictor.
Thus, the predictor and the corrector for the non-self-starting Heun method have
truncation errors of the same order. In addition to updating the accuracy of
predictor, this fact has additional benefits related to the analysis of
error like himself to elaborate in the next section.
Error estimation: If the predictor and the corrector of a multipoint method are
of the same order, the local truncation error can be estimated during the course
from a calculation. This is a tremendous advantage, as it establishes a criterion for the
step size adjustment.

The local truncation error for the predictor is estimated using equation ec.9.
The estimated error can be combined with the estimate of from predictor step
to give:

Portillo Contreras Misael ISC Group 1 Page 43

 Unit VI 'Solution of Differential Equations'

ec.11

Using a similar procedure, the estimated error for the corrector can be
combine with the result of the corrector to give

Equation ec.10 can be subtracted from equation ec.11 to give:

Where is E now among y Now, if you divide the equation by 5 and you
Rearranging the result we have:

ec.12

Note that the right side of equations ec. 6 and ec. 12 are identical, with the
exception of the third derivative argument. If there is no appreciable variation
Regarding the interval in question, we can assume that the right side is equal.
and therefore, the left sides should be equivalent, as in:

ec.13

Portillo Contreras Misael ISC Group 1 Page 44

 Unit VI 'Solution of Differential Equations'

Thus, we arrive at a relationship that can be used to estimate the error of

truncation by step based on two amounts, which are routine
by-products of calculus.
Solution. In =1, the predictor gives 5.607005 and the corrector gives 6.360865. These
values can be substituted in equation ec.13 to give:

Which is well compared to the exact error:

In =2, the predictor gives 13.44346 and the trajectory gives 15.30224, which is used
to calculate:

That is also favorably compared with the error

exactly,

Portillo Contreras Misael ISC Group 1 Page 45

 Unit VI 'Solution of Differential Equations'

Conclusion

Through the previously seen methods, I understand that they are of great help.
in the application of different areas such as electronics and bacteriology
such as probability among others. Therefore, we can say that.
Numerical methods and their applications in their different methods are greatly
help in the life of the engineer as well as the graduate.

Portillo Contreras Misael ISC Group 1 Page 46

 Unit VI 'Solution of Differential Equations'

Bibliography
Numerical Methods for Engineers; Chapra Steven C. and Canalé Raymond
5th Ed.

Information about: Chapters 25, 26, 28

Pages: 713 to 846

Portillo Contreras Misael ISC Group 1 Page 47