Samer Adeeb

Introduction to Numerical Analysis: Ordinary Differential Equations

Ordinary Differential Equations (ODEs) are equations that relate the derivatives of one or more smooth functions with an independent variable x\in\mathbb{R}. The term “ordinary” indicates that there is only one independent variable appearing in the equations.

Examples of Ordinary Differential Equations

ODEs appear naturally in almost all engineering applications. Here are some examples:

Newton’s Second Law of Motion

One of the most ubiquitously used ordinary differential equations is Newton’s second law of motion, which relates the second derivative of the position of a particle (i.e., the acceleration) to the applied force on the particle. The relationship has the form:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=F(x,t) \]

where m is the mass of the particle. The time t is the independent variable, while the dependent variable is the position x(t). Newton’s equation adopts different forms depending on the application. For example, if the applied force is the combined action of gravity and a spring with a constant k>0 (Figure 1), then the equation has the form:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=mg-kx \]

The explicit solution to this equation has the following form:

    \[ x=\frac{mg}{k}+C_1\cos\left(\frac{\sqrt{k}t}{\sqrt{m}}\right)+C_2\sin\left(\frac{\sqrt{k}t}{\sqrt{m}}\right) \]

where C_1 and C_2 are constants that depend on the initial conditions of the problem. If the initial displacement and velocity are equal to zero, i.e., x(0)=0 and x'(0)=0, then, the explicit solution has the form:

    \[ x=\frac{mg-mg \cos\left(\frac{\sqrt{k}t}{\sqrt{m}}\right)}{k} \]

The following code utilizes the built-in DSolve function in Mathematica to solve the above simple equations. Also, the Manipulate built-in function is used to visualize the effect of varying m and k on the resulting vibration of the system.

View Mathematica Code
Clear[x, g, k, m]
(*General Solution*)
a = DSolve[m*x''[t] == m*g - k*x[t], x, t];
Print["General Solution="]
x = x[t] /. a[[1]]

(*Solution with initial velocity and displacement =0*)
Clear[m, g, x, k]
a = DSolve[{m*g - k*x[t] == m*x''[t], x[0] == 0, x'[0] == 0}, x, t];
Print["Solution with initial velocity and displacement =0"]
x = x[t] /. a[[1]]

Manipulate[g = 1;
 Graphics3D[{Cuboid[{-1, -1, -5}, {0, 0, (g*m - g*m*Cos[Sqrt[k/m]*t])/k}]}, Axes -> True, PlotRange -> {{-1, 0}, {-1, 0}, {-5, 2*5/1}},  AxesOrigin -> {0, 0, 0}, ViewVertical -> {0, 0, -1}], {t, 0, 4*Pi*Sqrt[5]}, {m, 1, 5}, {k, 1, 5}]


Figure 1. Mass - Linear Elastic Spring System

Figure 1. Mass – Linear Elastic Spring System

If the external forces include a damping force (a force that is a function of the velocity and that is always opposite to the direction of motion), then, the equation has the form:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=mg-kx-\beta \frac{\mathrm{d}x}{\mathrm{d}t} \]

where \beta>0 is the damping coefficient.
The explicit solution to this equation has the following form:

    \[ x=\frac{mg}{k}+C_1e^{\frac{t}{2}\left(\frac{-\beta-\sqrt{\beta^2-4km}}{m}\right)}+C_2e^{\frac{t}{2}\left(\frac{-\beta+\sqrt{\beta^2-4km}}{m}\right)} \]

where C_1 and C_2 are constants that depend on the initial conditions of the problem. If the initial displacement and velocity are equal to zero, i.e., x(0)=0 and x'(0)=0, then, the explicit solution has the form:

    \[ x=-\frac{m e^{-\frac{t \left(\sqrt{\beta^2-4 k m}+\beta\right)}{2 m}} \left(\beta \left(e^{\frac{t \sqrt{\beta^2-4 k m}}{m}}-1\right)+\sqrt{\beta^2-4 k m} \left(e^{\frac{t \sqrt{\beta^2-4 k m}}{m}}-2 e^{\frac{t \left(\sqrt{\beta^2-4 k m}+\beta\right)}{2 m}}+1\right)\right)}{2 k \sqrt{\beta^2-4 k m}} \]

The following code utilizes the built-in DSolve function in Mathematica to solve the above simple equations. Also, the Manipulate built-in function is used to visualize the effect of varying m, k, and \beta on the resulting vibration of the system.

View Mathematica Code
Clear[x, g, k, m]
(*General Solution*)
a = DSolve[m*x''[t] == m*g - k*x[t] - beta*x'[t], x, t];
Print["General Solution (damped system)="]
x = x[t] /. a[[1]]

(*Solution*)
Clear[x, g, k, m, beta]
a = DSolve[{m*x''[t] == m*g - k*x[t] - beta*x'[t], x[0] == 0, x'[0] == 0}, x, t];
Print["Solution with initial velocity and displacement =0"]
x = x[t] /. a[[1]];
x = Simplify[x]
Manipulate[g = 1;Graphics3D[{Cuboid[{-1, -1, -5}, {0, 0, -((E^(-(((beta + Sqrt[beta^2 - 4 k m]) t)/(2 m)))m (beta (-1 + E^((Sqrt[beta^2 - 4 k m] t)/m)) + (1 + E^((Sqrt[beta^2 - 4 k m] t)/m)-2 E^(((beta + Sqrt[beta^2 - 4 k m]) t)/(2 m))) Sqrt[beta^2 - 4 k m]))/(2 k Sqrt[beta^2 - 4 k m]))}]},  Axes -> True, PlotRange -> {{-1, 0}, {-1, 0}, {-5, 10}},  AxesOrigin -> {0, 0, 0}, ViewVertical -> {0, 0, -1}], {t, 0,  2*4*Pi*Sqrt[5]}, {m, 1, 5}, {beta, 0, 5}, {k, 1, 5}]

Beam Deflection

The deflection of an Euler-Bernoulli beam follows the following equation:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}x^4}=q \]

where E is the beam’s Young’s modulus, I is the cross sectional moment of inertia, y is the beam’s deflection (the dependent variable), q is the applied distributed load on the beam, and x is the position along the beam’s length (the independent variable). For a fixed-ends beam; if the deflection and rotation at both ends of a beam of length L is equal to zero, i.e. y(0)=y(L)=y'(L)=y'(0)=0, and if q is constant, then, the deflection of the beam has the form:

    \[ y=\frac{qx^4-2Lqx^3+qL^2x^2}{24EI} \]

Growth and Decay

If x is the number of organisms at any time t, and if the growth of the number depends on the number that is already there, then, the ODE that would describe this uninhibited growth is given by:

    \[  \frac{\mathrm{d}x}{\mathrm{d}t}=kx \]

Where k>0 is a positive constant. The time t is the independent variable, while the dependent variable is the number of organisms x(t). The general solution of this equation has the form:

    \[ x=C_1e^{kt} \]

where C_1 is a constant depending on the boundary conditions. If the initial condition is such that x(0)=A then, the solution has the form:

    \[ x=Ae^{kt} \]

Similarly, if rate of decrease (decay) of a quantity is a function of the amount available, then, the uninhibited decay can be described by the equation:

    \[  \frac{\mathrm{d}y}{\mathrm{d}t}=-hy \]

where h>0 is a positive constant. The time t is the independent variable, while the dependent variable is the quantity y(t). Radioactive decay follows this ODE. The general solution of this equation has the form:

    \[ y=C_1e^{-ht} \]

where C_1 is a constant depending on the boundary conditions. If the initial condition is such that y(0)=A then, the solution has the form:

    \[ y=Ae^{-ht} \]

The following tool plots the quantities x(t) or y(t) when A=1 for different values of the exponent:

View Mathematica Code
Clear[x]
a = DSolve[x'[t] == k*x[t], x, t];
x = x[t] /. a[[1]]

Clear[x]
a = DSolve[{x'[t] == k*x[t], x[0] == A}, x, t]
x = x[t] /. a[[1]]

Manipulate[A = 1; r = Grid[{{"k=", k}}, Spacings -> 0]; Grid[{{Plot[A*E^(k*t), {t, 0, 20}, AxesLabel -> {"t", "x(t) or y(t)"}, ImageSize -> Medium]}, {r}}], {k, -0.5, 0.5}]

Electrical Capacitors and Inductors

Simple electric circuits are composed of three components: resistances, capacitors, and inductors. The relationship between the voltage drop V_R across a resistance, and the current I_R passing through it is given by:

    \[ V_R=R I_R \]

where R is the constant of proportionality, namely the value of the electric resistance. A capacitor on the other hand is an element that has the following relationship between the current passing through it I_C and the rate of change of the voltage drop across it V_C:

    \[ I_C=C \frac{\mathrm{d}V_C}{\mathrm{d}t} \]

where C is the constant of proportionality, namely the capacitance. An inductor is an element that has the following relationship between the voltage drop across it V_I and the rate of change of the current that passes through it I_I:

    \[ V_I=L \frac{\mathrm{d}I_I}{\mathrm{d}t} \]

Where L is the constant of proportionality, namely the inductance. When the elements are combined, ODEs are generated with t being the independent variable, and V and I being the dependent variables.

Classification of ODEs

Order

The order of an ODE is the order of the highest derivative appearing in the equation. For example, the following equation (Newton’s equation) is a second-order ODE:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=mg-kx \]

while the beam equation is a fourth-order ODE:

    \[ EI\frac{\mathrm{d}^4y}{\mathrm{d}x^4}=q \]

Linear vs. Nonlinear

An ODE is linear if it can be written as the linear combination:

    \[ \frac{\mathrm{d}^ny}{\mathrm{d}x^n}+\sum_{i=0}^{n-1} a_i(x)\frac{\mathrm{d}^iy}{\mathrm{d}x^i}=f(x) \]

where a_i(x) and f(x) are functions of the independent variable x. All the examples above are considered linear ODEs.
The following are two examples of nonlinear ODEs with x being the dependent variable and t being the independent variable:

    \[ \begin{split} \frac{\mathrm{d}x}{\mathrm{d}t}&=\frac{k}{x}\\ \frac{\mathrm{d}^2x}{\mathrm{d}t^2}&=\sin{x} \end{split} \]

Homogeneous vs. Nonhomogeneous

A homogeneous ODE is an equation whose every term contains either the dependent variable or one of its derivatives. For example, Newton’s second law applied to a spring without the gravitational term is a linear homogeneous ODE:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}&=-kx\\ \]

If one of the non-zero terms in the ODE contains neither the dependent variable nor any of its derivatives, the equation is nonhomogeneous. As an example, if the gravitational term is added to Newton’s second law, then, the equation becomes nonhomogeneous:

    \[ m \frac{\mathrm{d}^2x}{\mathrm{d}t^2}&=mg-kx\\ \]

The free term (the term devoid of the dependent variables and any of its derivatives) is called the source term.

A linear homogeneous ODE has the following form:

    \[ \frac{\mathrm{d}^ny}{\mathrm{d}x^n}+\sum_{i=0}^{n-1} a_i(x)\frac{\mathrm{d}^iy}{\mathrm{d}x^i}=0 \]

while a linear nonhomogeneous ODE has the form:

    \[ \frac{\mathrm{d}^ny}{\mathrm{d}x^n}+\sum_{i=0}^{n-1} a_i(x)\frac{\mathrm{d}^iy}{\mathrm{d}x^i}=f(x) \]

The term f(x) is called the source term. The solution to a linear homogeneous equation is called the complementary solution y_c while the solution when a source term appears in the equation is the sum of the complementary solution y_c and the particular solution y_p (which is particular to the source term f(x)). For example, the solution to Newton’s second law without the gravitational term has the form:

    \[ x=C_1\cos\left(\frac{\sqrt{k}t}{\sqrt{m}}\right)+C_2\sin\left(\frac{\sqrt{k}t}{\sqrt{m}}\right) \]

while the solution when the gravitational term appears has the form:

    \[ x=\frac{mg}{k}+C_1\cos\left(\frac{\sqrt{k}t}{\sqrt{m}}\right)+C_2\sin\left(\frac{\sqrt{k}t}{\sqrt{m}}\right) \]

The particular solution is the term \frac{mg}{k}.

IVP vs. BVP

Depending on the boundary conditions, an ODE can be classified as either an Initial Value Problem (IVP) or a Boundary Value Problem (BVP). An initial value problem is an ODE given with initial conditions of the dependent variable and its derivative at a particular value of the independent variable. This usually applies to dynamic systems whose independent variable is time. For example, Newton’s second law of motion is an initial value problem because the initial value at t=0 of the displacement x(0) and the velocity x'(0) are required to reach a solution. The initial values lead to a particular path that is a function of time x(t).

