Problem Set 1

Topics

  • Initial value problems
  • Explicit and implicit time integration
  • Stability regions

P1.1 Forward Euler on a scalar IVP

Consider the initial value problem \[\begin{cases} y' = -2y, \\ y(0) = 1. \end{cases}\]
  1. Find the exact solution.
  2. Write the Forward Euler update formula with time step \(\tau\).
  3. Compute approximations of \(y(1)\) using \(\tau = 0.5\), \(0.25\) and \(0.1\).
  4. Compare the numerical results with the exact solution.
  5. Comment on the observed convergence as \(\tau\) decreases.

Solution

Forward Euler method

Pair-wise ratios of the final errors (highlighted in bold in the tables below).
Time-steps Pairwise ratio
\(\tau = 0.25\) and \(0.125\) \(2.068\)
\(\tau = 0.125\) and \(0.0625\) \(2.040\)

See the accompanying Excel spreadsheet for detailed calculations.

\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.25\) \(0.6065\) \(0.5\) \(0.1065\)
\(2\) \(0.5\) \(0.3679\) \(0.25\) \(0.1179\)
\(3\) \(0.75\) \(0.2231\) \(0.125\) \(0.09813\)
\(4\) \(1\) \(0.1353\) \(0.0625\) \(\mathbf{0.07284}\)
\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.125\) \(0.7788\) \(0.75\) \(0.02880\)
\(2\) \(0.25\) \(0.6065\) \(0.5625\) \(0.04403\)
\(3\) \(0.375\) \(0.4724\) \(0.4219\) \(0.05049\)
\(4\) \(0.5\) \(0.3679\) \(0.3164\) \(0.05147\)
\(5\) \(0.625\) \(0.2865\) \(0.2373\) \(0.04920\)
\(6\) \(0.75\) \(0.2231\) \(0.1780\) \(0.04515\)
\(7\) \(0.875\) \(0.1738\) \(0.1335\) \(0.04029\)
\(8\) \(1\) \(0.1353\) \(0.1001\) \(\mathbf{0.03522}\)
\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.0625\) \(0.8825\) \(0.875\) \(0.007497\)
\(2\) \(0.125\) \(0.7788\) \(0.7656\) \(0.01318\)
\(3\) \(0.1875\) \(0.6873\) \(0.6699\) \(0.01737\)
\(4\) \(0.25\) \(0.6065\) \(0.5862\) \(0.02035\)
\(5\) \(0.3125\) \(0.5353\) \(0.5129\) \(0.02235\)
\(6\) \(0.375\) \(0.4724\) \(0.4488\) \(0.02357\)
\(7\) \(0.4375\) \(0.4169\) \(0.3927\) \(0.02417\)
\(8\) \(0.5\) \(0.3679\) \(0.3436\) \(0.02427\)
\(9\) \(0.5625\) \(0.3247\) \(0.3007\) \(0.02399\)
\(10\) \(0.625\) \(0.2865\) \(0.2631\) \(0.02343\)
\(11\) \(0.6875\) \(0.2528\) \(0.2302\) \(0.02265\)
\(12\) \(0.75\) \(0.2231\) \(0.2014\) \(0.02171\)
\(13\) \(0.8125\) \(0.1969\) \(0.1762\) \(0.02067\)
\(14\) \(0.875\) \(0.1738\) \(0.1542\) \(0.01956\)
\(15\) \(0.9375\) \(0.1534\) \(0.1349\) \(0.01842\)
\(16\) \(1\) \(0.1353\) \(0.1181\) \(\mathbf{0.01727}\)

Midpoint rule

Pair-wise ratios of the final errors (highlighted in bold in the tables below).
Time-steps Pairwise ratio
\(\tau = 0.25\) and \(0.125\) \(4.051\)
\(\tau = 0.125\) and \(0.0625\) \(4.013\)

See the accompanying Excel spreadsheet for detailed calculations.

