Numerical Differential Equations

Table of Contents

Notation

  • numerical_differential_equations_9dfb630910b0b3d2a61cb5f1673b9adaa08fe150.png almost always referes to the Lipschitz constant
  • numerical_differential_equations_ded6419ca754335b5d0d87283cd3621fd98424a6.png refers to LHS of a first-order system
  • numerical_differential_equations_3a68c3d2288d0943ca05fe542aa996519de6ff2a.png denotes a "sampled" value of numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png at numerical_differential_equations_e9e5146ccb36f651f748ef01a5d643bf0c07b17d.png, using our approximation to numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png
  • numerical_differential_equations_5bf24d8ad956ea1754976973a976bb7378ea20f0.png refers to a flow-map, with numerical_differential_equations_4e21bb983a723f0e8d1600c96cd9460643f4a766.png assuming numerical_differential_equations_afeafb2dace20d4696684bc27cf0f53550093371.png (as we can always do)
  • numerical_differential_equations_c64efccc7f804dbbbd616861932b89862c8cb925.png refers to the discrete approx. of the flow-map
  • numerical_differential_equations_b689cba8d7566f6adaf605a844e193a27e155078.png refers to the spatial domain, i.e. numerical_differential_equations_56923cdb42130fbc5f38c0a7381b9760d37ccfd7.png for some numerical_differential_equations_eb5d809ed7c492fae7d4927a6fc9a5e22f9b3831.png
  • numerical_differential_equations_e0a1be52c4af9683961ac1fd6da1b4f16cbbff44.png denote the space of real polynomials of degree numerical_differential_equations_54c94faafad0ff0c733bf57e1cea31c9391c2fee.png
  • numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png denotes the i-th abscissa point

Definitions

Continuous mechanics

A system of differential equations for which numerical_differential_equations_684f34a94654106dd5afb679c9688b420cda9800.png can be written as a function of numerical_differential_equations_a3a7f43f807b9e381fc50e0fab140c0df0a03e17.png only is an autonomous-differential-equation (often referred to as continuous mechanics):

numerical_differential_equations_9c27fcf6f0a9ee9697c9ec3243ef652cbac307f8.png

Euler's method

We start by taking the Taylor expansion of the solution numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png of the system:

numerical_differential_equations_eb731462c969b26a59b0703fce6cbb11893ae2a1.png

for some numerical_differential_equations_9bf71855f29538fb11330e947b9859d81113167a.png.

We have the ODE

numerical_differential_equations_3ea466cd3975df5c40554799d380b8189ea7646c.png

and the steps are defined as

numerical_differential_equations_549c600a18974ceb6599fd80b1ef5c30fc45cb15.png

for

numerical_differential_equations_abbf3768e80360804507e5c5365d49046469c6a6.png

where numerical_differential_equations_28a6c3e2b6553638465490c5c46622598b5b766a.png denotes our approximation for numerical_differential_equations_28d5a996b3922fba856f77de053129c309f8290a.png.

What we're doing in Euler's method is taking the linear approximation to numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png at each point numerical_differential_equations_e9e5146ccb36f651f748ef01a5d643bf0c07b17d.png, but in an interative manner, using our approximation for numerical_differential_equations_a3a7f43f807b9e381fc50e0fab140c0df0a03e17.png, at the previous step numerical_differential_equations_f5974766b622a636862b5f99c8be0e3ea23ee4ea.png to evaluate the point numerical_differential_equations_f5974766b622a636862b5f99c8be0e3ea23ee4ea.png.

This is because we can evaluate numerical_differential_equations_be87a63a3c9d3f0c778664e325625665e0ba90ea.png, since numerical_differential_equations_77a96d97cccd8732a088e0a643ad2a099bfdb193.png but we of course DON'T know numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png.