A boundary value problem is an ODE given with boundary conditions at different points. This usually applies to static systems whose independent variable is position. For example, the Euler-Bernoulli beam deflection equation is a boundary value problem. The boundary conditions of the displacement y, the rotation y', and the third and fourth derivatives, are usually given at boundary points on the beam which will then dictate the equilibrium deflection of the beam as a function of position y(x).

Solution Methods for IVPs

(Explicit) Euler Method

The Euler method is one of the simplest methods for solving first-order IVPs. Consider the following IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) \]

Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, we can use a Taylor approximation to estimate the value of x at t=t_{i+1}, namely x(t_{i+1}) with h=t_{i+1}-t_i:

    \[ x(t_{i+1})=x_i+\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{t_i} (h) +\mathcal{O}(h^2) \]

Substituting the differential equation into the above equation yields:

    \[ x(t_{i+1})=x_i+F(x_i,t_i) (h) +\mathcal{O}(h^2) \]

Therefore, as an approximation, an estimate for x(t_{i+1}) can be taken as x_{i+1} as follows:

    \[ x(t_{i+1})\approx x_{i+1}=x_i+F(x_i,t_i) (h) \]

Using this estimate, the local truncation error is thus proportional to the square of the step size with the constant of proportionality related to the second derivative of x, which is the first derivative of the given IVP:

    \[ E_{\mbox{local}}=\mathcal{O}(h^2) \]

If the errors from each interval are added together, with n=\frac{L}{h} being the number of intervals and L the total length t_n-t_0, then, the total error is:

    \[ E=\mathcal{O}(h^2) \frac{L}{h}=\mathcal{O}(h)  \]

There are many examples in engineering and biology in which such IVPs appear. Check this page for examples.

The following Mathematica code provides a procedure whose inputs are

  • The differential equation as a function of the dependent variable x and the independent variable t.
  • The step size h.
  • The initial value t_0 and the final value t_{max} of the independent variable.
  • The initial value x(t_0)=x_0.

The procedure then carries on the Euler method and outputs the required data vector.
View Mathematica Code

EulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] = xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2,  n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)
function1[y_, t_] := 0.05 y
function2[x_, t_] := 0.015 x
EulerMethod[function1, 35, 1.0, 0, 50] // MatrixForm
EulerMethod[function2, 35, 1.0, 0, 50] // MatrixForm

Example 1

The Canadian population at t=0 (current year) is 35 million. If 2.5% of the population have a child in a given year, while the death rate is 1% of the population, what will the population be in 50 years? What if 6% of the population have a child in a given year and the death rate is kept constant at 1%, what will the population be in 50 years?

Solution

The rate of growth is directly proportional to the current population with the constant of proportionality k=0.025-0.01=0.015. If x represents the population in millions and t represents the time in years, then the IVP is given by:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=0.015 x \]

With the initial condition of x(0)=35, the exact solution to this differential equation is:

    \[ x=35e^{0.015t} \]

If the birth rate is 6%, then k=0.06-0.01=0.05 and the exact solution of the population follows the following equation:

    \[ x=35e^{0.05t} \]

The Euler method provides a numerical solution as follows. Using a step size of h=t_{i+1}-t_i=1, we can set t_0=0, t_1=1, t_2=2, \cdots, t_{50}=50 years. The value of x_0 is given as 35 million. When k=0.015, the value of x_1 in millions is given by:

    \[ x_1=x_0+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_0,t=t_0}=35+0.015\times 35=35.525 \]

The value of x_2 is given by:

    \[ x_2=x_1+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_1,t=t_1}=35.525+0.015\times 35.525=36.0579 \]

Iterating further, the values of x at each time point can be obtained. The following is the Microsoft Excel table showing the values of t_i in years and the corresponding values of the population x_i in millions:
IVPe1
When k=0.015, the numerical procedure produces a very good approximation for the population at t=50 years. The error between the prediction of the exact solution to the differential equation and the numerical value at t=50 is given in millions by:

    \[ E=35e^{(0.015\times 50)} - 73.6835=74.095-73.6835=0.4 \]

The relative error in this case is given by:

    \[ E_r=\frac{0.4}{74.095}=0.005 \]

For the second case, when k=0.05, the value of x_1 in millions is given by:

    \[ x_1=x_0+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_0,t=t_0}=35+0.05\times 35=36.75 \]

Proceeding iteratively, the following is the Microsoft Excel table showing the values of t_i in years and the corresponding values of the population x_i in millions:
IVPe2
Compared to the previous case, when k=0.05, the numerical procedure is not as good. The error between the prediction of the exact solution to the differential equation and the numerical value at t=50 is given in millions by:

    \[ E=35e^{(0.05\times 50)} - 401.359=426.387-401.359=25.0283 \]

The relative error in this case is given by:

    \[ E_r=\frac{25.0283}{426.387}=0.06 \]

The following is the graph showing the exact solution versus the data using the Euler method for both cases. The Mathematica code used is given below. Notice how the error in the Euler method increases as t increases.
IVP2
View Mathematica Code

EulerMethod[fp_, x0_, h_, t0_, tmax_] := 
  (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] =  xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

Clear[x, xtable]
k = 0.015;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]

Clear[x, xtable]
k = 0.05;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]

The following tool illustrates the effect of choosing the step size h on the difference between the numerical solution (shown as dots) and the exact solution (shown as the blue curve) for the case when k=0.05. The smaller the step size, the more accurate the solution is. The higher the value of h, the more the numerical solution deviates from the exact solution. The error increases away from the initial conditions.

Example 2

This example is adopted from this page. A lake of volume 10,000 m^3 has an initial concentration of a particular pollutant of 5 parts per billion (ppb). This pollutant is due to a pesticide that is no longer available in the market. The volume of the lake is constant throughout the year, with daily water flowing into and out of the lake at the rate of q=100+50\cos{(2\pi t/365)} \frac{m^3}{day}. As the pesticide is no longer used, the concentration of the pollutant in the surrounding soil is given in ppb as a function of t(days) by c_{in}=5e^{\left(\frac{-2t}{1000}\right)} which indicates that the concentration is following an exponential decay. This concentration can be assumed to be that of the pollutant in the fluid flowing into the lake. The concentration of the pollutant in the fluid flowing out of the lake is the same as the concentration in the lake. What is the concentration of this particular pollutant in the lake after two years?

Solution

The first step in this problem is to properly write the differential equation that needs to be solved. If \alpha is the total amount of pollutant in the lake, then, the rate of change of \alpha with respect to t is equal to the total amount of pollutant entering into the lake per day. Each day, the total amount of pollutant entering the lake is given by c_{in}\times q, while the total amount of pollutant exiting the lake is given by c\times q. Therefore:

    \[ \frac{\mathrm{d}\alpha}{\mathrm{d}t}=\left(100+50\cos{(2\pi t/365)}\right)\left(c_{in}-c\right) \]

The concentration is equal to the total amount divided by the volume of the lake, therefore, the differential equation in terms of c is given by:

    \[ \frac{\mathrm{d}c}{\mathrm{d}t}=\frac{\left(100+50\cos{(2\pi t/365)}\right)}{10000}\left(c_{in}-c\right)=\frac{\left(100+50\cos{(2\pi t/365)}\right)}{10000}\left(5e^{\left(\frac{-2t}{1000}\right)}-c\right) \]

Notice that finding an exact solution to the above differential equation is not an easy task and a numerical solution would be the preferred solution in this case. Taking h=t_{i+1}-t_i=1 day, t_0=0, and c_0=5 ppb the Euler method can be used to find the concentration at c=730 days. You can download the Microsoft Excel file example2.xlsx showing the calculations and the produced graph.
Taking h=1 day, we will show the calculations for the estimates c_1 and c_2. At t_1=1 day:

    \[\begin{split} c_1&=c_0+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_0,t=t_0}\\ &=5+1\left(\frac{\left(100+50\cos{(2\pi 0/365)}\right)}{10000}\left(5e^{\left(\frac{-2(0)}{1000}\right)}-5\right)\right)\\ &=5 \end{split} \]

At t_2=2 days:

    \[\begin{split} c_2&=c_1+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_1,t=t_1}\\ &=5+1\left(\frac{\left(100+50\cos{(2\pi 1/365)}\right)}{10000}\left(5e^{\left(\frac{-2(1)}{1000}\right)}-5\right)\right)\\ &=4.9999 \end{split} \]

Carrying on produces the values of c at the remaining time intervals up to t=2(365) days. The following is the output produced by the Mathematica code below.
IPV1

View Mathematica Code
EulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] =  xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)
f[t_] := 100 + 50*Cos[2 Pi*t/365];
cin[t_] := 5 E^(-2 t/1000);
fp[c_, t_] := f[t] (cin[t] - c)/10000;
Data = EulerMethod[fp, 5, 1.0, 0, 365*2];
ListPlot[Data, AxesLabel -> {"time (days)", "Concentration ppb"}]

Example 3

In this example we investigate the issue of stability of numerical solution of differential equations using the Euler method. Consider the IVP of the form:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=-x \]

If x_0=1, compare the two solutions taking a step size of h=0.1, and a step size of h=3.

Solution

The exact solution to the given equation is given by:

    \[ x=e^{-t} \]

Using an explicit Euler scheme, the value of x_{i+1} can be obtained as follows:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x_i}=(1-h)x_i=(1-h)(1-h)x_{i-1} \]

Therefore:

    \[ x_{i+1}=(1-h)^{i+1}x_0 \]

When h=0.1, the term (1-h)^{i+1}=0.9^{i+1} is bounded for any value of i which means that the value of x_{i+1} is always bounded. However, if h=3, then, the term (1-h)^{i+1}=(-2)^{i+1} oscillates wildly leading to an instability in the numerical scheme!
The following tool illustrates the difference between the obtained result for t\in[0,10]. Vary the value of h and try to identify the point above which the obtained solution (black dots) starts oscillating wildly around the exact solution (blue curve).


 

Backward (Implicit) Euler Method

Consider the following IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) \]

Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, we can use a Taylor approximation to relate the value of x at t=t_{i+1}, namely x(t_{i+1}) with h=t_{i+1}-t_i. However, unlike the explicit Euler method, we will use the Taylor series around the point x(t_{i+1}), that is:

    \[ x_i=x(t_{i+1})-\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{t_{i+1}} (h) +\mathcal{O}(h^2) \]

Substituting the differential equation into the above equation yields:

    \[ x_i=x(t_{i+1})-F(x_{i+1},t_{i+1}) (h) +\mathcal{O}(h^2) \]

Therefore, as an approximation, an estimate for x(t_{i+1}) can be taken as x_{i+1} as follows:

    \[ x(t_{i+1})\approx x_{i+1}=x_i+F(x_{i+1},t_{i+1}) (h) \]

Using this estimate, the local truncation error is thus proportional to the square of the step size with the constant of proportionality related to the second derivative of x, which is the first derivative of the given IVP. The backward Euler method is termed an “implicit” method because it uses the slope \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) at the unknown point x_{i+1}, namely: F(x_{i+1},t_{i+1}). The developed equation can be linear in x_{i+1} or nonlinear. Nonlinear equations can often be solved using the fixed-point iteration method or the Newton-Raphson method to find the value of x_{i+1}. While the implicit scheme does not provide better accuracy than the explicit scheme, it comes with additional computations. However, one advantage is that this scheme is always stable.

The following Mathematica code adopts the implicit Euler scheme and uses the built-in FindRoot function to solve for x_{i+1}. The code is similar to the code provided for the explicit scheme except for the line to calculate x_{i+1}.
View Mathematica Code

ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
 xtable = Table[0, {i, 1, n}];
 xtable[[1]] = x0;
 Do[  s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}];  xtable[[i]] = x /. s[[1]],  {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

function1[y_, t_] := 0.05 y
function2[x_, t_] := 0.015 x
ImplicitEulerMethod[function1, 35, 1.0, 0, 50] // MatrixForm
ImplicitEulerMethod[function2, 35, 1.0, 0, 50] // MatrixForm
Example 1

Repeat Example 1 above using the implicit Euler method.

Solution

When k=0.015, the IVP is given by:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=0.015 x \]

Given a value for x_i at t_i, and taking h=t_{i+1}-t_i=1, the estimate for x_{i+1} can be calculated according to the formula:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(0.015x_{i+1}) \]

Rearranging yields:

    \[ x_{i+1}=\frac{x_i}{1-0.015} \]

With x_0=35 million and t_0=0 years, the estimate for the population at t=1 years using the implicit Euler method is given by:

    \[ x_1=\frac{35}{1-0.015}=35.533 \]

at t=2 years, the value of x_2 in millions is given by:

    \[ x_2=\frac{35.533}{1-0.015}=36.074 \]