\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.25\) \(0.6065\) \(0.6\) \(0.006531\)
\(2\) \(0.5\) \(0.3679\) \(0.36\) \(0.007879\)
\(3\) \(0.75\) \(0.2231\) \(0.216\) \(0.007130\)
\(4\) \(1\) \(0.1353\) \(0.1296\) \(\mathbf{0.005735}\)
\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.125\) \(0.7788\) \(0.7778\) \(0.001023\)
\(2\) \(0.25\) \(0.6065\) \(0.6049\) \(0.001592\)
\(3\) \(0.375\) \(0.4724\) \(0.4705\) \(0.001859\)
\(4\) \(0.5\) \(0.3679\) \(0.3660\) \(0.001929\)
\(5\) \(0.625\) \(0.2865\) \(0.2846\) \(0.001877\)
\(6\) \(0.75\) \(0.2231\) \(0.2214\) \(0.001753\)
\(7\) \(0.875\) \(0.1738\) \(0.1722\) \(0.001592\)
\(8\) \(1\) \(0.1353\) \(0.1339\) \(\mathbf{0.001416}\)
\(n\) \(t _ n\) \(q ( t _ n )\) \(q _ n\) \(\epsilon _ n\)
\(0\) \(0\) \(1\) \(1\) \(0\)
\(1\) \(0.0625\) \(0.8825\) \(0.8824\) \(0.0001440\)
\(2\) \(0.125\) \(0.7788\) \(0.7785\) \(0.0002541\)
\(3\) \(0.1875\) \(0.6873\) \(0.6870\) \(0.0003363\)
\(4\) \(0.25\) \(0.6065\) \(0.6061\) \(0.0003957\)
\(5\) \(0.3125\) \(0.5353\) \(0.5348\) \(0.0004364\)
\(6\) \(0.375\) \(0.4724\) \(0.4719\) \(0.0004622\)
\(7\) \(0.4375\) \(0.4169\) \(0.4164\) \(0.0004758\)
\(8\) \(0.5\) \(0.3679\) \(0.3674\) \(0.0004798\)
\(9\) \(0.5625\) \(0.3247\) \(0.3242\) \(0.0004763\)
\(10\) \(0.625\) \(0.2865\) \(0.2860\) \(0.0004670\)
\(11\) \(0.6875\) \(0.2528\) \(0.2524\) \(0.0004533\)
\(12\) \(0.75\) \(0.2231\) \(0.2227\) \(0.0004364\)
\(13\) \(0.8125\) \(0.1969\) \(0.1965\) \(0.0004172\)
\(14\) \(0.875\) \(0.1738\) \(0.1734\) \(0.0003964\)
\(15\) \(0.9375\) \(0.1534\) \(0.1530\) \(0.0003748\)
\(16\) \(1\) \(0.1353\) \(0.1350\) \(\mathrm{0.0003528}\)

P1.2 Stability contrast: Forward vs Backward Euler

Consider the IVP \[\begin{cases} y' = -50 y, \\ y(0) = 1. \end{cases}\]
  1. Write the Forward Euler update formula.
  2. For which values of \(\tau\) does the method remain stable?
  3. Write the Backward Euler update formula.
  4. Show that the method is unconditionally stable.
  5. Compute the numerical solution up to \(t=1\) using \(\tau = 0.1\) for both methods.
  6. Compare the results and comment on stability.

P1.3 Accuracy comparison: Euler vs Midpoint

Consider the IVP \[\begin{cases} y'(t) = \sin(t) - y(t), \\ y(0) = 1. \end{cases}\]
  1. Write the update formulas for:
    • Forward Euler method
    • Midpoint rule
  2. Compute the numerical solution on \([0,2]\) using \(\tau = 0.2\).
  3. Plot the numerical solutions together with the exact solution.
  4. Rank the methods according to accuracy.
  5. Identify which methods are first-order and which are second-order accurate.

P1.4 Mechanical system

A mass–spring–damper system is governed by \[ m \ddot{x} + c \dot{x} + k x = 0, \] with parameters \(m = 1\), \(c = 0.2\), \(k = 1\).

  1. Rewrite the second-order equation as a first-order system.
  2. Apply the following time-integration schemes:
    • Forward Euler
    • Backward Euler
    • Midpoint rule
    with initial conditions \(x(0) = 1\), \(x'(0) = 0\).
  3. Plot displacement and velocity versus time.
  4. Discuss how each numerical method affects the energy decay of the system.

P1.5 Nonlinear IVP

Consider the nonlinear initial value problem \[\begin{cases} y' = y(1 - y), \\ y(0) = 0.1. \end{cases}\]
  1. Apply the forward Euler method to the numerical solution of this problem.
  2. Consider applying the backward Euler method: what type of equation must be solved at each time step?