EA3: Systems Dynamics
Better algorithms for numerical integration
Sridhar Krishnaswamy
1
Better Algorithms for Numerical Integration
Let us consider a general differential equation involving some function y(t):
DE:
dy
dt
=
f
(
y
)
(1)
along with the initial conditions IC:
y(0) = y
0
.
Here the
rate
_
function
f(y) is given to us as a function or a graph or a table of values, and
the initial value y
0
is given to us as a number.
What we want to do is calculate the variation
of y with time t;,ie obtain y(t) for subsequent times.
If you have been following the
forward Euler method that we have used so far, what we have been doing is to estimate the
value of y at the next instant by adding a correction to the current value of y.
In the
forward Euler method, the correction term is obtained rather simplistically.
We assume that
the new value of y after a time
∆
t can be obtained by simply multiplying the rate at which y
is changing
now
by the time interval
∆
t.
This will be exactly correct only if over the time
interval
∆
t, the rate at which y is changing stays constant. If not, then clearly we will be
either over estimating or underestimating the new value of y (see fig. 1).
This is why we
need to keep
∆
t very small.
But this can become computationally tiresome. Also, it did not
work well for the springblock oscillator no matter how small a time step we chose.
Can
we do better? Well to find the answer to this, first let us ask why we chose the
current
rate
of change of y (ie the current slope of the curve y(t))?
For one thing, that is all we really
know.
But can we use the current rate to get a preliminary estimate for the next y, and then
figure out a better value for the rate at which y is changing? We can try!
Forward Euler method:
Suppose at step i we have a current_value y
i
.
(I will use y
i
to mean
y(t
i
) etc as a convenient shorthand notation.) We feed the current_value to the rate_function
to get the current_rate k
1
= f( y
i
) . And then we estimate the next_value y
i+1
through:
k
1
=
f
(
y
i
)
… current_rate using current_value in rate_function
y
i
+
1
=
y
i
+
dt
*
k
1
… next_value=current_value+timestep*current_rate
{
Always read equations using meaningful words.
}
This preview has intentionally blurred sections. Sign up to view the full version.
View Full Document
EA3: Systems Dynamics
Better algorithms for numerical integration
Sridhar Krishnaswamy
2
Average Rate method:
Now what if we use the current_rate k
1
to get a
preliminary
estimate
for the next_value y
i+1
as in the forward Euler method, but then use this preliminary
next_value to get an estimate for the rate of change of y at the
other
end; ie obtain the value
of the end_rate f(y
i+1
) = f( y
i
+ dt*k
1
). We can then make a new and better estimate for y
i+1
by using the
average
rate at which y changes over this interval dt.
That is:
k
1
=
f
(
y
i
)
…current_rate
y
i
+
1
(
approx
)
=
y
i
+
dt
*
k
1
…first approximation for next_value
k
2
=
f
(
y
i
+
1
(
approx
)
)
…estimated end_rate using approximate next_value
y
i
+
1
=
y
i
+
dt
*
k
1
+
k
2
2
…better next_value using average_rate.
This is the end of the preview.
Sign up
to
access the rest of the document.
 Fall '08
 KRISHNASWAMY
 Numerical Analysis, Runge–Kutta methods, Numerical ordinary differential equations, Sridhar Krishnaswamy

Click to edit the document details