Proceeding iteratively, the following is the Excel table showing the values of the population up to t=50 years.
Implicit1
When k=0.015, the implicit method provides an error similar to that of the explicit method:

    \[\begin{split} E&=35e^{(0.015\times 50)}-74.517=74.095-74.517=-0.4\\ E_r&=\frac{-0.4}{74.095}=-0.005 \end{split} \]

When k=0.05, the IVP is given by:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=0.05 x \]

Given a value for x_i at t_i, and taking h=t_{i+1}-t_i=1, the estimate for x_{i+1} can be calculated according to the formula:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(0.05x_{i+1}) \]

Rearranging yields:

    \[ x_{i+1}=\frac{x_i}{1-0.05} \]

With x_0=35 million and t_0=0 years, the estimate for the population at t=1 years using the implicit Euler method is given by:

    \[ x_1=\frac{35}{1-0.05}=36.842 \]

at t=2 years, the value of x_2 in millions is given by:

    \[ x_2=\frac{36.842}{1-0.05}=38.781 \]

Proceeding iteratively, the following is the Excel table showing the values of the population up to t=50 years.
Implicit11
Again, the error produced by the implicit method is comparable to that produced by the explicit method shown above:

    \[\begin{split} E&=35e^{(0.05\times 50)}-454.871=426.387-454.871=-28.483\\ E_r&=\frac{-28.483}{426.387}=-0.07 \end{split} \]

The following is the graph showing the exact solution versus the data using both the implicit (blue dots) and the explicit (black dots) Euler method for both values of k.
Implicit12
View Mathematica Code

EulerMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xtable[[i]] = xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

ImplicitEulerMethod[fp_, x0_, h_, t0_,  tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

Clear[x, xtable]
k = 0.015;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Data2 = ImplicitEulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data], Blue, Point[Data2]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]

Clear[x, xtable]
k = 0.05;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := k*x;
Data = EulerMethod[fp, 35, 1.0, 0, 50];
Data2 = ImplicitEulerMethod[fp, 35, 1.0, 0, 50];
Plot[x, {t, 0, 50}, Epilog -> {PointSize[Large], Point[Data], Blue, Point[Data2]}, AxesOrigin -> {0, 0}, AxesLabel -> {"Time (years)", "Population (Millions)"}]

The following tool illustrates the effect of choosing the step size h on the difference between the numerical solution obtained using the implicit methods (shown as dots) and the exact solution (shown as the blue curve) for the case when k=0.05. The smaller the step size, the more accurate the solution is. The higher the value of h, the more the numerical solution deviates from the exact solution. The error increases away from the initial conditions.


 

Example 2

Repeat Example 2 above using the implicit Euler method.

Solution

Recall that the IVP for this problem is given by:

    \[ \frac{\mathrm{d}c}{\mathrm{d}t}=\frac{\left(100+50\cos{(2\pi t/365)}\right)}{10000}\left(5e^{\left(\frac{-2t}{1000}\right)}-c\right) \]

Taking h=1 day, c_{i+1} can be predicted for a given value of c_i using the implicit Euler method as follows:

    \[ c_{i+1}=c_i+h\frac{\mathrm{d}c}{\mathrm{d}t}\bigg|_{c=c_{i+1},t=t_{i+1}}=c_i+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}\left(5e^{\left(\frac{-2t_{i+1}}{1000}\right)}-c_{i+1}\right) \]

Rearranging:

    \[ c_{i+1}=\frac{c_i+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}\left(5e^{\left(\frac{-2t_{i+1}}{1000}\right)}\right)}{1+\frac{\left(100+50\cos{(2\pi t_{i+1}/365)}\right)}{10000}} \]

With t_0=0 and c_0=5 ppb, c_1 is calculated as:

    \[ c_1=\frac{5+0.015\left(5e^{\left(\frac{-2(1)}{1000}\right)\right)}}{1+0.015}=4.9999 \]

Proceeding iteratively produces the following graph with results almost identical to those produced by the explicit method. You can download the Microsoft Excel file example2implicit.xlsx showing the calculations.
implicitexample2

View Mathematica Code
ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] :=
 (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

f[t_] := 100 + 50*Cos[2 Pi*t/365];
cin[t_] := 5 E^(-2 t/1000);
fp[c_, t_] := f[t] (cin[t] - c)/10000;
Data = ImplicitEulerMethod[fp, 5, 1, 0, 365*2];
ListPlot[Data, AxesLabel -> {"time (days)", "Concentration ppb"}]
Example 3

We will repeat Example 3 above to illustrate that the implicit Euler method is always stable.
Consider the IVP of the form:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=-x \]

If x_0=1, compare the two solutions taking a step size of h=0.1, and a step size of h=3.

Solution

Using an implicit Euler scheme, the value of x_{i+1} can be obtained as follows:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x_{i+1}}=x_i-h x_{i+1}\Rightarrow x_{i+1}=\frac{x_i}{(1+h)}=\frac{x_{i-1}}{(1+h)(1+h)} \]

Therefore:

    \[ x_{i+1}=\frac{x_0}{(1+h)^{i+1}} \]

h is a positive number. Therefore, the term \frac{1}{(1+h)^{i+1}} is bounded for every i and therefore, x_{i+1} is bounded. The following tool illustrates the difference between the obtained result for t\in[0,10]. Vary the value of h and notice how the obtained solution (black dots) compares to the exact solution (blue curve). While the accuracy is lost with higher values of h, the solution maintains numerical stability.


 

Example 4

Consider the following nonlinear IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=tx^2+2x \]

with the initial condition x(0)=-5. Find the values of x for t\in[0,5] taking h=0.1 using the implicit Euler method.

Solution

Using Mathematica, the exact solution to the differential equation given can be obtained as:

    \[ x=\frac{-20e^{2t}}{9-5e^{2t}+10e^{2t}t} \]

The implicit Euler scheme provides the following estimate for x_{i+1}:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1},t=t_{i+1}}=x_i+h(t_{i+1}x_{i+1}^2+2x_{i+1}) \]

Since x_{i+1} appears on both sides, and the equation is nonlinear in x_{i+1}, therefore, the Newton-Raphson method will be used at each time increment to find the value of x_{i+1}!
Setting x_0=-5, h=0.1, t_0=0, and t_1=0.1 yields:

    \[ x_1=-5+0.1(0.1 x_1^2+2x_1) \]

Solving the above nonlinear equation in x_1 using the Newton-Raphson method yields:

    \[ x_1=-5.82576 \]

Similarly, with t_2=0.2, the equation for x_2 is:

    \[ x_2=-5.82576+0.1(0.2x_2^2+2x_2) \]

Using the Newton-Raphson method:

    \[ x_2=-6.29235 \]

Proceeding iteratively produces the following table:
implicit42
The following is the plot of the exact solution (blue line) with the data obtained using the implicit Euler method. The Mathematica code is given below.
Implicit4

View Mathematica Code
ImplicitEulerMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[s = FindRoot[x == xtable[[i - 1]] + h*fp[x, t0 + (i - 1) h], {x, xtable[[i - 1]]}]; xtable[[i]] = x /. s[[1]], {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = ImplicitEulerMethod[fp, -5.0, 0.1, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]},AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm

Heun’s Method

Heun’s method provides a slight modification to both the implicit and explicit Euler methods. Consider the following IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) \]

Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, we have:

    \[ \int_{x_i}^{x_{i+1}}\!\frac{\mathrm{d}x}{\mathrm{d}t}\,\mathrm{d}t=\int_{t_i}^{t_{i+1}}\!F(x,t)\,\mathrm{d}t \]

Therefore:

    \[ x_{i+1}-x_i=\int_{t_i}^{t_{i+1}}\!F(x,t)\,\mathrm{d}t \]

If the trapezoidal rule is used for the right-hand side with one interval we obtain:

    \[ x_{i+1}-x_i=h\frac{F(x_{i+1},t_{i+1})+F(x_{i},t_{i})}{2}+\mathcal{O}(h^3) \]

where the error term is inherited from the use of the trapezoidal integration rule. Therefore an estimate for x_{i+1} is given by:

    \[ x_{i+1}=x_i+h\frac{F(x_{i+1},t_{i+1})+F(x_{i},t_{i})}{2}+\mathcal{O}(h^3) \]

Using this estimate, the local truncation error is thus proportional to the cube of the step size. If the errors from each interval are added together, with n=\frac{L}{h} being the number of intervals and L the total length t_n-t_0, then, the total error is:

    \[ E=\mathcal{O}(h^3) \frac{L}{h}=\mathcal{O}(h^2)  \]

which is better than both the implicit and explicit Euler methods. Note that Heun’s method is essentially a slight modification to the Euler’s method in which the slope used to calculate the value of x_{i+1} at the next time point is used as the average slope at x_i and at x_{i+1}. Heun’s method can be implemented in two ways. One way is to use the slope at x_i to calculate an initial estimate x_{i+1}^{(0)}. Then, the estimate for x_{i+1} would be calculated based on the slopes at x_i and x_{i+1}^{(0)}. Alternatively, the Newton-Raphson method or the fixed-point iteration method can be used to solve directly for x_{i+1}. We will only adopt the first way. The following Mathematica code presents a procedure to utilize Heun’s method. The procedure is essentially similar to the ones presented in the Euler methods except for the step to calculate the new estimate x_{i+1}.
View Mathematica Code

HeunMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xi0 = xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h]; 
   xtable[[i]] = xtable[[i - 1]] + h (fp[xtable[[i - 1]], t0 + (i - 2) h] +fp[xi0, t0 + (i - 1) h])/2, {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

Example

Solve Example 4 above using Heun’s method.

Solution

Heun’s method is implemented by first assuming an estimate for x_{i+1}^{(0)} based on the explicit Euler method:

    \[ x_{i+1}^{(0)}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i},t=t_{i}} \]

Then, the estimate for x_{i+1} is calculated as:

    \[ x_{i+1}=x_i+h\frac{\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1}^{(0)},t=t_{i+1}}+\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i},t=t_{i}}}{2} \]

Setting x_0=-5, t_0=0, t_1=0.1, and h=0.1, an initial estimate for x_1^{(0)} is given by:

    \[ x_1^{(0)}=-5+h(t_0x_0^2+2x_0)=-5+0.1(0-10)=-6 \]

Using this information, the slope at x_1^{(0)} can be calculated and used to estimate x_1:

    \[ x_1=-5+\frac{h}{2}\left(t_0x_0^2+2x_0 +t_1x_1^{(0)}^2+2x_1^{(0)} \right)=-5.92 \]

Similarly, an initial estimate for x_2^{(0)} is given by:

    \[ x_2^{(0)}=-5.92+h(t_1x_1^2+2x_1)=-5.92+0.1(0.1(-5.92)^2-2\times 5.92)=-6.7535 \]

Using this information, the slope at x_2^{(0)} can be calculated and used to estimate x_2:

    \[ x_2=-5.92+\frac{h}{2}\left(t_1x_1^2+2x_1 +t_2x_2^{(0)}^2+2x_2^{(0)} \right)=-6.556 \]

Proceeding iteratively gives the values of x_i up to t=5. The Microsoft Excel file Heun.xlsx provides the required calculations.

The following graph shows the produced numerical data (black dots) overlapping the exact solution (blue line). The Mathematica code is given below.
H1

View Mathematica Code
HeunMethod[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xi0 = xtable[[i - 1]] + h*fp[xtable[[i - 1]], t0 + (i - 2) h]; 
   xtable[[i]] =  xtable[[i - 1]] +   h (fp[xtable[[i - 1]], t0 + (i - 2) h] +  fp[xi0, t0 + (i - 1) h])/2, {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)
Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = HeunMethod[fp, -5.0, 0.1, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]}, AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm

The following tool provides a comparison between the explicit Euler method and Heun’s method. Notice that around t\in[0,0.8], the function x varies highly in comparison to the rest of the domain. The biggest difference between Heun’s method and the Euler method can be seen when h=0.1 around this area. Heun’s method is able to trace the curve while the Euler method has higher deviations.


 

The Midpoint Method

Similar to Heun’s method, the midpoint method provides a slight modification to both the implicit and explicit Euler methods. Consider the following IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) \]

Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, we have:

    \[ \int_{x_i}^{x_{i+1}}\!\frac{\mathrm{d}x}{\mathrm{d}t}\,\mathrm{d}t=\int_{t_i}^{t_{i+1}}\!F(x,t)\,\mathrm{d}t \]

Therefore:

    \[ x_{i+1}-x_i=\int_{t_i}^{t_{i+1}}\!F(x,t)\,\mathrm{d}t \]

If the midpoint rule of the rectangle method of numerical integration is used for the right-hand side with one interval we obtain:

    \[ x_{i+1}-x_i=hF\left(\frac{x_{i+1}+x_i}{2},\frac{t_{i+1}+t_{i}}{2}\right)+\mathcal{O}(h^3) \]