y = @(t, y0) (y0 * exp(t) ./ (1 - y0 + y0*exp(t))
y0 = 0.2; T=0:0.5:5; plot(T, y(T, t0));

Butcher table

In matrix notation:

numerical_differential_equations_47c63bea8bff1f08cc1eadc423f4c2efe51a5c6c.png

or elementwise:

numerical_differential_equations_771bbdba976592a01aeb39c7cd11ef4f9a9ac4be.png

Example of RK method with 4 stages:

numerical_differential_equations_fec850f044fc55a35bab747d503608f8fae74f09.png

Flow map

Consider the IVP:

numerical_differential_equations_09d048f8e5aed7a6be0ae26a1a6af5451fc0cb8b.png

which is assumed to have a unique solution on numerical_differential_equations_1e317d7bcb9876e4394628a1f5d060094895e80c.png.

The flow map numerical_differential_equations_8c0703d4a8de77b812dda8f64341c716bfc87f61.png is then an approximation to a specific solution (a trajectory), i.e. starting at some numerical_differential_equations_0428ec35c043ef5ebae3dea7001a167a9479ecc9.png such that numerical_differential_equations_60abf10db8f93ae14038fecf27f309b802cd069e.png and numerical_differential_equations_7d5ae6189e30100fc3d7a7e9ca76afb88477993e.png, which we write

numerical_differential_equations_652db22c744c55473a17adfd8c9a1e9832749da7.png

Observe that it's NOT a discrete approximation, but is continuous.

numerical_differential_equations_ffaff49aa5aaeb5da449b4298eedfa261b5dd3c1.png

Theorems

Lipschitz convergence

Suppose a given sequence of nonnegative numbers satisfies

numerical_differential_equations_324036b7b4bb039106e73f8f29fd59722c50a3cd.png

where numerical_differential_equations_370a3fd0dccdae44011a90f2651694251a0e67c5.png are nonnegative numbers numerical_differential_equations_92dd68dd38b5c71f89979ef90d5ffb6bc78fcd52.png, then, for numerical_differential_equations_981119e89c40f9a4a29cd8d3422c6741ea356dea.png

numerical_differential_equations_bb0feca9d2ec8c9a76fe90dc5125547cb2a5cb2f.png

Local Existence / Uniqueness of Solutions

Suppose numerical_differential_equations_30fd9abf0e747323eefc3e3afe1370bed0878232.png is:

  • continuous
  • has continuous partial derivatives wrt. all components of the dependent variable numerical_differential_equations_a3a7f43f807b9e381fc50e0fab140c0df0a03e17.png in a neighborhood of the point numerical_differential_equations_b6cea1488fec77943e5af798a442a58283a4e884.png

Then there is an interval numerical_differential_equations_972f3de10d44e70a225d2aa68f7a69d6a0f523b1.png and unique function numerical_differential_equations_a3a7f43f807b9e381fc50e0fab140c0df0a03e17.png which is continuously differentiable on numerical_differential_equations_6c03018fd4253198b850bc6339e99699b567a798.png s.t. it satisfies

numerical_differential_equations_2df60634d8440580a8eec6aeec930e80ba18189d.png

Methods

Euler method

Trapezoidal rule

The Trapezoidal rule computes the average change on the interval numerical_differential_equations_ec7f3c691f3d5ab3765d62a0b073b377de2fd154.png, allowing more accurate approximation to the integral, and the update rule is defined as follows:

numerical_differential_equations_c3cdfc5270f99ae61d4448ea02264b46e7faab06.png

Heun-2 (implicit approx. to Trapezoidal rule)

The Heun-2 update rule is an implicit approximation to Trapezoidal rule, where we've replaced the explicit computation of the endpoint numerical_differential_equations_a66a26096072fe3205af9cb8af296804797469f5.png which depend on numerical_differential_equations_6ddbdcab3e6b9c5be1795e98eef069edfef5847e.png, with a function of the current estimate numerical_differential_equations_28a6c3e2b6553638465490c5c46622598b5b766a.png:

numerical_differential_equations_0ddb1e49bc7ea03b5cb28d8b3bbb7c4653be5504.png

Convergence of one-step methods

Notation

  • numerical_differential_equations_36fcfe9ff814229278ef9930b497953c2aacfc61.png
  • numerical_differential_equations_50458081061f5c8b16fad4de9489afac7ffe7574.png
  • numerical_differential_equations_9441deb9842e2645a3c72b57bc780e8ec802456c.png
  • numerical_differential_equations_f121be31db8202ae375e62053ed4614312536927.png denotes the stage values, s.t. numerical_differential_equations_bc6a09e98ca36b04933b0924097c21d65b8dd36b.png

Approximation

In this chapter we focus on the autonomous differential equation

numerical_differential_equations_ab09b44a7b4a845a7c06e727e45a2260ef9dc544.png

but the non-autonomous differential equation

numerical_differential_equations_b865832b9713e3e8d4bf4f963f9e60197d5eb638.png

can be treated using similar methods.

The global error after numerical_differential_equations_f07b932b12c91bca424fd428d383495413b5b6a9.png time steps is the difference between the discrete approx. and the exact solution, i.e.

numerical_differential_equations_508f2b797bf157c6299dab042e74804e7402b722.png

A method is said to be convergent if, for every numerical_differential_equations_b3731d37cd447bd6c31809075a6be43b3d0b04ec.png,

numerical_differential_equations_ce0ed2f608758edb217e24db2c7d5250e8677a1f.png

We have the Taylor expansion given by

numerical_differential_equations_a4c512aa41d425caa5fd0dc2aad983ac99d95b22.png

for some "perturbation" numerical_differential_equations_e10db6797deed3c9f8f1edc4735630ee6de2ac0e.png about numerical_differential_equations_eb5d809ed7c492fae7d4927a6fc9a5e22f9b3831.png.

For a scalar function, assuming numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png is numerical_differential_equations_da8e75455462b818a9b164371ccdd9c0de0995fb.png times continuously differentiable (i.e. numerical_differential_equations_3b7c19d695bb859898bcdc95d407914d4356080b.png), Taylor's theorem states that

numerical_differential_equations_b2ca57ba2acd3f8638e90109c0b45daba42c34e2.png

Thus, the norm of the last (remainder) term is bounded on numerical_differential_equations_dc02381052a0b6a5c2a7c2f846bfcb951a794854.png, and we have

numerical_differential_equations_503d7058710e7b0d4236ef2faa6bb212bf58de43.png

for some constant numerical_differential_equations_e20dcda5f035650122343c61053a7c3ad6acacaa.png, and thus we write

numerical_differential_equations_c6a7a362469b950d5cc96b4580e05e39e7a3e05c.png

in general.

For a vector function (numerical_differential_equations_83637c30a7f380547c199d77745ca90aa4e5ca73.png), Taylor's theorem still holds componentwise, but in general the mean value will be attained at a different numerical_differential_equations_bfa837ba4bef06233f12915139774c147ebffe00.png for each component.

Local error of a numerical method is the difference between the (continuous) flow-map numerical_differential_equations_8c0703d4a8de77b812dda8f64341c716bfc87f61.png and its (discrete) approx. numerical_differential_equations_23e1b2ca370f895a2dfac77cad0a62e25e479ddf.png:

numerical_differential_equations_c59b983e55503c61c1824c8a91eeca3000981ac2.png

Which measures how much error is introduced in a single timestep of size numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png.

A method is said to be consistent if and only if the local error satisfies

numerical_differential_equations_a075a8390fb719b43dd663e39f8226d4117f5e2b.png

where numerical_differential_equations_e20dcda5f035650122343c61053a7c3ad6acacaa.png is a constant that depends on numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png and its derivatives, and numerical_differential_equations_813342791ebbfd19426b72f6a4b92a237d8049f2.png.

A method is said to be stable if and only if it satisfies an h-dependent Lipschits condition on numerical_differential_equations_b689cba8d7566f6adaf605a844e193a27e155078.png (the spatial domain)

numerical_differential_equations_5bfc5ef3418236176916d60e3ddbbef84fcb3666.png

where numerical_differential_equations_bf75946b4aa9467f1228df9d5d1f4306f520f57c.png is not necessarily the Lipschitz constant for the vector field.

Given a differential equation and a generalized one-step method numerical_differential_equations_2b7ec4e4c22bbd4a8c2092e4ff5187397b12ad35.png which is consistent and stable, the global error satisfies

numerical_differential_equations_1e77dac533fbc5c479d3761d963fe3c8d39166fa.png

Let the error be

numerical_differential_equations_2266c2d0c762811db1ecd29bd0005a2c813682eb.png

so

numerical_differential_equations_c1313b803b1d742041b58ab347afba10d6aa54ba.png

Then, adding and subtracting by numerical_differential_equations_b33ed7e4b005280edcf5dd6fc8fe5959a5433a86.png,

numerical_differential_equations_f99625ef46d70634eaab83709a1e3e52cdb2c3ba.png

Then, using the fact that the method is consistent and stable, we write

numerical_differential_equations_b7ac7957cf5f2d121d678f8466b91b6e3308e261.png

Using this theorem, yields

numerical_differential_equations_bdc2981f707c14ae0509fd43cc9b7d208413645f.png

where numerical_differential_equations_75f8efa80ac85446160efa7b47d57694ad1c3053.png. Finally, since numerical_differential_equations_22e859a650c6cf44c978ef741a08a2aacc102196.png and therefore numerical_differential_equations_e07273a5ddf0f7281f8b4973632024a77ea56f4d.png, if we assume the initial condition is exact (numerical_differential_equations_490d6b4cee752d2eb67ae898b0852cd7aaa3a392.png), we get the uniform bound

numerical_differential_equations_9fd70d2b0b4b1b7129fb06fffc7a6a283a51f64e.png

which proves convergence at order numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png.

A really good way of estimating the order of convergence numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png for the maximum global error wrt. numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png, is to plot numerical_differential_equations_eaac7646cb022aa963b396c0e936a271d8f5e34f.png vs. numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png on a log-log plot since

numerical_differential_equations_925ab04885a8801a251745062a65ba1dff9e03b2.png

Thus, numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png is the slope of numerical_differential_equations_e919a2737a9264078e33d9e2b2896fee45543fb7.png.

problem-2loglog.png

Figure 1: Observe that for smaller values of numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png we have numerical_differential_equations_45b45d58a4ae4130d828ee8475091cc1a733b764.png, i.e. numerical_differential_equations_eaac7646cb022aa963b396c0e936a271d8f5e34f.png decreases wrt. numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png, while for larger values of numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png we have numerical_differential_equations_be3793ad99d01feafd8da7b7ab380df0513155ac.png, i.e. decrease in numerical_differential_equations_eaac7646cb022aa963b396c0e936a271d8f5e34f.png is not even on the order of 1.

Example: proving convergence of Euler's method

Consider a compact domain numerical_differential_equations_346c703b166f8d27443db85216d64c4c0922845e.png and suppose numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png is a smooth vector field on numerical_differential_equations_b689cba8d7566f6adaf605a844e193a27e155078.png and has a Lipschitz constant numerical_differential_equations_9dfb630910b0b3d2a61cb5f1673b9adaa08fe150.png on numerical_differential_equations_b689cba8d7566f6adaf605a844e193a27e155078.png (since numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png is smooth, we can take numerical_differential_equations_2b96c03f8a37650d97f1fffc0ba38d5aba0afd11.png).

Then, since

numerical_differential_equations_d5b3bcb7cd6f3284b894d0257608f38721f5b111.png

the numerical flow map is Lipschitz with numerical_differential_equations_39f68ac849463f90f58a54714eb53d4e4dfbe643.png.

The exact solution satisfies

numerical_differential_equations_5e324f75fac657566dfbee317240c544f8675dd3.png

Therefore the local error is

numerical_differential_equations_43a8b7722c3023802f0d8d24d6547b781d5a67f9.png

and we can apply the convergence theorem for one-step methods with numerical_differential_equations_133788d5759a43698bd20076363e9ba52506a87f.png to show that Euler's method is convergent with order numerical_differential_equations_45b45d58a4ae4130d828ee8475091cc1a733b764.png.

Construction of more general one-step methods

Usually one aims to approximate the integral between two steps:

numerical_differential_equations_eb792a0e73cf6cd8e21daee28d47aab6a0b4968d.png

A numerical method to approximate a definite integral of one independent variable is known as numerical quadrature rule.

Example: Trapezoidal rule

The approximation rule becomes:

numerical_differential_equations_ee1ba249f522a966e9071df0cc2648737baac292.png

which results in the trapezoidal rule numerical method:

numerical_differential_equations_c71d2c7250bb950ac73976d63749d2cf08a62492.png

This is an implicit method, since we need to solve for numerical_differential_equations_6ddbdcab3e6b9c5be1795e98eef069edfef5847e.png to obtain the step.

Polynomial interpolation

Given a set of abscissa points numerical_differential_equations_0332c31b33e898a40b20b56cd371240af6356ae2.png and the corresponding data numerical_differential_equations_edf9042581e4c2a56d558da9de0c045515030e9e.png, there exists a unique polynomial numerical_differential_equations_251efcee8462f78b41e4654661a12e15ff27425f.png satisfying

numerical_differential_equations_87514fbf7bc1aa5fdcfc2ed13f2788d324c6873f.png

called the interpolating polynomial.

The Lagrange interpolating polynomials numerical_differential_equations_d64a53d0453a8368cf1cfb9f2408fa6c47443191.png for a set of abscissae are defined by

numerical_differential_equations_847baac8b29a96b0d47212c65fdb77d16362b8f6.png

Observe that numerical_differential_equations_3c351480f1a78cb5c7b6a4d8dd95f2ae7856fea5.png then has the property

numerical_differential_equations_76345a482b1a00eb553aae5e6f1d659561c490ee.png

Thus, the set of Lagrange interpolating polynomials form a basis for numerical_differential_equations_e0a1be52c4af9683961ac1fd6da1b4f16cbbff44.png. In this basis, the interpolating polynomial numerical_differential_equations_f8bec95b572243552399d6897a7cf76ffb3f9118.png assumes the simple form

numerical_differential_equations_6ffcc523de3c70e13f2d89133119222902de382e.png

We also know from Stone-Weierstrass Theorem that the polynomials are dense in numerical_differential_equations_e56a35cc80724c0d5f285ca269e95eee31cb7064.png for some compact metric space numerical_differential_equations_0207be880056b9a69e22e729dd37bced29cd174a.png (e.g. numerical_differential_equations_9f3373f62a3202ea9c8c771e8432b7f069ccf631.png).

Also, if it's not 100% clear why numerical_differential_equations_eac4f26e804bbd0d7edf1a1153db5e20be140433.png form a basis, simply observe that if numerical_differential_equations_30729d4c5ad30750dea6f9c40baf9508e908e80e.png then

numerical_differential_equations_56bb656d2cc335feaa61fb383784ef3ec31e5dad.png

and if numerical_differential_equations_a41c04a3f3dafda05d9e0c209b6ed8613c844475.png, then we'll multiply by a factor of numerical_differential_equations_8c5ea6850e8ec6c0c3437869dead68abd74771f8.png, hence

numerical_differential_equations_a2dbc4f2db5e757f14b851f77671d8cd598866dc.png

as wanted. Thus we have numerical_differential_equations_e841e92007e2b0d36c4f5c972a96001fd1289af9.png lin. indep. elements, i.e. a basis.

Numerical quadrature

We approximate the definite integral of numerical_differential_equations_bed561b338628a088ae69301d4bba9fd83a70bd2.png based on the interval numerical_differential_equations_14ad9fc7ad19fd0d1dd65feedc0d5b8c171b6486.png by exactly integrating the interpolating polynomial of order numerical_differential_equations_df5b12b9e164239a2d8e0837f798409ebcd29b6c.png based on numerical_differential_equations_9493b9776211ef51195c1379ad5f5540a04e566e.png points numerical_differential_equations_71272bc3ecd8a31d4b7a8adf145999d6391d71fb.png.

The points numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png are known as quadrature points.

We make the following observation:

numerical_differential_equations_a394dffd1ea69e3a6b7288d13c6502b137508ff1.png

We're interested in approximating the integral of numerical_differential_equations_4b13e3e4e74f6cb24c1f3a97899fc69c442035ad.png on the interval numerical_differential_equations_be9cc608ab30d84d7f8d8595d21933a78ae69595.png, thus

numerical_differential_equations_c72362284ae9da5861af7cd8c62000af5a3b3cfa.png

using the change of variables numerical_differential_equations_b30bfa53070a8a3f43e094b82a3853a7895f6aad.png. With the polynomial approx. we then get

numerical_differential_equations_9331b18dbce988a0dbb21e8da991a8bf5775f2b5.png

where we've let

numerical_differential_equations_7630ab46d810c28e7d9f6bb8721e64df29380fe1.png

The integral-approximation above is called a quadrature formula.

By construction, a quadrature formula using numerical_differential_equations_9493b9776211ef51195c1379ad5f5540a04e566e.png distinct abscissa points will exactly integrate any polynomial in numerical_differential_equations_fda6846e9d9d4e86a0249d4a83170869f1d1cf6e.png (remember that numerical_differential_equations_f74052048d67c9e5264d4eda9679c41ed6ebbe2a.png are integrals).

We can do better; we say that quadrature rule has order numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png if it exactly integrates any polynomial numerical_differential_equations_15a5947693ae6d0d3d895ee4a316de6ff225389b.png. It can be shown that numerical_differential_equations_5da89a289291ee26d84a7ad21fc1709954b97e3c.png always hold, and, for optimal choice of numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png, we have numerical_differential_equations_88d1db2925daa9358a20898ff0661988201de22c.png.

That is, if we want to exactly integrate some polynomial of order numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png, then we only need half as many quadrature points numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png!

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

x = np.arange(10)
y = x + 6 * x ** 2 + 0.4 * x ** 3 - 0.01 * x ** 5


def polynomial_fit(abscissae, y):
    s = abscissae.shape[0]
    zeros = np.zeros((s, s - 1))

    for i in range(s):
        if i == 0:
            zeros[i] = x[i + 1:]
        elif i == s - 1:
            zeros[i] = x[:i]
        else:
            zeros[i] = np.hstack((x[:i], x[i + 1:]))

    denominators = np.prod(abscissae.reshape(-1, 1) - zeros, axis=1) ** (-1)

    def p(x):
        return np.dot(np.prod(x - zeros, axis=1) * denominators, y)

    return p

p = polynomial_fit(x, y)

xs = np.linspace(0, 10, 100)
ys = [p(xs[i]) for i in range(xs.shape[0])]
plt.plot(xs, ys)
plt.suptitle("Polynomial fit")
plt.title("$x + 6x^2 + 0.4x^3 - 0.01x^5$")
_ = plt.scatter(x, y, c="orange")

polynomial_fit.png

One-step collocation methods

Let numerical_differential_equations_0cdbbd4eb99fa1feb3cbe4fd03daf0725d1b7328.png be the collocation polynomial of degree numerical_differential_equations_9493b9776211ef51195c1379ad5f5540a04e566e.png, satisfying

numerical_differential_equations_20516f439413b17a61cdf44eae45a9efdfc46034.png

Then we attempt to approximate the gradient as a polynomial, letting numerical_differential_equations_6625e5cf1a1bd08466cff630f2a877021b330253.png:

numerical_differential_equations_db7f5e794812b07bdcc20e8b6b6953c9ed74be1b.png

where we've let

numerical_differential_equations_5f0b735463a47c052c7537e3f52e7a711037c257.png

Then by the Fundamental Theorem of Calculus, and making the substitution numerical_differential_equations_6625e5cf1a1bd08466cff630f2a877021b330253.png:

numerical_differential_equations_ac4a75037ff3b290c1daf229ec5135ad3dc50281.png

And, using our Lagrange polynomial approximation:

numerical_differential_equations_b3e39c26672f7d3f688f36ac45de0562f7d5b5b5.png

we get

numerical_differential_equations_79ab5b04a96ff8297b89cf859c15d1990fea2aa8.png

Aaand finally, numerical_differential_equations_b98e8e8f2d058554b264030dd455afcb5dcd3a94.png as stated earlier, hence

numerical_differential_equations_c728379e5832b675b6a1decef2e90ec8d959578f.png

for all numerical_differential_equations_5d69a6e9ce33582e94dfa5df7e85940b962a2904.png. And since numerical_differential_equations_5c04d82169398cf1e8be090be69367e76c261baa.png, we get

numerical_differential_equations_bffb1963181b1824c526c54bd962ba70e765d0bc.png

which gives us a system of (non-linear) equations we need to solve simultaneously.

And therefore, our next step numerical_differential_equations_64474c660548e5020014aef2973303a5ba21b6d3.png is given by

numerical_differential_equations_f2d4b06137b90514048c85656023e03e8cad9275.png

where we've let

numerical_differential_equations_7630ab46d810c28e7d9f6bb8721e64df29380fe1.png

Thus,

numerical_differential_equations_a89f5d05dd6b4bae259a17f26816949529ddd3af.png

Let numerical_differential_equations_0cdbbd4eb99fa1feb3cbe4fd03daf0725d1b7328.png be the collocation polynomial of degree numerical_differential_equations_9493b9776211ef51195c1379ad5f5540a04e566e.png, satisfying

numerical_differential_equations_20516f439413b17a61cdf44eae45a9efdfc46034.png

and let numerical_differential_equations_e3c2acd8fda16fe4a737b744feab5cbfe8ecaf48.png be the values of the (as of yet undetermined) interpolating polynomial at the nodes:

numerical_differential_equations_8151e9a082e3be42a2f146840641deaa27d3ee87.png

Then we attempt to approximate the gradient as a polynomial, letting numerical_differential_equations_6625e5cf1a1bd08466cff630f2a877021b330253.png:

numerical_differential_equations_db7f5e794812b07bdcc20e8b6b6953c9ed74be1b.png

The following is something I was not quite getting to grips with, at first.

Now we make the following observation:

numerical_differential_equations_f24941b817d95c7db6723cc6a4d4549824856f33.png

by the Fundamental Theorem of Calculus. Thus, making the substitution numerical_differential_equations_6625e5cf1a1bd08466cff630f2a877021b330253.png, gives

numerical_differential_equations_3ba974a530346124735f352d548f7836abdba55a.png

And, using our Lagrange polynomial approximation:

numerical_differential_equations_b3e39c26672f7d3f688f36ac45de0562f7d5b5b5.png

we get

numerical_differential_equations_79ab5b04a96ff8297b89cf859c15d1990fea2aa8.png

Aaand finally, numerical_differential_equations_b98e8e8f2d058554b264030dd455afcb5dcd3a94.png as stated earlier, hence

numerical_differential_equations_c728379e5832b675b6a1decef2e90ec8d959578f.png

for all numerical_differential_equations_5d69a6e9ce33582e94dfa5df7e85940b962a2904.png.

Let us denote

numerical_differential_equations_b2e2cfb463304f3c597dc773854ed3fa2746e8c4.png

Then, a collocation method is given by

numerical_differential_equations_310ea7c787ac416bdae1e02215e6a81b3bcbb3cf.png

where one first solves the coupled sd-dimensional nonlinear system numerical_differential_equations_e3c2acd8fda16fe4a737b744feab5cbfe8ecaf48.png, and the update numerical_differential_equations_6ddbdcab3e6b9c5be1795e98eef069edfef5847e.png explicitly.

Important: here we're considering numerical_differential_equations_71e6b2f46ab6edf1e466465f10ffdfe2e8204352.png, rather than the standard numerical_differential_equations_684f34a94654106dd5afb679c9688b420cda9800.png, but you can do exactly the same but by plugging in numerical_differential_equations_82853626800c21e3b5d26351355e64bd3e86929c.png in the condition for numerical_differential_equations_51a7fb173967bfa12d6080495a8320a173104531.png instead of what we just did.

Comparing to interpolating polynomial, numerical_differential_equations_a8186ddb0dfdf0e1e9b80877afa28ef5f21d24f3.png becomes the numerical_differential_equations_5b75907baddc1f20479fae69727ff4819121a7ed.png sort of.

Using the collocation methods we obtain a piecewise continuous solution, or rather, a continuous approximation of the solution numerical_differential_equations_a01e4d2c2be3d717c85a6023ebafbec516e66e33.png on each interval numerical_differential_equations_ec7f3c691f3d5ab3765d62a0b073b377de2fd154.png.

In terms of order of accuracy, the optimal choice is attained by using so-called Guass-Legendre collocation methods and placement of nodes at the roots of a shifted Legendre polynomial.

Runge-Kutta Methods

A natural generalization of collocation methods is obtained by:

  • allow coefficients numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png, numerical_differential_equations_fba4537003dea0ac5aea2814cb6a3fac3cc1c7f7.png and numerical_differential_equations_4b87dc8c0f4546682a9531cb72f9bda044b83000.png to take on arbitrary values (not necessarily related to the quadrature rules)
  • no longer assume numerical_differential_equations_b628764d0635802d2b9bfc6478845b6323e92518.png to be distinct

We then get the class of Runge-Kutta methods!

We introduce the stage values numerical_differential_equations_7a9c1f2eb8260cbc3a6cb21f1c3eb4f357e9dd01.png, ending up with what's called the Runge-Kutta methods:

numerical_differential_equations_baf9256f51515b6d6fe6531cacf402c7bc2f2ed1.png

where we can view numerical_differential_equations_7a9c1f2eb8260cbc3a6cb21f1c3eb4f357e9dd01.png as the intermediate values of the solution at numerical_differential_equations_a3a7f43f807b9e381fc50e0fab140c0df0a03e17.png at time numerical_differential_equations_bc997069e7845a097fc7c1a9d72354ec84e83342.png.

For which we use the following terminology:

  • numerical_differential_equations_9493b9776211ef51195c1379ad5f5540a04e566e.png the number of stages of the Runge-Kutta method
  • numerical_differential_equations_f74052048d67c9e5264d4eda9679c41ed6ebbe2a.png are the weights
  • numerical_differential_equations_4b87dc8c0f4546682a9531cb72f9bda044b83000.png are the internal coefficients

Euler's method and trapezoidal rule are Runge-Kutta methods, where

numerical_differential_equations_b768254f28cf14c65bdce6a36b15a6aff78ab1e0.png

gives us the Euler's method.

If the matrix numerical_differential_equations_2d14b9d49b12c938696bf9b3754a651fec2ea2be.png is strictly lower-triangular, then the method is explicit.

Otherwise, the method is implicit, as the we might have "circular" dependency between two stages numerical_differential_equations_7a9c1f2eb8260cbc3a6cb21f1c3eb4f357e9dd01.png and numerical_differential_equations_7f258efac274f9dc04bd3e033fb7d9b7302eb446.png, numerical_differential_equations_a41c04a3f3dafda05d9e0c209b6ed8613c844475.png.

Examples of Runge-Kutta methods

Runge-Kutta method with 4 stages - "THE Runge-Kutta method"

numerical_differential_equations_7a3b2243b7a124b9f9d00c9a4a5ecc29e0a688a2.png

Which often are represented schematically in a Butcher table.

Accuracy of Runge-Kutta methods

Notation
Stuff
  • In general, we cannot have an exact expression for the local error numerical_differential_equations_7f829f6228cad80b5f32d9eafd616686b6bb3103.png
  • Can approximate this by computing Taylor expansion wrt. numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png

Consider the case of continuous mechanics, the derivates can be related directly to the solution itself by using the differential equation:

numerical_differential_equations_b2b0128394bf8b8fef4d104cf240e3493bcc9df9.png

Then, we can write the flow map as

numerical_differential_equations_6b458daddd071a075bb8dee08240d9d3fa4e5745.png

Order conditions for Runge-Kutta methods

With

numerical_differential_equations_d96107a416ab8dacedc0bddb42749e5ea2ca295b.png

we have

numerical_differential_equations_9a8b747cb00dae81b7ad93efb57b3a8ec2fdb512.png

For the method to have order of numerical_differential_equations_45b45d58a4ae4130d828ee8475091cc1a733b764.png, we must have

numerical_differential_equations_4231cbadf3f7adf9bd20e9201e0801cec7c8082b.png

since we will then have

numerical_differential_equations_64a3843ff8783720f929d5dd55e17eec59726184.png

and thus, the method will be of order numerical_differential_equations_45b45d58a4ae4130d828ee8475091cc1a733b764.png. This is the first what is termed as order conditions.

Then, to get the second order condition, we simply compute the second order derivative of numerical_differential_equations_0e3a01b1869d914b6d3ee46d8f6de341debf1176.png:

numerical_differential_equations_3b3be47d9d627e84fe89caaadcfc3fb060dc8951.png

Substituting in the expressions for numerical_differential_equations_e3c2acd8fda16fe4a737b744feab5cbfe8ecaf48.png and numerical_differential_equations_7a9c1f2eb8260cbc3a6cb21f1c3eb4f357e9dd01.png, we obtain

Variable stepsize

Notation

  • numerical_differential_equations_13bc26d1905131cfb066cc647837476d39433727.png and numerical_differential_equations_d2fe3246943fcd682d44656a1dd8cc165e5f984e.png are two different methods for numerical approximations to the flow map
  • numerical_differential_equations_13bc26d1905131cfb066cc647837476d39433727.png is order numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png (local error numerical_differential_equations_2591f42cb08c519be0667aa5d14c8314a391f384.png)
  • numerical_differential_equations_d2fe3246943fcd682d44656a1dd8cc165e5f984e.png is order numerical_differential_equations_0d2cbc4c152441edc328e6dc119bac015f4ab42a.png (local error numerical_differential_equations_b87c11bafdd32f7803ab6e11652e3d5fb72f1b12.png)
  • numerical_differential_equations_9275764a36cfa783f322d0258a386d4a3b0867f4.png

Stuff

Difference in an estimate of the local error by two different methods

numerical_differential_equations_5fc6e2f8f1dc8829f9cb3da98b0d330fa3976f92.png

And we assume the local error satisfies

numerical_differential_equations_2b490947907597b8828a532a1e74c7f82fdb9fe6.png

i.e. that numerical_differential_equations_2591f42cb08c519be0667aa5d14c8314a391f384.png is the dominating factor.

Error control

Assumptions:

  • numerical_differential_equations_13bc26d1905131cfb066cc647837476d39433727.png is a convergent one-step method
  • Desire local error at each step less than a given numerical_differential_equations_de2062a31663505a89681779aa24b7a555498684.png
  • numerical_differential_equations_21554d84f8abcf2f7d44d2526ce26cfe10314555.png estimates the local error at a given point numerical_differential_equations_b1d58f875e7a17df47fe1b0499e55f50d9a0eb34.png when we take a step size of numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png

Idea:

  • Estimate the local error with numerical_differential_equations_9dfb630910b0b3d2a61cb5f1673b9adaa08fe150.png
  • If numerical_differential_equations_ecf9b693ab704a73fe1dbe058928dd94bda97f6e.png then decrease numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png and redo step
  • If numerical_differential_equations_a9a31d242ae61ebba22011d7911cce58157b9299.png then keep step; may want to increase numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png

Question: how much do we decrease / increase numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png?

Suggestion: raise/lower timestep by a factor numerical_differential_equations_195548b99a9939528c7a54b1136fd2f39dd22dee.png

  • Can be inefficient in the case where we need to do this guessing multiple times, can instead do:

    numerical_differential_equations_cbf0867ae5500b70010e3c5907c0aaa9da419f63.png

    for some numerical_differential_equations_3680099aad3d20fd7f3b82ceb4080bfd2c3acfec.png, e.g. numerical_differential_equations_a21550c8e2b34a9de19e7a4f8a13c8f749eddfd5.png, i.e. numerical_differential_equations_339c704fc9bf0c2f6780a0a0e68b51e128bb7f74.png is a decaying factor for the "guess" of the stepsize.

Stability of Runge-Kutta Methods

Notation

  • numerical_differential_equations_e470da60c1191889804b5aca80404b884e0cd6c4.png denotes the set of fixed points of the flow map numerical_differential_equations_9d0a50f0aae06c250b35249a4a5b9f4e038b9f6c.png
  • numerical_differential_equations_4d4e7bcbba999b50bcebd85350518b71112a3643.png denotes the set of fixed points for the numerical estimation of the flow map
  • numerical_differential_equations_cc816f175ac6403febb6e1766efbd4bab5309e0f.png denotes the spectral radius of a matrix numerical_differential_equations_d80b633218f04f7a25f6a9dd7e4617a38e88f263.png
  • numerical_differential_equations_76f6ebb93e725a151a75229319b2b57af193e656.png is the ratio between two polynomials numerical_differential_equations_c64796e31b0c60d6de085f3d8856c916c47eb91f.png and numerical_differential_equations_7c7dd9689e54f37527c3218030f62b7fdcc01249.png called the stability function

Stuff

An equilibrium point numerical_differential_equations_688e6e5747f484cd8a4e57dfea18e23bd317e255.png of the scalar differential equation

numerical_differential_equations_ebe36757f2fd129b1e392db6f5b954a709eb3314.png

is a point for which

numerical_differential_equations_bdf1c699d41710879d9b33796a6aa35ee3b735f2.png

Fixed points which are NOT equilibrium points are called extraneous fixed points.

These are points such that

numerical_differential_equations_ffa9d5974f16c0ddf2b398e2f86d4e49af569443.png

Stability of Equilibrium points

Suppose that numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png in

numerical_differential_equations_9c27fcf6f0a9ee9697c9ec3243ef652cbac307f8.png

is numerical_differential_equations_92bd40c1917561cb3f61e41c484b1f835756c6db.png, i.e. numerical_differential_equations_feac5f1bc304042c6a7f37a53462e1840455ec0c.png, and has equilibrium point numerical_differential_equations_688e6e5747f484cd8a4e57dfea18e23bd317e255.png.

If the eigenvalues of

numerical_differential_equations_1bf03b030287825f251441c35b2a12897df3a408.png

all lie strictly in the left complex half-plane, then the equilibrium point numerical_differential_equations_688e6e5747f484cd8a4e57dfea18e23bd317e255.png is asymptotically stable.

If numerical_differential_equations_31783bf9ca68d36cf431288f11cfe6cf60a14ef3.png has any eigenvalue in the right complex half-plane, then numerical_differential_equations_688e6e5747f484cd8a4e57dfea18e23bd317e255.png is an unstable point.

Stability of Fixed points of maps

We say a numerical method is A-stable if the stability region includes the entire left half-plane, i.e. numerical_differential_equations_a4a87316f3c46331a99cc293d8168b155a9e532c.png is a contained in the stability region.

A A-stable numerical method has the property that it's stable at the origin.

I was wondering why numerical_differential_equations_311bb02290e1c19d78b19f9e962065c5c6c0be7e.png? I believe it's because we're looking at solutions which are exponential, i.e. numerical_differential_equations_4fbd4c8a9e26558df5c255c0f416f1750392580a.png, hence if numerical_differential_equations_a1a667e2b5d55e02d95239882100e02bc5e002f1.png, then the exponential will have modulus smaller than 1:

numerical_differential_equations_aa702a45746abb496f34dba1290f09a00012876f.png

if numerical_differential_equations_311bb02290e1c19d78b19f9e962065c5c6c0be7e.png, for numerical_differential_equations_8a27571bd1409a3a4bfab6e0de2f41ad3d12d1de.png.

We say a numerical method is L-stable if it has the properties:

  1. it's A-stable
  2. numerical_differential_equations_45b6604053ccb7f1f15c0951da46af716b6c0e7a.png as numerical_differential_equations_56bbb8ec3f5d3bcf542ae21e549434ec2372f8b5.png

Stability of Numerical Methods

Stability functions
  • Scalar case

    Consider the scalar first-order diff. eqn

    numerical_differential_equations_4837c8ea367061bc521768b4dd6bfd23197a358f.png

    When applied to this diff. eqn., many methods, including all explicit Runge-Kutta methods, have the form

    numerical_differential_equations_7365ea4fe250f286dfe63bda1a36ac19ee62a21d.png

    for some polynomial numerical_differential_equations_f8bec95b572243552399d6897a7cf76ffb3f9118.png.

    More generally, all RK methods have the form

    numerical_differential_equations_02a8bc339490ddd74811f97d108879c79e740ff3.png

    where numerical_differential_equations_77a66189e3eddfb016541a3663fe397662008bd9.png is rational polynomial, i.e.

    numerical_differential_equations_cc97821a5e08053b4e012c4e2ce9b2e06e546b85.png

    Then we call numerical_differential_equations_6dd4ec21f92a21e77e43374b4c7fd1cc8c2234c6.png the stability function of the RK method.

    numerical_differential_equations_f85ef8288dfbf56b0c47ae36d5609dcab4f355be.png

    numerical_differential_equations_6fa656a2ad2885c399345f9a8ec2949a3867104c.png, then we can write this sum in a matrix form:

    numerical_differential_equations_1bb0b89313549420ad93c7468968cf75ad6d41a7.png

    which has the solution

    numerical_differential_equations_c8036cb2f36d4c61a82459113a8fc95d832e6b64.png

    and therefore, finally, we can write

    numerical_differential_equations_c7e0475ef9bfebda03dcec46c70b8d9f96c78340.png

    numerical_differential_equations_927c6d81ea8dca10f2b1642789e7fd452e4c2aaa.png

    which gives us the stability function numerical_differential_equations_6dd4ec21f92a21e77e43374b4c7fd1cc8c2234c6.png for RK:

    numerical_differential_equations_5972ead1bd112985beabb404e335b8e16ea1a661.png

Stability of Numerical Methods: Linear case

Linear system of ODEs

numerical_differential_equations_d2f25dd8d326b07a0ab1ff8a42de40419415b52d.png

where numerical_differential_equations_e0e8c4da2f4377b71762ac6ef5cdd8187c4b0be3.png is diagonalisable. Can write solution as

numerical_differential_equations_e127d42b21657c8c2f22357ef4d7f57c93785867.png

numerical_differential_equations_d2f25dd8d326b07a0ab1ff8a42de40419415b52d.png

where numerical_differential_equations_aa45b78ca4fcfb1128b81983bc69cebaddb6486d.png is a diagonalisable matrix.

Further, let an RK method be given with stability function numerical_differential_equations_6dd4ec21f92a21e77e43374b4c7fd1cc8c2234c6.png. The origin is stable for numerical method is stable if and only if

numerical_differential_equations_4ccc4469e873759d77db997261f650173ed2773a.png

Linear Multistep Methods

Notation

  • numerical_differential_equations_53d38fed0a51bcb2f1ae8f0904e7ff3f946e2ffd.png
  • numerical_differential_equations_98fcfcfeb7955153183967280e6969b658c2f23b.png where numerical_differential_equations_6b97edc8dd09f1ede9c25803d5eeb6e0c84beac7.png

Definitions

A linear k-step method is defined as

numerical_differential_equations_317e1e61f76c1a8a803276c689e31d19a8410c91.png

where numerical_differential_equations_c9aff27c6b7e229eaf939cbbb05a6065cbdeef1a.png and either numerical_differential_equations_5b620bdf48cd9f102d7094545c6077a641236238.png or numerical_differential_equations_a416a2eaaef624443e062a92f37b129f92cff15c.png.

The coefficients are usually normalized such that either

numerical_differential_equations_b31c1f31d33cba88b85392cabcfc24f6b9cbd7f8.png

In implementation is's assumed that the values numerical_differential_equations_93d6594e498bf552fb4039c68dea0fe75dc20475.png are already computed, so that numerical_differential_equations_94f18eeba7de458cbdf3c2ac548facd7cec2ec33.png is the only unknown in the formula.

If numerical_differential_equations_73810f242dd319bfb19919bac8ef54e74ec63003.png is non-zero the method is implicit, otherwise it's explicit.

Examples

numerical_differential_equations_e7a86163f39fa21d4a2ed66946369cdeb900ef42.png method - generalization of one-step methods

The numerical_differential_equations_e7a86163f39fa21d4a2ed66946369cdeb900ef42.png multistep method is of the form

numerical_differential_equations_7e4f92fd02d73079b068b721e267236d19c537f7.png

which is a linear k-step method with coefficients numerical_differential_equations_ac696b897ace457137e349ec1f59edc78a69c79f.png, numerical_differential_equations_0c4321b2a4bccc5b824b1d42c91481a25fdd1fcf.png, numerical_differential_equations_381392986f015d5f651df6d9807ac5424a057da8.png and numerical_differential_equations_94cddfc881c57ad2f0fc08c1bc771f6c15582dc4.png.

Examples:

Leapfrog

Leapfrog is an explicit two-step method (numerical_differential_equations_d1ce0769b32ee0f544eb7d7cfd83ff873317145f.png) given by

numerical_differential_equations_439ec9feb7a4c963fa91dd96b1ff4f409fff4a35.png

i.e. numerical_differential_equations_ac696b897ace457137e349ec1f59edc78a69c79f.png, numerical_differential_equations_48d9dbf83cdfae48281c4e5607fde8ff4cea638b.png, numerical_differential_equations_6b8ece6c53fafb3ef5c44631f24c02d341f722fd.png and numerical_differential_equations_ce04cdc7f8911ec05f89adfd9ea2984c71191e2b.png.

Adams methods

The class of Adams methods have

  • numerical_differential_equations_cce8af81a0435125e0b579d5b04aaf10b6793787.png
  • numerical_differential_equations_ec3e03de70af40992dab0ecd9f8743167ea56832.png
  • numerical_differential_equations_285350b2e63fcf67cdb9bc0454161c054c0ab2a3.png for numerical_differential_equations_c46b4d193e5e4fbb81a19f5f66241e14361add27.png

If we also have the additional numerical_differential_equations_53c785d82703c2ab95ed93c933606b7eecde4b7c.png, we have explicit Adams methods, known as Adams-Bashforth methods, e.g.

numerical_differential_equations_8fa4afd9a51bb8cd7f6635103316ac53899fdc72.png

where numerical_differential_equations_0fa1476e4e961926940fcb2117f14487599c06eb.png.

Adams-Moulton methods are implicit, i.e. numerical_differential_equations_9689f13c5a56e41e53cd8fb181af82fb86ee4f99.png.

Backward differentiation formulae (BDF)

The Backward differentiation formulae (BDF) are a class of linear multistep methods satisfying

numerical_differential_equations_897e78072ec424328f1e98d1b78c26aef48ab483.png

and generalizing backward Euler.

E.g. the two-step method (BDF-2) is

numerical_differential_equations_a2aa7453c09dd00a855e803ecb567a061a167531.png

Order of accuracy

Associated with the linear multistep method are the polynomials

numerical_differential_equations_02260f97b2a8f073be16841f8754dfc14c950e8e.png

Let numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png be the exact solution at times numerical_differential_equations_9f9f8efc59745fbcd492ccb1168629b643ca227a.png for numerical_differential_equations_e0cff938190a47873c4117407d7157cd594edb7c.png.

Then, subsituting numerical_differential_equations_03f5c1abbdf9cf9665afac4d6d9fe988f91b7b20.png into the linear k-step method expression

numerical_differential_equations_73d7b37fe22fb7b3880406cf46be944958c8f617.png

(which is actually the residual accumulated in the numerical_differential_equations_0600a9f2278fdce30f15b32e3430569ec8bd03e7.png th step, but for notational convenience we will denote it numerical_differential_equations_e64178a349205cceadd9f36788e2cd3ac0ff3f8a.png).

A linear multistep method has order of consistency numerical_differential_equations_fefe9e556d399665a26a37824ec578cbffb0cabe.png if (with numerical_differential_equations_e64178a349205cceadd9f36788e2cd3ac0ff3f8a.png being the residual)

numerical_differential_equations_373edf49bd5852db85fcd585c533b91f6ca1d6c6.png

for all sufficiently smooth numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png, or (it can be shown) equivalently have the following properties:

  • Coeffiicents numerical_differential_equations_799fca826aaa60e554f69bd8f737276dc3af499f.png and numerical_differential_equations_be5b8329c4c3afbfdfa3ea12837c89609bfdbc5e.png satisfy

    numerical_differential_equations_5e6d3be4561e3c4e5dffa03522f83f4c1ee9ab66.png

  • The polynomials numerical_differential_equations_3b2fbc1e8fe1103b19ba84179700d0bb4913c850.png and numerical_differential_equations_82bc9b57fa58e0079fe68566314cb97df9c15dc9.png satisfy

    numerical_differential_equations_728f38e6dcd087f7320ffda654b519ae43bc76fb.png

  • The polynomials numerical_differential_equations_3b2fbc1e8fe1103b19ba84179700d0bb4913c850.png and numerical_differential_equations_82bc9b57fa58e0079fe68566314cb97df9c15dc9.png satisfy

    numerical_differential_equations_39d0f993e3e2e1580ba66447a322e0358bedfaa3.png

Unlike the single-step methods, e.g. Euler method, consistency is insufficient to ensure convergence!

Convergence

A linear multistep method is said to satisfy the root condition, if all roots numerical_differential_equations_5f2ba5d98a1201e7d39a0c689133d800b4680423.png of

numerical_differential_equations_21ffc53c5344ea64e0883eccc69e574ac8886d9c.png

lie in the unit disc, i.e. numerical_differential_equations_94ffbfe8d5d26ae248051bb7e81d7323d6eaec8e.png, and any root on the boundary, i.e. numerical_differential_equations_01dd09c8dcd0202c5deaaf50abcf8ca9a85963ca.png, has algebraic multiplicity one.

Otherwise, the modulus of the solution grows in time.

Suppose a linear multistep method is quipped with a starting procedure satisfying

numerical_differential_equations_2aaa3fc1de5a80e12bb673fe84f26c89d4a7b50e.png

Then the method converges to the exact solution of

numerical_differential_equations_aea280f77c00e3df5007df14546c2b72661bd6f9.png

on a fixed interval as numerical_differential_equations_d580fd7dcc4f0f0ed2e3ab37a90b197e7e19d698.png if and only if it has order of accuracy numerical_differential_equations_813342791ebbfd19426b72f6a4b92a237d8049f2.png and satisfies the root condition.

Weak Stability and the Strong Root Condition

A convergent multistep method has to be:

  1. Consistent: so numerical_differential_equations_024d03a02a9d88c55286d6e5c947fd215b9412d9.png and numerical_differential_equations_3f7f0da4f26304827d3901786629ad4f232df9fa.png where numerical_differential_equations_9ab766fe2c081b82a304866d50f951fadbdb5704.png and numerical_differential_equations_f2eebbd62e41be9d5a2457b4bd291dd096896884.png are the characteristic polynomials
  2. Initialized with a convergent starting procedure: i.e. a method to compute numerical_differential_equations_77d159763a83ea66f5b921aa213d3a4bd8da1a0a.png for the first multistep iteration, to get started, satisfying

    numerical_differential_equations_71a28d91e4eb914be7282c91e5fd30b750aa19b7.png

  3. Zero-stable: i.e. it satisfies the root condition; the roots of numerical_differential_equations_3eb17318fa26ae708e30a41f656a6f50985a8fe8.png are in the unit disc and simple if on the boundary.

If they're not simple, we get

numerical_differential_equations_caf815b544b5b6bfb747f3411a998a272eb3b178.png

as we iterate, which is a problem!

We refer to a multistep method that is zero-stable but which has some extraneous roots on the boundary of the unit disk as a weakly stable multistep method.

Same as root condition, but all roots which are not numerical_differential_equations_1946721e0e622f1cf94dc1d6d823cb61287f1567.png are within the unit disk.

Asymptotic Stability

The stability region numerical_differential_equations_3b4672f4daa50549084f26685be4992faceeb7ec.png of a linear multistep method is the set of all points numerical_differential_equations_9c15196dd07b1add486b8b54592e74bfe946ed95.png such that all roots numerical_differential_equations_5f2ba5d98a1201e7d39a0c689133d800b4680423.png of the polynomial equation

numerical_differential_equations_a6b45b702c30d6897ade986740d72784523b913d.png

lie on the unit disc numerical_differential_equations_94ffbfe8d5d26ae248051bb7e81d7323d6eaec8e.png, and those roots with modulues one are simple.

On the boundary of the stability region numerical_differential_equations_3b4672f4daa50549084f26685be4992faceeb7ec.png, precisely one root has modulus one, say numerical_differential_equations_f4d31d2b946535648ee75858f56002cf8af70278.png. Therefore an explicit representation for the boundary of numerical_differential_equations_3b4672f4daa50549084f26685be4992faceeb7ec.png is:

numerical_differential_equations_f622455711f7d9eb18036b73bd8dfa853213d95e.png

This comes from applying the multistep method to the diff. eqn. (Dahlquist test equation)

numerical_differential_equations_4837c8ea367061bc521768b4dd6bfd23197a358f.png

Then

numerical_differential_equations_4596488a61d704008200fd287c65aa13d453b19d.png

Letting numerical_differential_equations_98fcfcfeb7955153183967280e6969b658c2f23b.png, we get

numerical_differential_equations_53df994cfe44b5ea1fbf74b6a4fb7e15c41ebd82.png

For any numerical_differential_equations_9c15196dd07b1add486b8b54592e74bfe946ed95.png, this is a recurrence relation / linear difference with characteristic polynomial

numerical_differential_equations_b128ea381f41b2a29c95f7d1ab66216ef723aa35.png

The stability region is numerical_differential_equations_8a27571bd1409a3a4bfab6e0de2f41ad3d12d1de.png s.t. all roots of numerical_differential_equations_b1249e0adf552ec5309f3e6c6eec6badb45590f1.png lie in the unit disc and those on the boundary are simply, which is just what we stated above.

Geometric Integration

Focused on constructing methods for specific problems to account for properties of the system, e.g. symmetry.

Notation

First integrals

For some special functions numerical_differential_equations_7c06d00663ab38a8fa27deb6d2bc3ee00181ce60.png we may find that numerical_differential_equations_6c03018fd4253198b850bc6339e99699b567a798.png is constant along every solution to an /autonomous differential equation

numerical_differential_equations_2613b854b7739e9ec81125d0a3eaa7a5bcd19791.png

Such a function numerical_differential_equations_6c03018fd4253198b850bc6339e99699b567a798.png is called a first integral.

This closely related to the conservation we in Hamiltonian physics.

Let

numerical_differential_equations_ebe36757f2fd129b1e392db6f5b954a709eb3314.png

have a quadratic first integral. A Runge-Kutta method preserves this quadratic first integral if

numerical_differential_equations_088136fb4e867b63827a197e72e9063af2888606.png

Splitting methods

  • Useful for constructing methods with specific properties

Suppose we wish to solve the differential equation

numerical_differential_equations_257fefbb1910024ee3cc06c4d85df3cde1d9ad4c.png

where each of the two differential equations

numerical_differential_equations_e15e023220b32429c0acf10904af9459797539f0.png

happen to be completely integrable (i.e. analytically solvable).

Then, any point in the phase space, the vector field can be written as two compontents numerical_differential_equations_866ead89ec63f64bf1f7216d2aac4c2c37c2c921.png and numerical_differential_equations_8523abf12e4b850e630a76804e0f6241b0b93067.png. We first step forward along one tangent, and then step forward in the second, each time for a small timestep numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png.

That is, we start from some point numerical_differential_equations_2cf02d5943c9689ad064b553c5a546d1a23c9e9a.png and solve first the IVP

numerical_differential_equations_0e593189f24d959acbc1c4355ce2f022171d5862.png

for time numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png, which takes us to the point numerical_differential_equations_ce3154aba7f159166f98ebb7993c3ff756985fd6.png where numerical_differential_equations_8695baa5ba4d5b19604a9593a10a84def7b6564b.png represents the exact solution of the differential equations on vector field numerical_differential_equations_866ead89ec63f64bf1f7216d2aac4c2c37c2c921.png.

Next, we start from numerical_differential_equations_544e971b7c542f60831d1bfc5142af2754ad68f3.png and solve the IVP

numerical_differential_equations_28b692af98ea56d7878099fdb021c2c9c85162b5.png

for time numerical_differential_equations_3807b8495f01004f1d88370e2f4fad2f5032db06.png, taking us to the point numerical_differential_equations_8e2a015972142b02b665a30de0dc32ec2e0d5238.png.

Denoting the flow maps of the vector fields numerical_differential_equations_866ead89ec63f64bf1f7216d2aac4c2c37c2c921.png and numerical_differential_equations_8523abf12e4b850e630a76804e0f6241b0b93067.png as numerical_differential_equations_bdf5541528333026a0ac5d62e0678ddb5f0d3779.png and numerical_differential_equations_e87b42649b3621c1aeaa8db71706b0024ff12c74.png we have just computed

numerical_differential_equations_70efffce6a643581e89effbe15f3cfe4399849f3.png

Such a method is called a splitting method.

The splitting method numerical_differential_equations_73e21d43312a3ee844db87c50227b4e4a539387b.png where numerical_differential_equations_e72fc2adb55f3e441cd9901212964be4d6ca1981.png, has local error

numerical_differential_equations_42af48862b6a482feaad5ff2bc6f311c43719564.png

where numerical_differential_equations_18f478091698d2f132658854b5dee04a51e6ec51.png are the commutator of the vector fields numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png and numerical_differential_equations_bed561b338628a088ae69301d4bba9fd83a70bd2.png, which is another vector field

numerical_differential_equations_553e1af8fcd25925f618616eda380680dec2389b.png

More generally

The splitting method numerical_differential_equations_73e21d43312a3ee844db87c50227b4e4a539387b.png where numerical_differential_equations_e72fc2adb55f3e441cd9901212964be4d6ca1981.png, has local error

numerical_differential_equations_841a80071222a4c1fc177b5cae3390b242439eec.png

where numerical_differential_equations_18f478091698d2f132658854b5dee04a51e6ec51.png are the commutator of the vector fields numerical_differential_equations_cdd1cc131da6040eca078917132a377727053c44.png and numerical_differential_equations_bed561b338628a088ae69301d4bba9fd83a70bd2.png, which is another vector field

numerical_differential_equations_553e1af8fcd25925f618616eda380680dec2389b.png

If the fields commute, i.e. numerical_differential_equations_28e9815026d040c523d4e34faff7b5719606e0f2.png, then the splitting is exact.