The fixed step-size h midpoint rule for numerical integration is yn+1 = yn + hΦ(tn, yn), with k1(tn, yn) = f(tn, yn), k2(tn, yn) = f(tn + 1
computer science
Description
1. The fixed step-size h midpoint rule for numerical integration is yn+1 = yn + hΦ(tn, yn), with k1(tn, yn) = f(tn, yn), k2(tn, yn) = f(tn + 1 2h, yn + 1 2hk1), Φ(tn, yn) = k2(tn, yn). Write a code for numerical integration of the IVP (2), and verify that it is second-order accurate. Use the vector ∞-norm to measure errors. You will have to decide on an appropriate value for the final time, ideally neither too small nor too large. Document your results to convince the reader that your code is correct.
2. The midpoint rule with adaptive step size is yn+1 = yn + hnΦ(tn, yn; hn), where
k1(tn, yn) = f(tn, yn), k2(tn, yn; hn) = f(tn + 1
2hn, yn + 1
2hnk1), Φ(tn, yn, hn) = k2(tn, yn; hn).
The method allows for dynamical adjustment of the step-size hn as needed. Given yn and hn,
compute a prediction y∗ for yn+1. If the error in y∗ is small enough, then set yn+1 = y∗ and try a
larger time-step hn+1 next time. If the error in y∗ is too large, then reject it and start over with a
smaller hn. The local truncation error in y∗ is estimated by comparing two different advancements
of the solution: one by the forward Euler method, and the other by the mid-point rule.
Introduce automatic error control as follows. First, define
where, for example, yi are the components of the state vector yn ∈ Rd at the current time-step.
The request for relative accuracy in the defintion of δ may cause troubles when some components
yi are close to zero; in this case replace |yi| by |y¯i| = max(|yi|, 0.001), or similar. A step is accepted
if δ ≤ tol, in which case the next step is hn+1 = hn min "
1.5,
#tol/(1.2δ)
$
. At step is rejected if
δ > tol, in which case it must be recomputed with the smaller step hn ← hn min "
0.1,
#tol/(1.2δ)
$
.
Write a program implementing the midpoint rule with automatic step size control that can
be applied to a system of differential equations. Store the results so that they can be processed
afterwards, in order to make a table of the results, and/or curves showing y(t) versus t, or, say for
a two-component system, y2 versus y1. Use your code to study the IVPs (1) and (2).