where the error term is inherited from the use of I_2 in the rectangle rule. Therefore an estimate for x_{i+1} is given by:

    \[ x_{i+1}=x_i+hF\left(\frac{x_{i+1}+x_i}{2},\frac{t_{i+1}+t_{i}}{2}\right)+\mathcal{O}(h^3) \]

Using this estimate, the local truncation error is thus proportional to the cube of the step size. If the errors from each interval are added together, with n=\frac{L}{h} being the number of intervals and L the total length t_n-t_0, then, the total error is:

    \[ E=\mathcal{O}(h^3) \frac{L}{h}=\mathcal{O}(h^2)  \]

which is similar to Heun’s method but better than both the implicit and explicit Euler methods. Note that the midpoint method is essentially a slight modification to the Euler’s method in which the slope used to calculate the value of x_{i+1} at the next time point is used as the slope at the average point x_{i+1/2}=\frac{x_{i+1}+x_i}{2}. The midpoint method can be implemented in two ways. One way is to use the slope at x_i to calculate an initial estimate x_{i+1/2}^{(0)}. Then, the estimate for x_{i+1} would be calculated based on the slope at x_{i+1/2}^{(0)}. Alternatively, the Newton-Raphson method or the fixed-point iteration method can be used to solve directly for x_{i+1}. We will only adopt the first way. The following Mathematica code presents a procedure to utilize the midpoint method. The procedure is essentially similar to the ones presented in the Euler methods except for the step to calculate the new estimate x_{i+1}.
View Mathematica Code

MidPointMethod[fp_, x0_, h_, t0_, tmax_] := 
(n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xi120 = xtable[[i - 1]] + h/2*fp[xtable[[i - 1]], t0 + (i - 2) h]; 
   xtable[[i]] = xtable[[i - 1]] + h (fp[xi120, t0 + (i - 1-1/2) h]), {i, 2, n}];
  Data = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data)

Example

Solve Example 4 above using the midpoint method.

Solution

The midpoint method is implemented by first assuming an estimate for x_{i+1/2}^{(0)} based on the explicit Euler method:

    \[ x_{i+1/2}^{(0)}=x_i+\frac{h}{2}\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i},t=t_{i}} \]

Then, the estimate for x_{i+1} is calculated as:

    \[ x_{i+1}=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{x=x_{i+1/2}^{(0)},t=t_{i+1/2}} \]

Setting x_0=-5, t_0=0, t_1=0.1, and h=0.1, an initial estimate for x_{0.5}^{(0)} is given by:

    \[ x_{0.5}^{(0)}=-5+\frac{h}{2}(t_0x_0^2+2x_0)=-5+\frac{0.1}{2}(0-10)=-5.5 \]

Using this information, the slope at x_{0.5}^{(0)} can be calculated and used to estimate x_1:

    \[ x_1=-5+h\left(t_{0.5}x_{0.5}^{(0)}^2+2x_{0.5}^{(0)} \right)=-5.949 \]

Similarly, an initial estimate for x_{1.5}^{(0)} is given by:

    \[ x_{1.5}^{(0)}=-5.949+\frac{h}{2}(t_1x_1^2+2x_1)=-5.949+0.05(0.1(-5.949)^2-2\times 5.949)=-6.3667 \]

Using this information, the slope at x_{1.5}^{(0)} can be calculated and used to estimate x_2:

    \[ x_2=-5.949+h\left(t_{1.5}x_{1.5}^{(0)}^2+2x_{1.5}^{(0)}\right)=-6.614 \]

Proceeding iteratively gives the values of x_i up to t=5. The Microsoft Excel file MidPoint.xlsx provides the required calculations.

The following graph shows the produced numerical data (black dots) overlapping the exact solution (blue line). The Mathematica code is given below.
m1

View Mathematica Code
MidPointMethod[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[xi120 = xtable[[i - 1]] + h/2*fp[xtable[[i - 1]], t0 + (i - 2) h];
   xtable[[i]] =  xtable[[i - 1]] + h (fp[xi120, t0 + (i - 1 - 1/2) h]), {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)
Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = MidPointMethod[fp, -5.0, 0.1, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]}, AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm

The following tool provides a comparison between the explicit Euler method and the midpoint method. The midpoint method behaves very similar to Heun’s method. Notice that around t\in[0,0.8], the function x varies highly in comparison to the rest of the domain. The biggest difference between midpoint method and the Euler method can be seen when h=0.1 around this area. The midpoint method is able to trace the curve while the Euler method has higher deviations.


 

Runge-Kutta Methods

The Runge-Kutta methods developed by the German mathematicians C. Runge and M.W. Kutta are essentially a generalization of all the previous methods. Consider the following IVP:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t) \]

Assuming that the value of the dependent variable x (say x_i) is known at an initial value t_i, then, the Runge-Kutta methods employ the following general equation to calculate the value of x at t=t_{i+1}, namely x_{i+1} with h=t_{i+1}-t_i:

    \[ x_{i+1}=x_i+h\phi \]

where \phi is a function of the form:

    \[ \phi=\alpha_1 k_1+\alpha_2k_2+\alpha_3k_3+\cdots+\alpha_nk_n \]

\alpha_j are constants while k_j=F(x,t) evaluated at points within the interval [x_i,x_{i+1}] and have the form:

    \[\begin{split} k_1&=F(x_i,t_i)\\ k_2&=F(x_i+q_{11}k_1h,t_i+p_1h)\\ k_3&=F(x_i+q_{21}k_1h+q_{22}k_2h,t_i+p_2h)\\ &\vdots\\ k_n&=F(x_i+q_{(n-1),1}k_1h+q_{(n-1),2}k_2h+\cdots+q_{(n-1),(n-1)}k_{n-1}h,t_i+p_{n-1}h) \end{split} \]

The constants \alpha_j, and the forms of k_j are obtained by equating the value of x_{i+1} obtained using the Runge-Kutta equation to a particular form of the Taylor series. The k's are recurrence relationships, meaning k_1 appears in k_2, which appears in k_3, and so forth. This makes the method efficient for computer calculations. The error in a particular form depends on how many terms are used. The general forms of these Runge-Kutta methods could be implicit or explicit. For example, the explicit Euler method is a Runge-Kutta method with only one term with \alpha_1=1 and k_1=F(x_i,t_i). Heun’s method on the other hand is a Runge-Kutta method with the following non-zero terms:

    \[\begin{split} \alpha_1&=\alpha_2=\frac{1}{2}\\ k_1&=F(x_i,t_i)\\ k_2&=F(x_i+hk_1,t_i+h) \end{split} \]

Similarly, the midpoint method is a Runge-Kutta method with the following non-zero terms:

    \[\begin{split} \alpha_2&=1\\ k_1&=F(x_i,t_i)\\ k_2&=F(x_i+\frac{h}{2}k_1,t_i+\frac{h}{2}) \end{split} \]

The most popular Runge-Kutta method is referred to as the “classical Runge-Kutta method”, or the “fourth-order Runge-Kutta method”. It has the following form:

    \[ x_{i+1}=x_i+\frac{h}{6}\left(k_1+2k_2+2k_3+k_4\right) \]

with

    \[\begin{split} k_1&=F(x_i,t_i)\\ k_2&=F\left(x_i+\frac{h}{2}k_1,t_i+\frac{h}{2}\right)\\ k_3&=F\left(x_i+\frac{h}{2}k_2,t_i+\frac{h}{2}\right)\\ k_4&=F(x_i+hk_3,t_i+h) \end{split} \]

The error in the fourth-order Runge-Kutta method is similar to that of the Simpson’s 1/3 rule, with a local error term E_{\mbox{local}}=\mathcal{O}(h^5) which would translate into a global error term of E=\mathcal{O}(h^4). The method is termed fourth-order because the error term is directly proportional to the step size raised to the power of 4.

View Mathematica Code
RK4Method[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[k1 = fp[xtable[[i - 1]], t0 + (i - 2) h];
   k2 = fp[xtable[[i - 1]] + h/2*k1, t0 + (i - 1.5) h];
   k3 = fp[xtable[[i - 1]] + h/2*k2, t0 + (i - 1.5) h];
   k4 = fp[xtable[[i - 1]] + h*k3, t0 + (i - 1) h];
   xtable[[i]] = xtable[[i - 1]] + h (k1 + 2 k2 + 2 k3 + k4)/6, {i, 2,
     n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)

Example

Solve Example 4 above using the classical Runge-Kutta method but with h=0.4.

Solution

Recall the differential equation of this problem:

    \[ \frac{\mathrm{d}x}{\mathrm{d}t}=F(x,t)=tx^2+2x \]

Setting x_0=-5, t_0=0, h=0.4, the following are the values of k_1, k_2, k_3, and k_4 required to calculate x_1:

    \[\begin{split} k_1&=F(x_0,t_0)=F(-5,0)\\ &=0+2(-5)=-10\\ k_2&=F\left(x_0+\frac{h}{2}k_1,t_0+\frac{h}{2}\right)=F(-5+0.2(-10),0+0.2)\\ &=0.2(-5+0.2(-10))^2+2\times (-5+0.2(-10))=-4.2\\ k_3&=F\left(x_0+\frac{h}{2}k_2,t_0+\frac{h}{2}\right)=F(-5+0.2(-4.2),0+0.2)\\ &=0.2(-5+0.2(-4.2))^2+2\times (-5+0.2(-4.2))=-4.8589\\ k_4&=F(x_0+hk_3,t_0+h)=F(-5+0.4(-4.8589),0+0.4)\\ &=0.4(-5+0.4(-4.8589))^2+2\times(-5+0.4(-4.8589))=5.3981 \end{split} \]

Therefore:

    \[ x_1=x_0+\frac{h}{6}(k_1+2k_2+2k_3+k_4)=-5+\frac{0.4}{6}(-10-2\times 4.2-2\times 4.8589+5.3981)=-6.51465 \]

Proceeding iteratively gives the values of x_i up to t=5.2. The Microsoft Excel file RK4.xlsx provides the required calculations. The following Microsoft Excel table shows the required calculations:

RK42

The following Mathematica code implements the classical Runge-Kutta method for this problem with h=0.3. The output curve is shown underneath.
View Mathematica Code

RK4Method[fp_, x0_, h_, t0_, tmax_] := (n = (tmax - t0)/h + 1;
  xtable = Table[0, {i, 1, n}];
  xtable[[1]] = x0;
  Do[k1 = fp[xtable[[i - 1]], t0 + (i - 2) h];
   k2 = fp[xtable[[i - 1]] + h/2*k1, t0 + (i - 1.5) h];
   k3 = fp[xtable[[i - 1]] + h/2*k2, t0 + (i - 1.5) h];
   k4 = fp[xtable[[i - 1]] + h*k3, t0 + (i - 1) h];
   xtable[[i]] = xtable[[i - 1]] + h (k1 + 2 k2 + 2 k3 + k4)/6, {i, 2,
     n}];
  Data2 = Table[{t0 + (i - 1)*h, xtable[[i]]}, {i, 1, n}];
  Data2)
Clear[x, xtable]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
fp[x_, t_] := t*x^2 + 2*x;
Data2 = RK4Method[fp, -5.0, 0.3, 0, 5];
Plot[x, {t, 0, 5}, Epilog -> {PointSize[Large], Point[Data2]},  AxesOrigin -> {0, 0}, AxesLabel -> {"t", "x"}, PlotRange -> All]
Title = {"t_i", "x_i"};
Data2 = Prepend[Data2, Title];
Data2 // MatrixForm
RK41

The following tool provides a comparison between the explicit Euler method and the classical Runge-Kutta method. The classical Runge-Kutta method gives excellent predictions even when h=0.5, at which point the explicit Euler method fails to predict anything close to the exact equation. The classical Runge-Kutta method also provides better estimates than the midpoint method and Heun’s method above.

Solving Systems of IVPs

Many engineering problems contain systems of IVPs in which more than one dependent variable is a function of one independent variable, say t. For example, chemical reactions of more than one component are usually described using systems of IVPs in which the initial conditions are known. Similarly, the interaction of biological growth of different species (predators and prey) is usually represented by a system of IVPs in which the initial conditions are given. Generally, these systems can be represented as:

    \[\begin{split} \frac{\mathrm{d}x_1}{\mathrm{d}t}&=F_1(x_1,x_2,\cdots,x_n,t)\\ \frac{\mathrm{d}x_2}{\mathrm{d}t}&=F_2(x_1,x_2,\cdots,x_n,t)\\ &\vdots\\ \frac{\mathrm{d}x_n}{\mathrm{d}t}&=F_n(x_1,x_2,\cdots,x_n,t) \end{split} \]

The independent variable is discretized such that t_0<t_1<t_2<\cdots<t_m with a constant step h=t_{i+1}-t_i. In such systems, the initial conditions, namely the values of x_1(t_0), x_2(t_0), x_3(t_0),\cdots,x_n(t_0), are given. The values of x_1,x_2,\cdots,x_n when t=t_{i} are denoted: x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)}. In this case, the values of the dependent variables at t=t_{i+1} can be obtained using the following system of equations:

    \[ \begin{split} x_1^{(i+1)}&=x_1^{(i)}+\phi_1h\\ x_2^{(i+1)}&=x_2^{(i)}+\phi_2h\\ &\vdots\\ x_n^{(i+1)}&=x_n^{(i)}+\phi_nh\\ \end{split} \]

