- Section 1.1 -
Direction Fields. On can interpret the
first order differential equation
$$\frac{dy}{dt}=f(t,y)$$
as a way to calculate the slope of the solution curve at various points.
This generates a direction field. This tells you (graphically)
what the solution curves are doing for various initial conditions.
For example, see the graphic below (click for full size image):
For the special case
$$\frac{dy}{dt}=f(y)$$
where $f(c)=0$ for c=constant, these values are called critical
points or equilibrium points and represent limiting values
of the integral curves as $t \rightarrow \infty$.
Make sure you can plot the direction field by hand, sketch solution
curves, and describe the behavior of integral curves.
- Section 1.2 -
When $f(t,y)$ is linear with constant coefficients,
that is for equations of the form
$$\frac{dy}{dt}=ay+b$$
One can solve by simple integration
$$ \int \frac{dy}{ay+b}=\int dt $$
This allows one to model simple population problems, and problems
involving Newton's Law of Cooling.
- Section 1.3 -
This is covered in the Matlab Project 1, but will not be on Exam 1.
- Section 2.1 -
For the general linear first order ordinary differential equation
$$\frac{dy}{dt}+p(t)y=g(t)$$
one can solve it by the method of integrating factors.
The integrating factor
$$I(t)=e^{\int p(t) dt}$$
Multiplying the differential equation by $I(t)$ and simplifying
gives
$$\frac{d{I(t)y(t)}}{dt}=I(t)g(t)$$
which can be integrated immediately, and then solved for $y(t)$.
Make sure you include the integration constant...
- Section 2.2 -
If the ordinary differential equation is nonlinear, and is of
the form
$$\frac{dy}{dt}=f(t)g(y)$$
it is called seperable and can be integrated
$$ \int \frac{dy}{g(y)} = \int f(t) dt $$
to find an implicit solution.
- Section 2.3 -
Many important physical problems can be modeled by first order
differential equations. We did two examples in class - interest,
and escape velocity. Make sure you review the details of these
two problems.
- Section 2.4 -
Differences between linear and nonlinear ordinary
differential equations.
| ODE |
$dy/dt+p(t)y=g(t)$
|
$dy/dt=f(t,y)$
|
Existence |
Existence and Uniqueness for any interval
on which $p(t)$ and $g(t)$ are continuous.
|
Continuity of $f$ and $\partial f/\partial y$
imply local existence and uniqueness.
|
| Type of Solution |
Eplicit (using integration factors) |
Implicit (or no analytical solution possible) |
| Numerical Solution |
Euler's method, fixed step size |
Implicit methods, variable step size |
| Graphical Solution, Direction Field |
Solution bounded if $p(t)$ and $g(t)$ are. |
Solutions can "blow up" at finite time |
Almost all differential equations modeling real physical problems
are non linear. Linear equations occur when processes are simplified,
or non linear terms are neglected. Singularities (i.e. infinite
behavior) does not occur in nature because terms that are small away
from the singularity begin to dominate near the singularity, which
necessitates a more complicted modeling process.
- Section 2.5 -
Autonomous Differential Equations. These nonlinear equations have
no explicit dependence on $t$ and are of the form
$$\frac{dy}{dx}=f(y)$$
If $y*$ is a constant and $f(y*)=$, then $y(t)=y*$ is a solution of
the differential equation (called a critical point) or
equilibrium solution. The direction fields are easy to
plot (since the slopes are constant on horizontal lines!)
One-species population models are of this form, including the
logistic equation
$$\frac{dy}{dx}=f(y)=r(1-\frac{y}{K})y$$
or logistic equation with threshold
$$\frac{dy}{dx}=f(y)=-r(1-\frac{y}{T})(1-\frac{y}{K})y$$
Equilibria ($f(y*)=0$) can be steady or unsteady.
depending on the sign of $f'(y*)$.
- Section 2.6 -
Exact Equations. If a nonlinear ode can be put in the form
$$M(x,y)dx+N(x,y)dy=0 (*)$$
and the mixed derivatives satisfy
$$M_y(x,y)=N_x(x,y)$$
then a function $\psi(x,y)$ exists with $\psi_x=M$ and $\psi_y=N$
and $\psi(x,y)=C$ is an implicit solution of (*). $\psi$ can be found
by integration $\int M(x,y)dx$ and $\int N(x,y)dy$.
- Section 3.1 -
Systems of two linear algebraic equations.
The system of equations
$$ a_{11}x_1+a_{12}x_2=b_1 $$
$$ a_{21}x_1+a_{22}x_2=b_2 $$
can be written as the matrix equation
`[[a_{11},a_{12}],[a_{21},a_{22}]]`
`[[x_1],[x_2]]`=`[[b_1],[b_2]]`
or Ax=b. If $\det(A)\ne 0$, then the inverse of $A$ exists,
and there is a unique solution
Solutions for a 2x2 system can be found by either elimination,
or row reduction. Cramer's rule can also be used.
An eigenvalue-eigenvector pair satisfies the relation
$$A x = \lambda x$$
which has non-trivial solutions if and only if $det(A-\lambda I)=0$.
This is called the characteristic equation, and has $N$ roots,
if $A$ is an $N\times N$ matrix. They may be real, complex, simple
or multiple roots.
- Section 3.2 -
Systems of Two First Order Linear Differential equations.
Such systems can be written as
$$ \frac{dx}{dt} = P(t) x + g(t) $$
where $P(t)$ is a $2\times 2$ matrix, and $x(t)$ and $g(t)$ are $2\times 1$
column vector. $x_1(t)$ and $x_2(t)$ are called state
variables, and the $x_1-x_2$ plane is called the phase
plane. A parametric solution curve $(x_1(t),x_2(t))$ is called a
trajectory or orbit.
If $g(t)=0$ the system is said to be homogeneous
There are several tools that are very helpful in analyzing behavior
- Direction Field - a vector field where the direction at each point
is given by the slope of the solution curve through that point.
- Phase Plane - parametric solution curves $(x_1(t),x_2(t))$ or
sometimes $(x_1(t),x_2(t),x_3(t))$.
- Component Plot - individual plots, e.g. $x_1(t)$ vs $t$.
- Section 3.3 -
Homogeneous Linear systems with constant coefficients.
Such systems can be written as
$$ \frac{dx}{dt} = [A] x(t) $$
and have solutions of exponential form $x(t)=\vec u e^{\lambda t}$
where $\vec u$ is a vector, and $[A]$ is a constant, square matrix.
Substituting this form into the differential equation, we have
$$ \lambda \vec u e^{\lambda t} = [A] \vec u e^{\lambda t} $$
or
$$ e^{\lambda t} \{ [A] \vec u - \lambda \vec u \} = 0 $$
which can only happen if
$$ [A] \vec u - \lambda \vec u $$
which means that $\vec u$ is an eigenvector of $[A]$ and
$\lambda$ is the associated eigenvalue.
The general solution can be written as
$$ x(t) = c_1 x_1(t) + c_2 x_2(t) $$
and the initial conditions yield the constraint
$$ x_0 = x(0) = c_1 x_1(0) + c_2 x_2(0) $$.
This can be written as a matrix equation
`[\vec x_1(0)|\vec x_2(0)][[c_1],[c_2]]=x_0`
This has a unique solution as long as the determinant of the
matrix `[\vec x_1(0)|\vec x_2(0)]` is non zero. This determinant
is called the Wronskiian of the vectors $\vec x_1, \vec x_2$.
For a $2\times 2$ matrix, there are at most two eigenvalues, which can
only occur as one of three cases
- real, distinct values
- real, repeated values
- complex conjugate pair
- Section 3.4 -
Complex eigenvalues.
In this case, the eigenvalues must occur as a conjugate pair
$\lambda = \mu \pm i \nu$. Using the identity,
$$ e^{\lambda t}=e^{(\mu \pm i \nu)t} =e^{\mu t}e^{\pm i \nu t} = e^{\mu t} (\cos(\nu t) \pm i \sin(\nu t)) $$
the eigenfunctions are $e^{\mu t} \sin(\nu t)$ and
$e^{\mu t} \cos(\nu t)$. The resulting phase plane portrain is a
periodic orbit (if $\mu=0$), and inward spiral (if $\mu<0$ ) and
outward spiral (if $\mu>0$ ).
- Section 3.5 -
Repeated Eigenvalues.
In this case, if $\lambda_1=\lambda_2=\lambda$ then the two eigenvectors
have the t-dependence $e^{\lambda t}$ and $te^{\lambda t}$
respectively. This results in an improper node.
Review table 3.5.1 on page 185, as well as figure 3.5.10 on page 187.
- Section 7.1 -
Autonomous systems and stability.
A system of the form
$$ \frac{dx}{dt} = f(x(t)) $$
which has no explicit time dependence is called autonomous.
Critical points are constant vectors $x_0$ satisfying $f(x_0)=0$.
the critical point is stable if every orbit $\phi(t)$
starting in the neighborhood of the critical point stays in a neighborhood.
$$ \| \phi(0) - x_0 \| < \delta => \| \phi(t) - x_0 \| < \epsilon $$
If $\lim \phi(t)=x_0$ then the critical point is asypmtotically
stable.
Nonlinear systems may posess more than one critical point, possibly even
an infinite number of critical points (e.g. nonlinear pendulum, page
477).
Curves which separate regions of differing behavior in the phase
plane are called separatrices (singular separatrix). In
general separatrices connect critical points, and divide the phase plane
into disjoint regions.
- Section 7.2 -
Almost Linear Systems.
Given an autonomous system of the form
$$ \frac{dx}{dt} = f(x) $$
We can use Taylor's theorem to expand $f(x)$ in the neighborhood of a
critical point $x_0$ (where $f(x_0)=0$). We have
$$ f(x) = f(x_0) + [J(x)](x-x_0) + ... $$
where $[J]$ is the Jacobian matrix of derivatives
`J(x) = [[ \frac{\partial f_1}{\partial x_1}, \frac{\partial
f_1}{\partial x_2}],[ \frac{\partial f_2}{\partial x_1}, \frac{\partial
f_2}{\partial x_2}]]`
If $\det J(x_0) \neq 0$, then the critical point is isolated. If
the jacobian is smooth enough, then the higher order terms are small
compared to the linear term, and the phase plane will be a small
perturbation of the linear behavior.
Adding a small term $g(t)$ to a linear system $dx/dt = Ax$
(with $det A \neq 0$) causes the critical point to move
continuously in the phase plane, and the eigenvalues to
vary continuously. Each of the cases remains the same, except for the
case of multiple eigenvalues which may split into unequal eigenvalues
(see table 7.2.2).
- Section 7.3 -
Competing species.
Systems of this type are quadratic, and consist of products of linear
terms.
$$ dx/dt = x(\epsilon_1 - \sigma_1 x - \alpha_1 y) $$
$$ dy/dt = y(\epsilon_2 - \sigma_2 y - \alpha_2 x) $$
There are at most 4 critical points, and the behavior can be completely
determined by calculating the jacobian (and eigenvalues) in terms of the
parameters.
- Section 7.4 -
Predator Prey.
These equations are a special case of the competing species equations,
with `\sigma_1=\sigma_2=0`.
$$ dx/dt = x(\epsilon_1 - \alpha_1 y) $$
$$ dy/dt = y(\epsilon_2 - \alpha_2 x) $$
Two of the critical points vanish, leaving to, one at the origin (0,0)
and another generating periodic orbits about the point
$(\epsilon_2/\alpha_2,\epsilon_1/\alpha_1)$.
- Section 4.1 -
Second Order Linear Equations
These equations are of the form $y''(t)=f(t,y,y')$ where $f$ is a linear
function of $y$, that is
$$ y''(t)+p(t)y'(t)+q(t)y(t)=g(t) $$
This can be transformed into a $2\times 2$ system
using the substitutions $x_1=y$ and $x_2=y'$.
The following three physical examples all reduce to a second order
linear equation:
- `m y''(t) + \gamma y'(t) + k y(t) = F(t)` [spring-mass system]
- `mL \theta''(t) + c \theta'(t) + mg \sin(\theta) = 0` [nonlinear
pendulum]
- `LI''(t) + R I'(t) + \frac{1}{C} I(t) = E(t)` [RLC circuit]
- Section 4.2 -
Theory of Second Order Linear Equations:
These are of the form
$$ y''(t)+p(t)y'(t)+q(t)y(t)=g(t) $$
If the coefficients $p(t),q(t),g(t)$ are all continuous on an interval
time $I=[a,b]$ (containing the initial time $t_0$) then a unique
solution of the initial value problem will exist, as long as the
coefficients remain continuous. One way to think about this is to state
that the solution will continue to exist as long as nothing "goes wrong"
with the coefficients in the problem.
Specializing to the homogeneous version of the ode
$$ y''(t)+p(t)y'(t)+q(t)y(t)=0 (1) $$
we have the Principle of superposition. If $y_1(t),y_2(t)$ are
two solutions of (1), then any linear combination is also
$$y_h(t)=c_1 y_1(t) + c_2 y_2(t) $$
This is the general solution of the homogeneous second order equation.
If the Wronskian
`W(y_1,y_2) = \det [[y_1(t_0),y_2(t_0)],[y_1'(t_0),y_2'(t_0)]]`
does not vanish, then there is a unique solution of the initial value
problem.
- Section 4.3 -
Linear Homogenous Equations with Constant Coefficients
Specializing to this case
$$ y''(t)+by'(t)+cy(t)=0 $$
we observe that solutions must be of the form $y(x)=e^{\lambda t}$.
Substituting this, we conclude that $\lambda$ must satisfy the
characteristic polynomial equation $Z(\lambda)=\lambda^2+b \lambda +
c=0$. The same values of $\lambda$ result from writing the
second order equation as a first order system.
As before, there are 3 cases to consider
- real, distinct roots
- real, repeated root
- complex conjugate root
To graph the behavior in the phase plane, it is easiest to consider the
equivalent $2\times 2$ system.
- Section 4.4 -
Characteristic Equations with Complex Roots
For linear constant coefficient differential equations of the form
a y''(t) + b y'(t) + c y(t) = 0
the solutions are of exponential form $e^{\lambda t}$,
where $\lambda$ is a root of the characteristic equation $a \lambda^2
+b\lambda +c=0$.
This has complex roots if the discriminant $b^2-4ac<0$. It will be
purely imaginary if, in addition, $b=0$.
To write the solutions in real terms, we use the Euler formula
$$exp((a+bi)t)=exp(at)(\cos(bt)+i \sin(bt))$$
We can take the real and imaginary parts as homongenous solutions
$$y_1(t)=exp(at)\cos(bt), y_2(t)=exp(at)\sin(bt)$$
Another equation that can have complex roots is the Caucy-Euler equation
$$ a t^2 y''(t) + b t y'(t) + c y(t) = 0 $$
This ode has solutions of the form $t^\lambda$. After substituting this
into the differential equation, we see that \lambda must satisfy the
equation
$$ a\lambda(\lambda-1)+b\lambda+c=0 $$
If $\lambda$ is complex, we write
$$ t^\lambda = t^{a+ib} = e^{\ln(t)(a+bi)} $$
$$ = e^{a\ln(t)}e^{ib\ln(t)} = t^a (cos(b\ln(t))+i\sin(b\ln(t))) $$
- Section 4.6 -
Method of Undetermined Coefficients.
For linear, constant coefficient, second order odes, if the forcing term
is polynomial, trigonometric, or exponential (or products thereof), one
can find particular solutions of special form.
If
$$ a y''(t) + b y'(t) + c y(t) = g(t)$$
then the following cases hold, in general:
- If $g(t)=e^{at}$ then $y_p(t)=Ae^{at}$
- If $g(t)=\sin(at),\cos(at)$ then $y_p(t)=A\sin(at)+B\cos(at)$
- If $g(t)=P_N(t)$ (a polynomial of degree N) then
$y_p(t)=Q_N(t)$ (a polynomial of degree N).
- if $g(t)=P_N(t)e^{at}\sin(bt)$ then
$y_p(t)=Q_N(t)e^{at}\sin(bt)+R_N(t)e^{at}\cos(bt)$, where $Q_N$ and
$R_N$ are polynomials of degree N as well.
NOTE: If the form of the particular solution includes
a homogeneous equation, one must multiply for an appropriate power
of $t$.
- Section 4.8 -
Variation of Parameters
For a second order linear differential equation, or a system of 2 first
order linear differential equations, knowledge of the two fundamental
(homogeneous) solutions allows one to construct a particular solution.
In the 2x2 matrix case
$$ x'(t) = [P(t)] x(t)+g(t) $$
where $x_1(t)$ and $x_2(t)$ are fundamental solutions of
$$ x'(t) = [P(t)] x(t) $$
and $[X]=[x_1|x_2]$ is the 2x2 matrix formed from the homogeneous
solutions, then the particular solution is formed from
$$[X(t)]u(t)$$
where the vector $u(t)$ satisfies the equation
$$u'(t)=[X(t)]^{-1}g(t)$$
Integrating this last equation, and multiplying by $[X(t)}$ gives the
particular solution.
In the case of a single second order equation
$$ y''(t)+p(t)y'(t)+q(t)y(t)=g(t)$$,
with the homogeneous solutions $y_1(t)$ and $y_2(t)$,
a particular solution will be of the form
$$y_p(t)=u_1(t)y_1(t)+u_2(t)y_2(t)$$
where $u_1(t)$ and $u_2(t)$ satisfy the system
$$ y_1(t)u_1'(t)+y_2(t)u_2'(t)=0 $$
$$ y_1'(t)u_1'(t)+y_2'(t)u_2'(t)=g(t)$$
Note: The equation must be in standard form, in order to
calculate the correct $g(t)$.
- Section 5.1 -
Definition of the Laplace Transform
The definition of the Laplace Transform is
$$F(s)=\int_0^\infty e^{-st}f(t)dt $$
This is well defined for any piecewise continuous functions $f(t)$
which grow at most linearl-exponentially, that is
$$|f(t)|\leq Ke^{at}$$
for some $K>0$. The laplace transform is linear. The Laplace transform
is a mapping from the physical domain (time) to the transform domain.
The laplace transform of any piecewise continuous, linear-exponentially
bounded function, must go to zero as $s \rightarrow \infty$.
- Section 5.2 -
Properties of the Laplace Transform
The Laplace transform has many properties, the major one being that it
maps differential equations in the physical domain to algebraic
equations in the transform domain.
- $L\{e^{ct}f(t)\}=F(s-c)$, i.e. exponentiation -> translation
- $L\{f'(t)\}=sL\{f(t)\}-f(0)$, i.e. differentiaion -> multiplication
by $s$ and incorporation of initial conditions.
- $L\{t^n f(t) \} = (-1)^n \frac{d}{ds} F(s)$, i.e. powers ->
derivatives
- Section 5.3 -
The Inverse Laplace Transform
If $F(s)=L\{f(t)\}$ then $f(t)=L^{-1}[F(s)]$. $f(t)$ is the
inverse laplace transform of $F(s)$. Inverting the laplace transform is
equivalent to recognizing which $f(t)$ gives rise to $F(s)$.
Because we cannot pre-compute every possible transform, only the major
ones are listed [Table 5.3.1 on page 325 of the text].
The major techniques used to simplfy tranforms are
(roughly in order used)
- Method of partial fractions
- completing the square
- differentiation in the transform plane
The method of partial fractions is outlined on page 328 of the text.
Completing the square ensures that the constants found in the partial
fraction decomposition match the constants appearing in the particular
solution.
- Section 5.4 -
Solving Differential Equations with Laplace Transforms
Solving a linear differential equation by the method of Laplace
transforms amounts to performing the following:
- Take the laplace transform of both side
- Incorporate initial conditions
- Solve for laplace transform of solution, $Y(s)$
- Simplify $Y(s)$ using method of partial fractions, completing the
square, etc.
- Find the inverse transform of each term.
- Section 5.5 -
Discontinuous functions, Periodic functions
If the forcing term is piecewise continuous, or periodic, special
methods have been developed. Shifting the forcing term by an amount
$t=c$ gives a forcing term of the form $u_c(t)f(t-c)$ where
$u_c(t)$ is the Heaviside step function
$$ u_c(t)=\{0, 0 < t < c ; 1, t\geq c \}$$
In this case
$$L\{u_c(t)f(t-c)\} = e^{cs}F(s) \}$$
Note: it is critical that the values of "c" match in
$u_c(t)f(t-c)$.
Generally, in break in smoothness (continuity of function or a
derivative) will intruduce a factor $u_c(t)$ in the solution.
Periodic functions: If $f(t+T)=f(t)$ for all $t>0$ then $f(t)$ is
periodic, with period $T$. (We assume $T$ is the smallest period).
For such functions,
$$L\{f(t)\} = \frac{\int_0^T e^{-st}f(t)dt}{1-e^{-sT}}$$
- Section 5.6 -
Differential Equations with discontinuous forcing functions
when solving differential equations with discontinuous forcing
functions, the critical first step is to write the forcing term
correctly! Some familiarity with infinite (geometric) series is assumed
(from Calc II, III).
- Section 5.7 -
Impulse functions
A unit force, applied over an arbitrarily small time interval, is called
an impulse function, \delta(t-t_0). It has the property that
\int_{-\infty}^{\infty} f(t) \delta(t-t_0) dt = f(t_0)
In particular, it has the transform
$ L[\delta(t-t_0)] = \int_0^{\infty} e^{-st} \delta(t-t_0) dt
= e^{-st_0} $, and $L[\delta(t)]=1$.
When an impulse function \delta(t-c) is present as a forcing term,
the transform of the solution will involve $e^{-cs}$ in the numerator,
hence there will be a term u_c(t) present in the solution.
- Section 5.8 -
Convolution Integrals.
The convolution operator f \star g \equiv \int_0^t f(t-s) g(s) ds is a linear
operator. The most important property is
$ L[f \star g] = L[f] L[g] = F(s)G(s) $
and
$ L^{-1}[F(s)G(s)] = f \star g $
When solving $ay''+by'+cy=g(t),y(0)=A,y'(0)=B$, we get
$
Y(x) = \frac{asA+aB+bA+G(s)}{as^2+bs+c}
=A\frac{as+b}{as^2+bs+c} + B\frac{a}{as^2+bs+c}
+H(s)G(s)
=A Y_1(s) + B Y_2(s) + H(s)*G(s)
$
where
$H(S)=1/(as^2+bs+c)$
This gives us
y(t) = A y_1(t) + B y_2(t) + H(t) \star g(t)
y_1(t) is a homogeneous solution satisfying y_1(0)=1 and y_1'(0)=0,
y_2(t) is a homogeneous solution satisfying y_2(0)=0 and y_2'(0)=1.
The particular solution y_p(t) = H(t) \star g(t) is a convolution,
and satisfies homogeneous initial conditions y_p(0)=0, y_p'(0)=0.
This is a general result for linear constant coefficient
differential equations.
- Section 2.7 -
Accuracy of Numerical Methods.
Euler's method for the initial value problem
y'(t)=f(t,y), y(t_0)=y_0
can be derived in many ways. We consider three.
- I. Tangent line approximation.
The slope of the tangent to the solution curve at a point t_n is given by
m=f(t_n,y(t_n)). Therefore, at the next point t_{n+1}, we have the
tangent line approximation y(t_{n+1}) \approx y(t_n)+f(t_n,y(t_n))*
(t_{n+1}-t_n). This gives us Euler's method:
y_{n+1} = y_n + f(t_n,y_n)*h
t_{n+1}=t_n+h
starting with the initial condition (t_0,y_0).
II. Integral formulation.
Since any solution of the initial value problem satisfies
the identity
\phi'(t) = f(t,\phi(t))
integrating both sides, and using the left hand rule, we get
\phi(t_{n+1}) - \phi(t_n) = \int_{t_n}^{t_{n+1}} f(s,\phi(s)) ds
\approx f(t_n,\phi(t_n))*(t_{n+1}-t_n).
This yields
\phi(t_{n+1}) = \phi(t_n) + f(t_n,\phi(t_n))*(t_{n+1}-t_n)
which is Euler's method:
y_{n+1} = y_n + f(t_n,y_n)*h
t_{n+1}=t_n+h
III. Taylor approximation, with remainder.
A generalization of the tangent line approximating, using polynomials,
is called Taylor series. At the point t_n, we have
\phi(t_{n+1}) = \phi(t_n) + \phi'(t_n)h + \phi''(t_n)\frac{h^2}{2}
+ ... + \phi^{(n)}(t_n) \frac{h^n}{n!} + ...
There is a version, finite terms with remainder.
\phi(t_{n+1}) = \phi(t_n) + \phi'(t_n)h + \phi''(t_n)\frac{h^2}{2}
+ ... + \phi^{(n)}(\xi_n) \frac{h^n}{n!}
where \xi_n is a point on the interval (t_n,t_{n+1}).
Setting n=2 in this formula, yields a version of Euler's method
\phi(t_{n+1}) = \phi(t_n) + f(t_n,\phi(t_n))h + \phi''(\xi_n)\frac{h^2}{2}
y_{n+1} = y_n + f(t_n,y_n)*h
Subracting these two equations, and setting the error E_n
= \phi(t_n) - y_n, we have
E_{n+1} = E_n + [f(t_n,\phi(t_n))-f(t_n,y_n)]h + \phi''(\xi_n)\frac{h^2}{2}
= E_n + \frac{\partial f}{\partial y}(t_n,z_n)*E_n*h+
\phi''(\xi_n)\frac{h^2}{2}
|E_{n+1}| \leq |E_n| + M*|E_n|*h + C*\frac{h^2}{2}
where M=max|\frac{\partial f}{\partial y}|, and C=max|\phi''()|.
This yeilds the fundamental inequality for errors
|E_{n+1}| \leq (1+Mh)*|E_n| + C*\frac{h^2}{2}
The term |E_n| is called the global truncation error. If we assume
that E_n=0, we get a single step error
e \leq C*\frac{h^2}{2}
where e is the local truncation error (ie single step, with
perfect starting values).
Euler's method has local truncation error of O(h^2), and a global
truncation error of O(h). To improve accuracy, we try to increase the
poewr of h. This can be done by using a higher order Taylor
series, or a better inegral approximation (see next section).
- Section 2.8 -
Improved Euler and Runge Kutta method.
If we use the trapezoidal method for approximating the integral:
\phi(t_{n+1}) - \phi(t_n) = \int_{t_n}^{t_{n+1}} f(s,\phi(s)) ds
\approx (1/2)[f(t_n,\phi(t_n))+f(t_{n+1},\phi(t_{n+1})]*(t_{n+1}-t_n).
we get an implicit numerical integration scheme
y_{n+1} = y_n + (h/2)*(f(t_n,y_n)+f(t_{n+1},y_{n+1}))
the update y_{n+1} is on both sides of the equation, so we either must
solve a nonlinear equation (computationally expensive), or we
approximate y_{n+1} on the right hand side. This leads to
$\tilde y_{n+1} = y_n + h*(f(t_n,y_n))$
$ y_{n+1} = y_n + (h/2)*(f(t_n,y_n)+f(t_{n+1},\tilde y_{n+1})) $
This is called a predictor-corrector
method, where the initial update
is given by a lower order scheme (called the predictor) and followed
by a more accurate scheme (corrector).
This is called the improved euler method, and is third order accurate
(locally) and second order accurate globally. This process can be
continued, using better integration methods.
An alternative is to develop a more accurate integrator, called
the Runge Kutta method. It uses a sequence of updates to the slopes
according to
k_{n1} = f(t_n,y_n) [Euler]
k_{n2} = f(t_n+h/2,y_n+k_{n1}*h/2)
k_{n3} = f(t_n+h/2,y_n+k_{n2}*h/2)
k_{n4} = f(t_n+h,y_n+k_{n3}*h)
k_n = (1/6)*(k_{n1}+2k_{n2}+2k_{n3}+k_{n4})
y_{n+1} = y_n + h*k_n
t_{n+1} = t_n + h
This method is O(h^5) locally, and O(h^4) globally. It is one of the
most common numerical schemes for nonlinear differential equations.