where \phi_1,\phi_2,\cdots,\phi_n depend on the method used. For example, if the explicit Euler method is used, then:

    \[\begin{split} \phi_1&=F_1(x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)},t_i)\\ \phi_2&=F_2(x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)},t_i)\\ &\vdots\\ \phi_n&=F_n(x_1^{(i)},x_2^{(i)},\cdots,x_n^{(i)},t_i) \end{split} \]

Alternatively, if the classical Runge-Kutta method is used, then:

    \[\begin{split} \phi_1&=\frac{1}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ \phi_2&=\frac{1}{6}(k_{21}+2k_{22}+2k_{23}+k_{24})\\ &\vdots\\ \phi_n&=\frac{1}{6}(k_{n1}+2k_{n2}+2k_{n3}+k_{n4}) \end{split} \]

where k_{j1}, k_{j2}, k_{j3}, and k_{j4} are the k values associated with the dependent variable x_j and calculated according to the classical Runge-Kutta technique.

It should also be noted that higher-order IVPs can be solved by converting them into a system of first-order IVPs. For example, consider the IVP:

    \[ \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=F(x,t) \]

with the initial conditions x(0) and x'(0) given. Then, setting y=\frac{\mathrm{d}x}{\mathrm{d}t}, the equation can be written as a system of two IVPs:

    \[\begin{split} \frac{\mathrm{d}x}{\mathrm{d}t}&=y\\ \frac{\mathrm{d}y}{\mathrm{d}t}&=F(x,t) \end{split} \]

with initial conditions given for x and y.

Example

Consider the damped mechanical system with k=m=1, g=10, \beta=0.15 with the initial conditions x(0)=x'(0)=0. Find the numerical solution for the position x(t) and the velocity x'(t) for t\in[0,2\pi] using the explicit Euler method and the Runge-Kutta method. Use h=0.1.

Solution

The differential equation of the damped system is given by:

    \[ \frac{\mathrm{d}^2x}{\mathrm{d}t^2}=10-x-0.15\frac{\mathrm{d}x}{\mathrm{d}t} \]

The following exact solution can be obtained using the built-in DSolve function in Mathematica:

    \[\begin{split} x&=10e^{(-0.075t)}\left(e^{(0.075t)}-\cos{(0.997184t)-0.0752118\sin{(0.997184t)}}\right)\\ x'&=10.0282e^{(-0.075t)}\sin{(0.997184t)} \end{split} \]

To find the numerical solution, the equation will be converted into a system of 2 IVPs. Setting y=\frac{\mathrm{d}x}{\mathrm{d}t}, the equation can be converted into the following system:

    \[\begin{split} \frac{\mathrm{d}x}{\mathrm{d}t}&=y\\ \frac{\mathrm{d}y}{\mathrm{d}t}&=10-x-0.15y\\ \end{split} \]

with the initial conditions y(0)=x(0)=0. With h=0.1, the interval can be split into 2\pi/0.1=63 time steps. The time discretization will be such that t_0=0, t_1=0.1, t_2=0.2, up to t_{63}=6.3. The corresponding positions are given by: x_0, x_1, \cdots,x_{63} while the corresponding velocities are given by: y_0, y_1, \cdots,y_{63}.

Explicit Euler Method

Using the explicit Euler method, the values of x_{i+1} and y_{i+1} can be evaluated as:

    \[\begin{split} x_{i+1}&=x_i+h\frac{\mathrm{d}x}{\mathrm{d}t}\bigg|_{(x_i,y_i,t_i)}=x_i+hy_i\\ y_{i+1}&=y_i+h\frac{\mathrm{d}y}{\mathrm{d}t}\bigg|_{(x_i,y_i,t_i)}=y_i+h(10-x_i-0.15y_i) \end{split} \]

With the initial conditions x_0=y_0=0, the values of x_1 and y_1 at t_1=0.1 can be calculated as:

    \[\begin{split} x_1&=x_0+hy_0=0+0=0\\ y_1&=y_0+h(10-x_0-0.15y_0)=0+0.1(10-0-0)=1 \end{split} \]

Similarly, the values of x_2 and y_2 can be obtained as:

    \[\begin{split} x_2&=x_1+hy_1=0+0.1(1)=0.1\\ y_2&=y_1+h(10-x_1-0.15y_1)=1+0.1(10-0-0.15\times 1)=1.985 \end{split} \]

Proceeding iteratively, the values of x and y=x' up to t=6.3 can be obtained. The Microsoft Excel file EulerSystem.xlsx shows the obtained values.

Classical Runge-Kutta method

The two differential equations are designated F_1 and F_2 as follows:

    \[\begin{split} F_1(x,y,t)&=\frac{\mathrm{d}x}{\mathrm{d}t}=y\\ F_2(x,y,t)&=\frac{\mathrm{d}y}{\mathrm{d}t}=10-x-0.15y \end{split} \]

Using the classical Runge-Kutta method, the values of x_{i+1} and y_{i+1} can be evaluated as:

    \[\begin{split} x_{i+1}&=x_i+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ y_{i+1}&=y_i+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24}) \end{split} \]

where:

    \[\begin{split} k_{11}&=F_1(x_i,y_i,t_i)=y_i\\ k_{21}&=F_2(x_i,y_i,t_i)=10-x_i-0.15y_i\\ k_{12}&=F_1\left(x_i+\frac{h}{2}k_{11},y_i+\frac{h}{2}k_{21},t_i+\frac{h}{2}\right)=y_i+\frac{h}{2}k_{21}\\ k_{22}&=F_2\left(x_i+\frac{h}{2}k_{11},y_i+\frac{h}{2}k_{21},t_i+\frac{h}{2}\right)=10-(x_i+\frac{h}{2}k_{11})-0.15(y_i+\frac{h}{2}k_{21})\\ k_{13}&=F_1\left(x_i+\frac{h}{2}k_{12},y_i+\frac{h}{2}k_{22},t_i+\frac{h}{2}\right)=y_i+\frac{h}{2}k_{22}\\ k_{23}&=F_2\left(x_i+\frac{h}{2}k_{12},y_i+\frac{h}{2}k_{22},t_i+\frac{h}{2}\right)=10-(x_i+\frac{h}{2}k_{12})-0.15(y_i+\frac{h}{2}k_{22})\\ k_{14}&=F_1(x_i+hk_{13},y_i+hk_{23},t_i+h)=y_i+hk_{23}\\ k_{24}&=F_2(x_i+hk_{13},y_i+hk_{23},t_i+h)=10-(x_i+hk_{13})-0.15(y_i+hk_{23}) \end{split} \]

With the initial conditions x_0=y_0=0, the values of k_{ij} can be calculated as:

    \[\begin{split} k_{11}&=y_0=0\\ k_{21}&=10-x_0-0.15y_0=10\\ k_{12}&=y_0+\frac{h}{2}k_{21}=0+0.05\times 10=0.5\\ k_{22}&=10-(x_0+\frac{h}{2}k_{11})-0.15(y_0+\frac{h}{2}k_{21})\\ &=10-(0+0.05\times 0)-0.15(0+0.05\times 10)=9.925\\ k_{13}&=y_0+\frac{h}{2}k_{22}=0+0.05(9.925)=0.49625\\ k_{23}&=10-(x_0+\frac{h}{2}k_{12})-0.15(y_0+\frac{h}{2}k_{22})\\ &=10-(0+0.05\times 0.5)-0.15(0+0.05\times 9.925)=9.90056\\ k_{14}&=y_0+hk_{23}=0.1\times 9.90056=0.990056\\ k_{24}&=10-(x_0+hk_{13})-0.15(y_0+hk_{23})\\ &=10-(0+0.1\times 0.49625)-0.15(0+0.1\times 9.90056)=9.80187 \end{split} \]

Therefore:

    \[\begin{split} x_1&=x_0+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ &=0+\frac{0.1}{6}(0+2\times 0.5+2\times 0.49625+0.990056)=0.049709\\ y_1&=y_0+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24})\\ &=0+\frac{0.1}{6}(10+2\times 9.925+2\times 9.90056+9.80187)=0.990883 \end{split} \]

For x_2 and y_2, the values of k_{ij} can be calculated as:

    \[\begin{split} k_{11}&=y_1=0.990883\\ k_{21}&=10-x_1-0.15y_1=10-0.049709-0.15\times 0.990883=9.80166\\ k_{12}&=y_1+\frac{h}{2}k_{21}=0.990883+0.05\times 9.80166=1.480966\\ k_{22}&=10-(x_1+\frac{h}{2}k_{11})-0.15(y_1+\frac{h}{2}k_{21})\\ &=10-(0.049709+0.05\times 0.990883)-0.15(0.990883+0.05\times 9.80166)=9.678902\\ k_{13}&=y_1+\frac{h}{2}k_{22}=0.990883+0.05(9.678902)=1.47482\\ k_{23}&=10-(x_1+\frac{h}{2}k_{12})-0.15(y_1+\frac{h}{2}k_{22})\\ &=10-(0.049709+0.05\times 1.480966)-0.15(0.990883+0.05\times 9.678902)=9.65502\\ k_{14}&=y_1+hk_{23}=0.990883+0.1\times 9.65502=1.956385\\ k_{24}&=10-(x_1+hk_{13})-0.15(y_1+hk_{23})\\ &=10-(0.049709+0.1\times 1.47482)-0.15(0.990883+0.1\times 9.65502)=9.50935 \end{split} \]

Therefore:

    \[\begin{split} x_2&=x_1+\frac{h}{6}(k_{11}+2k_{12}+2k_{13}+k_{14})\\ &=0.049709+\frac{0.1}{6}(0.990883+2\times 1.480966+2\times 1.47482+1.956385)=0.197356\\ y_2&=y_1+\frac{h}{6}(k_{21}+2k_{22}+2k_{23}+k_{24})\\ &=0.990883+\frac{0.1}{6}(9.80166+2\times 9.678902+2\times 9.65502 +9.50935)=1.95719 \end{split} \]

Proceeding iteratively, the values of x and y=x' up to t=6.3 can be obtained. The Microsoft Excel file RK4System.xlsx shows the calculated up to the required time point.

Comparison

The following figure shows the exact solution (blue curve) and the numerical solution (black dots) for the position x and the velocity x' obtained using the explicit Euler method and the classical Runge-Kutta method. Using the explicit Euler method, the numerical solution deviates from the exact solution for higher values of t. However, the classical Runge-Kutta method provides very accurate predictions. The Mathematica code utilized to produce the curves is shown below.
systemsIVP

View Mathematica Code
Clear[x]
a = DSolve[{x''[t] == 10 - x[t] - 0.15 x'[t], x'[0] == 0, x[0] == 0}, x, t]
x = x[t] /. a[[1]]
xp = Simplify[D[x, t]]

Plot[x, {t, 0, 2 Pi}, AxesOrigin -> {0, 0}, AxesOrigin -> {0, 0}, AxesLabel -> {"time", "x (position)"}]
Plot[xp, {t, 0, 2 Pi}, AxesOrigin -> {0, 0}, AxesOrigin -> {0, 0}, AxesLabel -> {"time", "x' (velocity)"}]

EulerMethod2[fp1_, fp2_, x10_, x20_, h_, t0_, tmax_] := (
  n = (tmax - t0)/h + 1;
  x1table = Table[0, {i, 1, n}];
  x2table = Table[0, {i, 1, n}];
  x1table[[1]] = x10;
  x2table[[1]] = x20;
  Do[x1table[[i]] = x1table[[i - 1]] +   h*fp1[x1table[[i - 1]], x2table[[i - 1]], t0 + (i - 2) h];
   x2table[[i]] = x2table[[i - 1]] +  h*fp2[x1table[[i - 1]], x2table[[i - 1]], t0 + (i - 2) h], 
    {i, 2,     n}];
  Data = Table[{t0 + (i - 1)*h, x1table[[i]], x2table[[i]]}, {i, 1,  n}];
  Data)

RK4Method2[fp1_, fp2_, x10_, x20_, h_, t0_,   tmax_] := (
  n = (tmax - t0)/h + 1;
  x1table = Table[0, {i, 1, n}];
  x1table[[1]] = x10;
  x2table = Table[0, {i, 1, n}];
  x2table[[1]] = x20;
  Do[k11 = fp1[x1table[[i - 1]], x2table[[i - 1]], t0 + (i - 2) h];
   k21 = fp2[x1table[[i - 1]], x2table[[i - 1]], t0 + (i - 2) h];
   k12 = fp1[x1table[[i - 1]] + h/2*k11, x2table[[i - 1]] + h/2*k21,      t0 + (i - 1.5) h];
   k22 = fp2[x1table[[i - 1]] + h/2*k11, x2table[[i - 1]] + h/2*k21,      t0 + (i - 1.5) h];
   k13 = fp1[x1table[[i - 1]] + h/2*k12, x2table[[i - 1]] + h/2*k22,      t0 + (i - 1.5) h];
   k23 = fp2[x1table[[i - 1]] + h/2*k12, x2table[[i - 1]] + h/2*k22,      t0 + (i - 1.5) h];
   k14 = fp1[x1table[[i - 1]] + h*k13, x2table[[i - 1]] + h*k23,      t0 + (i - 1) h];
   k24 = fp2[x1table[[i - 1]] + h*k13, x2table[[i - 1]] + h*k23,      t0 + (i - 1) h];
   x1table[[i]] = x1table[[i - 1]] + h (k11 + 2 k12 + 2 k13 + k14)/6;
   x2table[[i]] = x2table[[i - 1]] + h (k21 + 2 k22 + 2 k23 + k24)/6, {i, 2, n}];
  Data2 = Table[{t0 + (i - 1)*h, x1table[[i]], x2table[[i]]}, {i, 1, n}];
  Data2)


fp1[x_, y_, t_] := y;
fp2[x_, y_, t_] := 10 - x - 0.15 y;
Data = EulerMethod2[fp1, fp2, 0, 0, 0.1, 0, 2 Pi];
DataPosition = Drop[Data, None, {3}];
DataVelocity = Drop[Data, None, {2}];
a1 = Plot[x, {t, 0, 6.3},   Epilog -> {PointSize[Large], Point[DataPosition]},    AxesOrigin -> {0, 0}, ImageSize -> Medium,    AxesLabel -> {"Time", "x (position)"}];
a2 = Plot[xp, {t, 0, 6.3},    Epilog -> {PointSize[Large], Point[DataVelocity]},    AxesOrigin -> {0, 0}, ImageSize -> Medium,    AxesLabel -> {"Time", "x' (velocity)"}];
a = Grid[{{a1, a2}}];
Grid[{{"Explicit Euler"}, {a}}]
Data = RK4Method2[fp1, fp2, 0, 0, 0.1, 0, 2 Pi];
DataPosition = Drop[Data, None, {3}];
DataVelocity = Drop[Data, None, {2}];
a1 = Plot[x, {t, 0, 6.3},    Epilog -> {PointSize[Large], Point[DataPosition]},    ImageSize -> Medium, AxesOrigin -> {0, 0},    AxesLabel -> {"Time", "x (position)"}];
a2 = Plot[xp, {t, 0, 6.3},    Epilog -> {PointSize[Large], Point[DataVelocity]},    ImageSize -> Medium, AxesOrigin -> {0, 0},    AxesLabel -> {"Time", "x' (velocity)"}];
a = Grid[{{a1, a2}}];
Grid[{{"Classical Runge-Kutta"}, {a}}]

The effect of changing the value of h on the resulting numerical solution for the system of IVPs given here can be illustrated using the following tool. Even for h=0.5 the classical Runge-Kutta produces very accurate predictions at which point the explicit Euler method produces very erroneous results.


 

Using Mathematica

The built-in NDSolve in Mathematica uses the Runge-Kutta methods to solve a system of IVPs and produces piecewise interpolating polynomials for the dependent variables. The following code solves some of the preceding examples and compares with the exact solution. Note that in those examples when the exact and numerical solutions are plotted together, they overlap!

View Mathematica Code
Clear[x, xn, xtable]
(*Example 1*)
k = 0.015;
a = DSolve[{x'[t] == k*x[t], x[0] == 35}, x, t];
x = x[t] /. a[[1]]

b = NDSolve[{xn'[t] == k*xn[t], xn[0] == 35}, xn, {t, 0, 50}]
xn = xn /. b[[1]]
Plot[{x, xn[t]}, {t, 0, 50}, PlotLegends -> {"Exact", "Numerical using Mathematica"}, AxesLabel -> {"t", "Population"}]

(*Example 2*)
Clear[cn]
b = NDSolve[{cn'[t] == (100 + 50*Cos[2 Pi*t/365])/10000*(5 Exp[-2 t/1000] - cn[t]), cn[0] == 5}, cn, {t, 0, 700}]
cn = cn /. b[[1]]
Plot[cn[t], {t, 0, 700}, AxesLabel -> {"t", "Concentration"}]


(*Example 4*)

Clear[xn, x]
a = DSolve[{x'[t] == t*x[t]^2 + 2*x[t], x[0] == -5}, x, t];
x = x[t] /. a[[1]]
b = NDSolve[{xn'[t] == t*xn[t]^2 + 2 xn[t], xn[0] == -5},  xn, {t, 0, 5}]
xn = xn /. b[[1]]
Plot[{x, xn[t]}, {t, 0, 5}, PlotRange -> All, PlotLegends -> {"Exact", "Numerical using Mathematica"}, AxesLabel -> {"t", "x"}]

(*Example of a system*)
Clear[xn, yn, x]
a = DSolve[{x''[t] == 10 - x[t] - 0.15 x'[t], x'[0] == 0, x[0] == 0},   x, t]
x = x[t] /. a[[1]]
xp = Simplify[D[x, t]]
b = NDSolve[{xn'[t] == yn[t], yn'[t] == 10 - xn[t] - 0.15 yn[t],    xn[0] == 0, yn[0] == 0}, {xn, yn}, {t, 0, 2 Pi}]
xn = xn /. b[[1]]
yn = yn /. b[[1]]

Plot[{x, xn[t]}, {t, 0, 2 Pi}, AxesOrigin -> {0, 0},  AxesOrigin -> {0, 0},  PlotLegends -> {"Exact", "Numerical using Mathematica"},  AxesLabel -> {"time", "x (position)"}]
Plot[{xp, yn[t]}, {t, 0, 2 Pi}, AxesOrigin -> {0, 0},  AxesOrigin -> {0, 0},  PlotLegends -> {"Exact", "Numerical using Mathematica"},  AxesLabel -> {"time", "x' (velocity)"}]

Solution Methods for BVPs

As mentioned above, BVP are ordinary differential equations whose boundary conditions are given at different values of the independent variable. Usually, these points are the extreme points of the region of interest, or the “boundaries” of the domain and hence the term: “boundary value problem”. Consider a BVP with y as the dependent variable and x as the independent variable. The boundary conditions would be given in terms of y or the derivatives of y at different values of the independent variable x. When the boundary conditions are given in terms of y, then, this type of boundary conditions is termed Dirichlet boundary condition. Otherwise, when the boundary conditions are given in terms of the derivatives of y, this type of boundary conditions is termed Neumann boundary condition. A BVP is usually an ODE that is at least of the second-order as equations of the first order require only one boundary condition which is usually an initial boundary condition. There are various methods to solve boundary value problems; among them are the following:

  1. Shooting method. In the shooting method the boundary value problem is solved as an initial value problem and the initial conditions are varied until the boundary conditions at the extreme ends are satisfied. For example, consider the BVP

        \[ \frac{\mathrm{d}^2y}{\mathrm{d}x^2}=f(x) \]

    with the boundary conditions:

        \[ y(0)=y_0\qquad y(L)=y_L \]

    Then, the system can be converted into an initial value problem with the following boundary conditions:

        \[ y(0)=y_0\qquad y'(0)=s \]

    After a solution is obtained, the value of y(L) is compared to y_L and s is varied and the solution is repeated until y(L)=y_L.

  2. Finite difference method. The finite difference method relies on approximating the BVP with difference equations by replacing the derivatives of the function with their finite difference approximations.
  3. Weak formulation methods. The weak formulation methods rely on converting the problem into a different “integral” form and finding the best solution that would satisfy the new form. Examples of such methods are the finite element analysis method or the weighted residuals methods.

The finite difference method is the only method that we will cover in these notes.

Finite Difference Method

Consider a BVP with the dependent variable y and the independent variable x with x\in[a,b]. The finite difference method in a BVP relies on discretizing the domain [a,b] into n intervals with a=x_0<x_1<x_2<\cdots<x_{n-1}<x_n=b with a constant step size h=x_{i+1}-x_i. The values of the dependent variable at the interval junctions are denoted by: y_0,y_1,y_2,\cdots,y_{n-1},y_n. The finite difference method seeks to find a numerical solution to each y_i by solving difference equations formed by replacing the derivatives in the BVP with a finite difference approximation scheme (usually the centred finite difference scheme). Dirichlet boundary conditions can be directly imposed on the known values of y_i while Neumann boundary conditions can be imposed as additional difference equations using an appropriate finite difference scheme. The following examples show how the method can be implemented.

Example 1

In this example, we first derive the heat equation and then attempt to solve it using the finite difference method.
Heat Equation Derivation: The differential equation describing thermal energy within objects (Steady State Heat Equation) is one example of a BVP that can be solved using the finite difference method. In this example, the steady state heat equation is first derived for a one-dimensional rod of length 10m connecting two heat sources of constant temperature as shown in the figure such that the temperature of the rod at the end x=0 is equal to T(0)=40^{\circ} C while the temperature at the other end x=10 is equal to T(10)=200^{\circ} C. The ambient temperature around the rod is equal to T_{out}=20^{\circ} C. The steady state heat equation can be derived by balancing the amount of heat energy per unit time entering and exiting through a control length of the rod. The amount of heat energy transferred to the ambient environment per unit time and per unit length is equal to

    \[ q_{\mbox{ambient}}=B(T(x)-T_{out}) \]

With B=1 units for this example. The amount of heat energy transferred through the rod per unit time is equal to:

    \[ q_{in}=k\frac{\mathrm{d}T}{\mathrm{d}x} \]

With k=-10 units in this example. The negative sign is essential as heat energy is transferred from areas of high temperature to areas of lower temperature. The steady state differential equation can be derived by equating the amount of energy entering to the amount of energy exiting per unit time considering a control length equal to \Delta x as shown in the figure below:

    \[ q_{in}=q_{\mbox{ambient}}+q_{out}= q_{\mbox{ambient}}+q_{in}+\frac{\mathrm{d}q}{\mathrm{d}x}\Delta x \]

Where a Taylor series with only two terms was used to approximate q_{out}. Therefore:

    \[ \frac{\mathrm{d}^2T}{\mathrm{d}x^2}=\frac{1}{10}\left(T(x)-20 \right) \]

The solution to this equation provides the temperature distribution (with T being the dependent variable) within the rod (with x as the independent variable).

temperature
Statement of Problem: Consider the BVP developed above with the Dirichlet boundary conditions T(0)=40 and T(10)=200.

  1. Find the exact solution of the differential equation.
  2. Use the finite difference method with n=4 intervals to find a numerical solution to the BVP.

Repeat the problem considering the two boundary conditions: the Dirichlet boundary condition T(0)=40 and the Neumann boundary condition

    \[ \frac{\mathrm{d}T}{\mathrm{d}x}\bigg|_{x=10}=-3 \]

Solution

Dirichlet Boundary Conditions
The exact solution can be directly obtained in Mathematica using the DSolve function. Note that an exact solution is not always available.

    \[ T(x)= -4.80593 \sinh (0.316228 x)+20. \cosh (0.316228 x)+20. \]

The numerical solution is obtained by first discretizing the domain into n=4 intervals with h=\frac{10}{4}=2.5 and x_0=0, x_1=2.5, x_2=5.0, x_3=7.5, x_4=10.0. The values of T at these points are denoted by T_0, T_1, T_2, T_3, T_4. The boundary conditions are given as Dirichlet type with T_0=40 and T_4=200. The remaining 3 unknowns can be obtained by utilizing the basic formulas of the centred finite difference scheme to replace the derivatives in the differential equation. These can be utilized to write 4 equations of the form:

    \[ \frac{T_{i+1}-2T_i+T_{i-1}}{h^2} =\frac{1}{10}\left(T_i-20\right) \]

This results in the following three equations (after the substitution T_0=40 and T_4=200):

    \[ \begin{split} \frac{T_2-2T_1+40}{6.25} &=\frac{1}{10}\left(T_1-20\right)\\ \frac{T_3-2T_2+T_1}{6.25} &=\frac{1}{10}\left(T_2-20\right)\\ \frac{200-2T_3+T_2}{6.25} &=\frac{1}{10}\left(T_3-20\right) \end{split} \]

Rearranging yields the following linear system of equations:

    \[ \left(\begin{matrix}-0.42&0.16&0\\0.16&-0.42&0.16\\0&0.16&-0.42 \end{matrix}\right)\left(\begin{array}{c}T_1\\T_2\\T_3\end{array}\right)=\left(\begin{array}{c}-8.4\\-2\\-34\end{array}\right) \]

Solving the above system yields:

    \[ \left(\begin{array}{c}T_1\\T_2\\T_3\end{array}\right)=\left(\begin{array}{c}43.1979\\60.8946\\104.15\end{array}\right) \]

Combined with T_0=40 and T_4=200, a list plot can be drawn to show the obtained numerical solution. The following figure shows the exact solution (blue curve) with the numerical solution (black dots). The figure shows that the numerical solution is almost identical to the exact solution.

BVP11
View Mathematica Code

ClearAll["Global`*"]
hp = 1/10;
n = 4;
L = 10 - 0;
h = L/n;
Ts = Table[Symbol["T" <> ToString[i]], {i, 0, n}]
T0 = 40;
Ts[[n + 1]] = 200;
Ts // MatrixForm
Equations =   Table[(Ts[[i + 1]] - 2 Ts[[i]] + Ts[[i - 1]])/h^2 ==     hp*(Ts[[i]] - 20), {i, 2, n}];
Equations // MatrixForm
Equations = Expand[Equations];
Equations // MatrixForm
N[Equations] // MatrixForm

unknowns = Table[Symbol["T" <> ToString[i]], {i, 1, n - 1}]
a = Solve[Equations, unknowns]
Ts = Ts /. a[[1]]
SolList = Table[{0 + (i)*h, Ts[[i + 1]]}, {i, 0, n}];
SolList // MatrixForm

a = DSolve[{T''[x] == hp*(T[x] - 20), T[0] == 40, T[10] == 200}, T, x]
T = T[x] /. a[[1]];
FullSimplify[N[T]]
Plot[T, {x, 0, 10},  AxesLabel -> {"distance x (m)", "Temperature T (Celsius)"},  AxesOrigin -> {0, 0}, Epilog -> {PointSize[Large], Point[SolList]}]

In the following tool, vary the number of intervals n to show the effect on the numerical solution:

Dirichlet and Neumann Boundary Conditions
The exact solution can be directly obtained in Mathematica using the DSolve function. Note that an exact solution is not always available.

    \[ T(x)= -20.7302 \sinh (0.316228 x)+20. \cosh (0.316228 x)+20. \]

The numerical solution is obtained by first discretizing the domain into n=4 intervals with h=\frac{10}{4}=2.5 and x_0=0, x_1=2.5, x_2=5.0, x_3=7.5, x_4=10.0. The values of T at these points are denoted by T_0, T_1, T_2, T_3, T_4. The Dirichlet boundary condition can be used to set T_0=40. The remaining 4 unknowns can be obtained using four equations. Three equations are composed by utilizing the basic formulas of the centred finite difference scheme to replace the derivatives in the differential equation at the three intermediate points. The fourth equation is composed by utilizing the backward finite difference scheme at x_4 and equate it to the Neumann boundary condition. The three difference equations at the intermediate points have the form:

    \[ \frac{T_{i+1}-2T_i+T_{i-1}}{h^2} =\frac{1}{10}\left(T_i-20\right) \]

This results in the following three equations (after the substitution T_0=40):

    \[ \begin{split} \frac{T_2-2T_1+40}{6.25} &=\frac{1}{10}\left(T_1-20\right)\\ \frac{T_3-2T_2+T_1}{6.25} &=\frac{1}{10}\left(T_2-20\right)\\ \frac{T_4-2T_3+T_2}{6.25} &=\frac{1}{10}\left(T_3-20\right) \end{split} \]

The fourth equation can be obtained using the backward finite difference at point x_4:

    \[ \frac{T_4-T_3}{2.5}=-3 \]

Combining the four equations and rearranging yields the following linear system of equations:

    \[ \left(\begin{matrix}-0.42&0.16&0&0\\0.16&-0.42&0.16&0\\0&0.16&-0.42&0.16\\0&0&-0.4&0.4 \end{matrix}\right)\left(\begin{array}{c}T_1\\T_2\\T_3\\T_4\end{array}\right)=\left(\begin{array}{c}-8.4\\-2\\-2\\-3\end{array}\right) \]

Solving the above system yields:

    \[ \left(\begin{array}{c}T_1\\T_2\\T_3\\T_4\end{array}\right)=\left(\begin{array}{c}28.3216\\21.8443\\16.5195\\9.01954\end{array}\right) \]

Combined with T_0=40, a list plot can be drawn to show the obtained numerical solution. The following figure shows the exact solution (blue curve) with the numerical solution (black dots). The figure shows that the numerical solution of the same differential equation with Neumann boundary conditions provides less accurate results than the solution when Dirichlet boundary conditions were used.

BVP12
View Mathematica Code

ClearAll["Global`*"]
hp = 1/10;
n = 4;
L = 10 - 0;
h = L/n;
Ts = Table[Symbol["T" <> ToString[i]], {i, 0, n}]
T0 = 40;
Ts // MatrixForm
Equations =   Table[(Ts[[i + 1]] - 2 Ts[[i]] + Ts[[i - 1]])/h^2 ==     hp*(Ts[[i]] - 20), {i, 2, n}];
NeumannEquation = ((Ts[[n + 1]] - Ts[[n]])/h == -3);
Equations = Append[Equations, NeumannEquation];
Equations // MatrixForm
Equations = Expand[Equations];
Equations // MatrixForm
N[Equations] // MatrixForm

unknowns = Table[Symbol["T" <> ToString[i]], {i, 1, n}]
a = Solve[Equations, unknowns];
Ts = N[Ts /. a[[1]]]
SolList = Table[{0 + (i)*h, Ts[[i + 1]]}, {i, 0, n}];
SolList // MatrixForm

a = DSolve[{T''[x] == hp*(T[x] - 20), T[0] == 40, T'[10] == -3}, T, x]
T = T[x] /. a[[1]];
FullSimplify[N[T]]
Plot[T, {x, 0, 10},  AxesLabel -> {"distance x (m)", "Temperature T (Celsius)"},  AxesOrigin -> {0, 0}, Epilog -> {PointSize[Large], Point[SolList]}]

The numerical solution approaches the exact solution as the number of intervals is increased. Vary n in the following tool and observe how it affects the numerical solution.

 

Example 2

Consider the following second-order BVP:

    \[ u''(x)+u'(x)+u(x)=1 \]

with the boundary conditions:

    \[ u(0)=1.5\qquad u(3)=2.5 \]

Use the finite difference method to find a numerical solution to the above equation by dividing the domain x\in[0,3] into n=5 intervals. Compare with the exact solution.

Solution

The exact solution can be directly obtained in Mathematica using the DSolve function. Note that an exact solution is not always available.

    \[ u=\frac{1}{2} e^{-\frac{x}{2}} \left(2 e^{\frac{x}{2}}+\cos \left(\frac{\sqrt{3} x}{2}\right)-\left(\cos \left(\frac{3 \sqrt{3}}{2}\right)-3 e^{3/2}\right) \csc \left(\frac{3 \sqrt{3}}{2}\right) \sin \left(\frac{\sqrt{3} x}{2}\right)\right) \]

The numerical solution is obtained by first discretizing the domain into n=5 intervals with h=\frac{3}{5}=0.6 and x_0=0, x_1=0.6, x_2=1.2, x_3=1.8, x_4=2.4, x_5=3.0. The values of u at these points are denoted by u_0, u_1, u_2, u_3, u_4, u_5. The boundary conditions are given as Dirichlet type with u_0=1.5 and u_5=2.5. The remaining 4 unknowns can be obtained by utilizing the basic formulas of the centred finite difference scheme to replace the derivatives in the differential equation. These can be utilized to write 4 equations of the form:

    \[ \frac{u_{i+1}-2u_i+u_{i-1}}{h^2}+\frac{u_{i+1}-u_{i-1}}{2h}+u_i=1 \]

This results in the following four equations (after the substitution u_0=1.5 and u_5=2.5):

    \[ \begin{split} \frac{u_{2}-2u_1+1.5}{0.36}+\frac{u_{2}-1.5}{1.2}+u_1&=1\\ \frac{u_{3}-2u_2+u_1}{0.36}+\frac{u_{3}-u_1}{1.2}+u_2&=1\\ \frac{u_{4}-2u_3+u_2}{0.36}+\frac{u_{4}-u_2}{1.2}+u_3&=1\\ \frac{2.5-2u_4+u_3}{0.36}+\frac{2.5-u_3}{1.2}+u_4&=1 \end{split} \]

Rearranging yields the following linear system of equations:

    \[ \left(\begin{matrix}-4.556&3.611&0&0\\1.944&-4.556&3.611&0\\0&1.944&-4.556&3.611\\0&0&1.944&-4.556\end{matrix}\right)\left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right)=\left(\begin{array}{c}-1.91667\\1\\1\\-8.02778\end{array}\right) \]

Solving the above system yields:

    \[ \left(\begin{array}{c}u_1\\u_2\\u_3\\u_4\end{array}\right)=\left(\begin{array}{c}7.649\\9.1183\\7.6615\\5.0323\end{array}\right) \]

Combined with u_0=1.5 and u_5=2.5, a list plot can be drawn to show the obtained numerical solution. The following figure shows the exact solution (blue curve) with the numerical solution (black dots). The maximum error in the numerical solution is obtained at x_2=1.2. At this location, the exact solution is given by u(1.2)=7.6803 while the numerical solution is given by u_2=9.11829 with a relative error of -19%:

    \[ E_r=\frac{7.68-9.12}{7.68}=-0.19 \]

BVP2

View Mathematica Code
ClearAll["Global`*"]
n = 5;
L = 3 - 0;
h = L/n;
us = Table[Symbol["u" <> ToString[i]], {i, 0, n}]
u0 = 1.5;
us[[n + 1]] = 2.5;
us // MatrixForm
Equations = Table[(us[[i + 1]] - 2 us[[i]] + us[[i - 1]])/h^2 + (us[[i + 1]] - us[[i - 1]])/2/h + us[[i]] == 1, {i, 2, n}];
Equations // MatrixForm
Equations = Expand[Equations];
Equations // MatrixForm
N[Equations] // MatrixForm

unknowns = Table[Symbol["u" <> ToString[i]], {i, 1, n - 1}]
a = Solve[Equations, unknowns]
us = us /. a[[1]]
SolList = Table[{0 + (i)*h, us[[i + 1]]}, {i, 0, n}];
SolList // MatrixForm;

Clear[u]
(*Exact Solution*)
b = DSolve[{u''[x] + u'[x] + u[x] == 1, u[0] == 15/10, u[3] == 25/10},   u, x]
u = FullSimplify[u[x] /. b[[1]], Reals]
Plot[u, {x, 0, 3}, PlotLegends -> {"Exact u"},  AxesLabel -> {"x", "u"},  Epilog -> {PointSize[Large], Point[SolList]},  PlotRange -> {{-0.1, 3.1}, {0, 10}}]
Print["u at x=1.2 is ", uexact = u /. x -> 1.2]
Print["u_2 is ", u2 = u2 /. a[[1]]]
Print["Relative Error is ", (uexact - u2)/uexact]

In a BVP, the accuracy of the solution depends on the accuracy of the finite difference scheme used. As shown in the following tool, by decreasing the value of h, the numerical solution gets closer to the exact solution.

 

Example 3

The differential equation describing the displacement of a hanging chain as a function of the weight per unit length w of the chain and the desired horizontal tension in the chain T_0 is given by:

    \[ \frac{\mathrm{d}^2y}{\mathrm{d}x^2}=\frac{w}{T_0}\sqrt{1+\left(\frac{\mathrm{d}y}{\mathrm{d}x}\right)^2} \]

Where y is the dependent variable being the chain curve and x is the independent variable being the horizontal distance. The solution to this equation, namely the function y(x) is called the catenary or the curve that describes the hanging chain. See this page for a detailed derivation and some of the properties of the curve.
If the boundary conditions are given by y(0)=0 and y(4)=0 and if w=1 and T_0=1 units, then use the finite difference method to find a numerical solution by dividing the domain x\in[0,4] into n=5 intervals. Compare with the exact solution.

Solution

Notice that this equation is a nonlinear second-order differential equation. Its exact solution describing the catenary of the chain can be obtained using the built-in Integrate function in Mathematica and has the form:

    \[ y= \cosh (2-x)-\cosh (2)=\frac{1}{2}\left(e^{2-x}+e^{x-2}-e^{2}-e^{-2}\right) \]

The numerical solution is obtained by first discretizing the domain into n=5 intervals with h=\frac{4}{5}=0.8 and x_0=0, x_1=0.8, x_2=1.6, x_3=2.4, x_4=3.2, x_5=4.0. The values of y at these points are denoted by y_0, y_1, y_2, y_3, y_4, y_5. The boundary conditions are given as Dirichlet type with y_0=0 and y_5=0. The remaining 4 unknowns can be obtained by utilizing the basic formulas of the centred finite difference scheme to replace the derivatives in the differential equation. These can be utilized to write 4 nonlinear equations of the form:

    \[ \frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=\sqrt{1+\left(\frac{y_{i+1}-y_{i-1}}{2h}\right)^2} \]

This results in the following four nonlinear equations (after the substitution y_0=0 and y_5=0):

    \[ \begin{split} \frac{y_2-2y_1+0}{0.64}&=\sqrt{1+\left(\frac{y_2-0}{1.6}\right)^2}\\ \frac{y_3-2y_2+y_1}{0.64}&=\sqrt{1+\left(\frac{y_3-y_1}{1.6}\right)^2}\\ \frac{y_4-2y_3+y_2}{0.64}&=\sqrt{1+\left(\frac{y_4-y_2}{1.6}\right)^2}\\ \frac{0-2y_4+y_3}{0.64}&=\sqrt{1+\left(\frac{0-y_3}{1.6}\right)^2} \end{split} \]

Notice that the resulting equations are nonlinear because the differential equation itself was nonlinear. The above system can be solved using the Newton-Raphson method with initial guesses of -0.1 for all of the unknowns which yields the following solution:

    \[ \left(\begin{array}{c}y_1\\y_2\\y_3\\y_4\end{array}\right)=\left(\begin{array}{c}-1.92863\\-2.62693\\-2.62693\\-1.92863\end{array}\right) \]

Combined with y_0=0 and y_5=0, a list plot can be drawn to show the obtained numerical solution. The following figure shows the exact solution (blue curve) with the numerical solution (black dots). The two solutions are very close to each other.
BVP3

View Mathematica Code
ClearAll["Global`*"]
w = 1.;
T = 1;
a = DSolve[{y''[x] == w/T*Sqrt[1 + (y'[x])^2], y[0] == 0, y[4] == 0},y, x]
y = FullSimplify[y[x] /. a[[1]]]
Plot[y, {x, 0, 4}]
L = 4 - 0
n = 5;
h = L/n;
ys = Table[Symbol["y" <> ToString[i]], {i, 0, n}];
y0 = 0;
ys[[n + 1]] = 0;
ys // MatrixForm
Equations =   Table[(ys[[i + 1]] - 2 ys[[i]] + ys[[i - 1]])/h^2 ==     w/T*Sqrt[1 + ((ys[[i + 1]] - ys[[i - 1]])/2/h)^2], {i, 2, n}];
Equations // MatrixForm
unknowns = Table[Symbol["y" <> ToString[i]], {i, 1, n - 1}];
Initvalues = Table[{unknowns[[i]], -0.1}, {i, 1, n - 1}];
a = FindRoot[Equations, Initvalues]
ys = ys /. a
SolList = Table[{0 + (i)*h, ys[[i + 1]]}, {i, 0, n}];
Plot[y, {x, -0.1, 4.1}, AxesLabel -> {"Cable displacement y", "x"}, PlotLegends -> {"Exact Solution"},  Epilog -> {PointSize[Large], Point[SolList]}]

As shown in the following tool, by decreasing the value of h, the numerical solution gets closer to the exact solution.

 

Example 4

Consider the following BVP

    \[ 7\frac{\mathrm{d}^2y}{\mathrm{d}x^2}-2\frac{\mathrm{d}y}{\mathrm{d}x}-y+x=0 \]

With the boundary conditions y(0)=5 and y(20)=8. Assuming h=2, use the finite difference method to find a numerical solution to the given BVP. Compare with the exact solution.

Solution

The exact solution can be directly obtained in Mathematica using the DSolve function.

    \[ y= x+7.00018 e^{-0.261204 x}-0.000178305 e^{0.546918 x}-2. \]

The numerical solution is obtained by first discretizing the domain into n=10 intervals with h=\frac{20}{10}=2 and x_0=0, x_1=2, x_2=4, x_3=6, x_4=8, x_5=10, x_6=12, x_7=14, x_8=16, x_9=18, x_{10}=20. The values of y at these points are denoted by y_0, y_1, y_2, y_3, y_4, y_5, y_6, y_7, y_8, y_9, y_{10}. The boundary conditions are given as Dirichlet type with y_0=5 and y_{10}=8. The remaining 9 unknowns can be obtained by utilizing the basic formulas of the centred finite difference scheme to replace the derivatives in the differential equation. These can be utilized to write 9 equations of the form:

    \[ 7\left(\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}\right)-2\left(\frac{y_{i+1}-y_{i-1}}{2h}\right)-y_i+x_i=0 \]

This results in the following four equations (after the substitution y_0=5 and y_{10}=8):

    \[ \begin{split} 7\left(\frac{y_2-2y_1+5}{4}\right)-\frac{y_2-5}{2}-y_1+x_1&=0\\ 7\left(\frac{y_3-2y_2+y_1}{4}\right)-\frac{y_3-y_1}{2}-y_2+x_2&=0\\ 7\left(\frac{y_4-2y_3+y_2}{4}\right)-\frac{y_4-y_2}{2}-y_3+x_3&=0\\ 7\left(\frac{y_5-2y_4+y_3}{4}\right)-\frac{y_5-y_3}{2}-y_4+x_4&=0\\ 7\left(\frac{y_6-2y_5+y_4}{4}\right)-\frac{y_6-y_4}{2}-y_5+x_5&=0\\ 7\left(\frac{y_7-2y_6+y_5}{4}\right)-\frac{y_7-y_5}{2}-y_6+x_6&=0\\ 7\left(\frac{y_8-2y_7+y_6}{4}\right)-\frac{y_8-y_6}{2}-y_7+x_7&=0\\ 7\left(\frac{y_9-2y_8+y_7}{4}\right)-\frac{y_9-y_7}{2}-y_8+x_8&=0\\ 7\left(\frac{8-2y_9+y_8}{4}\right)-\frac{8-y_8}{2}-y_9+x_9&=0 \end{split} \]

Rearranging yields the following linear system of equations:

    \[ \left(\begin{matrix}-4.5&1.25&0&0&0&0&0&0&0\\2.25&-4.5&1.25&0&0&0&0&0&0\\0&2.25&-4.5&1.25&0&0&0&0&0\\0&0&2.25&-4.5&1.25&0&0&0&0\\0&0&0&2.25&-4.5&1.25&0&0&0\\0&0&0&0&2.25&-4.5&1.25&0&0\\0&0&0&0&0&2.25&-4.5&1.25&0\\0&0&0&0&0&0&2.25&-4.5&1.25\\0&0&0&0&0&0&0&2.25&-4.5\end{matrix}\right)\left(\begin{array}{c}y_1\\y_2\\y_3\\y_4\\y_5\\y_6\\y_7\\y_8\\y_9\end{array}\right)=\left(\begin{array}{c}-13.25\\-4\\-6\\-8\\-10\\-12\\-14\\-16\\-28\end{array}\right) \]

Solving the above system yields:

    \[ \left(\begin{array}{c}y_1\\y_2\\y_3\\y_4\\y_5\\y_6\\y_7\\y_8\\y_9\end{array}\right)=\left(\begin{array}{c}4.19959\\4.51853\\5.50744\\6.89345\\8.50301\\10.2026\\11.824\\13.0018\\12.7231\end{array}\right) \]

Combined with y_0=5 and y_{10}=8, a list plot can be drawn to show the obtained numerical solution. The following figure shows the exact solution (blue curve) with the numerical solution (black dots). The numerical solution exactly overlaps the exact solution!

BVPlast
View Mathematica Code

ClearAll["Global`*"]
a = DSolve[{7 y''[x] - 2 y'[x] - y[x] + x == 0, y[0] == 5, y[20] == 8}, y, x]
y = FullSimplify[y[x] /. a[[1]]]
y = FullSimplify[N[y]]
Plot[y, {x, 0, 20}]
L = 20 - 0
n = 10;
h = L/n;
ys = Table[Symbol["y" <> ToString[i]], {i, 0, n}];
y0 = 5;
ys[[n + 1]] = 8;
ys // MatrixForm
Equations =   Table[7 (ys[[i + 1]] - 2 ys[[i]] + ys[[i - 1]])/h^2 -      2 (ys[[i + 1]] - ys[[i - 1]])/2/h - ys[[i]] + 0 + (i - 1)*h ==     0, {i, 2, n}];
Equations // MatrixForm
Expand[N[Equations]] // MatrixForm
unknowns = Table[Symbol["y" <> ToString[i]], {i, 1, n - 1}];
a = Solve[Equations, unknowns];
ys = N[ys /. a[[1]]]
SolList = Table[{0 + (i)*h, ys[[i + 1]]}, {i, 0, n}];
Plot[y, {x, -0.1, 20.1}, AxesLabel -> {"x", "y"},  PlotLegends -> {"Exact Solution"},  Epilog -> {PointSize[Large], Point[SolList]}]

Problems

  1. The current code given in the explicit Euler method assumes that h is chosen such that \frac{t_{max}-t_0}{h} is a natural number. Modify the code so that if this is not the case, the solution for the dependent variable at t_{max} is available.
  2. Consider the following IVP:

        \[ \frac{\mathrm{d}x}{\mathrm{d}t}=xt+4x \]

    with the initial condition x(0)=1.

    • Use the built-in DSolve function to find the exact solution and plot x for t\in[0,1].
    • Use the explicit Euler method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the implicit Euler method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use Heun’s method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the midpoint method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
    • Use the classical Runge-Kutta method with h=0.1 to find the values of x_i. Plot the exact solution overlapping the obtained numerical solution.
  3. Consider a mass-spring system with k=m=1 units, g=10 units with the initial conditions x(0)=x'(0)=0.
    • Use the built-in DSolve function to find the exact solution and plot the position x and velocity x' for t\in[0,2\pi].
    • Use the explicit Euler method with h=0.2 to find the values of x_i and x'_i. Plot the exact solutions overlapping the obtained numerical solutions.
    • Use the classical Runge-Kutta method with h=0.2 to find the values of x_i and x'_i. Plot the exact solution overlapping the obtained numerical solution.
  4. Consider the following IVP:

        \[ y''(t)+9y=0 \]

    with the initial conditions y(0)=1 and y'(0)=0. Use the implicit Euler method to find a numerical solution to the given IVP for t\in[0,5] with h=0.1. Compare with the exact solution.

  5. Consider the following IVP:

        \[ y''(t)+2y'(t)+6y=\cos{(2t)} \]

    with the initial conditions y(0)=1 and y'(0)=0. Use the classical 4th-order Runge-Kutta method to find a numerical solution to the given IVP for t\in[0,10] with h=0.1. Compare with the exact solution.

  6. Consider the following nonlinear system of IVPs:

        \[\begin{split} \frac{\mathrm{d}x}{\mathrm{d}t}&=10(y-x)\\ \frac{\mathrm{d}y}{\mathrm{d}t}&=x(28-z)-y\\ \frac{\mathrm{d}z}{\mathrm{d}t}&=xy-\frac{8}{3}z\\ \end{split} \]

    with the initial conditions x(0)=-10, y(0)=-10, and z(0)=30. Use the classical 4th-order Runge-Kutta method to find a numerical solution to the given IVP for t\in[0,20] with h=0.01. Then, using list plots, plot the three curves showing y versus x, z versus x, and y versus z.
    Note: The above system is an example of a Lorenz system of equations which can have chaotic solutions depending on the parameters used and the initial conditions.

  7. Use the finite difference method to find a numerical solution to the following BVP:

        \[ 4y''(x)-2y-x=0 \]

    with the boundary conditions y(0)=0 and y(6)=0. Use n=6 intervals, and compare with the exact solution.

  8. The deflection y as a function of the position x along the length of a simply supported Euler-Bernoulli beam uniformly loaded is given by:

        \[ EI\frac{\mathrm{d}^2y}{\mathrm{d}x^2}=\frac{wLx}{2}-\frac{wx^2}{2} \]

    Where E=200,000 MPa is the Young’s modulus, I=30,000 cm^4 is the moment of inertia, w=15 kN/m is the uniformly distributed load, and L=3 m is the length of the beam. The boundary conditions are given by y(0)=y(L)=0. Use the finite difference method to find a numerical solution to the deflection y using n=5 intervals. Compare with the exact solution.

Leave a Reply

Your email address will not be published. Required fields are marked *