I - An Embedded Runge-Kutta Methods

Task 1: Write out the formulas for k₁, . . . , k₄, y_n+1 and y^{^}_n+1.

Task 2: Explain how this specific method can be implemented so that only 3 function evaluations per step are required.

Task 3: Write out R(z) in the same way and state the values of the coefficients α₁, α₂ and α₃.

Task 4: Use Maple to make a plot of the polynomials R(z) and R^{^}(z) for z ∈ [-3, 1]. Use one graph to show both polynomials together.

Task 5: Use Maple to find the value z* such that |R(z)| < 1 for all values z ∈ (z*, 0). In the same way, use Maple to find the value z^* such that |R^{^}(z)| < 1 for all values z ∈ (z^*, 0). Hint: You can use the fsolve command.

Task 6: Suppose we use the RK2 method with a constant step size h for the system of ODEs

X' = Ax = XΛX^-1x, where X is an arbitrary 2 × 2 matrix and

Use the results from the previous task to find the upper bound on h required to obtain a stable method. Also find the upper bound for the RK3 method. Which of the two methods gives the most severe restriction on h?

II - Constant Step Size Implementation

Task 1: Implement a procedure that takes as arguments the parameters a, b, c ∈ R and z : R → R, the initial condition t₀, x₀, y₀ ∈ R, the final time t_e ∈ R and the number of steps N ∈ N and calculates an approximation to the solution of the Fitz Hugh-Nagumo model (1.1)-(1.2). The output of your procedure should be five arrays t[0..N], x[0..N], y[0..N], x^{^}[0..N], y^{^}[0..N] containing the grid points t_n and the numerical approximations x_n, y_n, x^{^}n and y^{^}n. It may be a good idea to base your implementation on several smaller procedures.

Task 2: Use the procedure from the previous task to calculate an approximate solution to the Fitz Hugh-Nagumo problem (1.1)-(1.2) for the parameter values a = 0.7, b = 0.8, c = 1, 3, 9 and z(t) = -0.4. Use t₀ = 0, x₀ = 1.2 and y₀ = -0.625 as the values for the initial condition, t_e = 32 and N = 400. For each value for c create three graphs:

• x_n and y_n versus t_n,
• the error estimate Ln =
  

  , versus t_n

• a phase plot: y_n versus x_n.

The error estimate is measured using the norm ||[^x_y]||₂ = √(x² + y²). Use a logarithmic scale for the error estimate L_n (see plots[logplot]).

It may be a good idea to write several procedures to solve this task. For example, you could write a procedure that takes as arguments N the arrays t[0..N], x[0..N], y[0..n], x^{^}[0..N] and y^{^}[0..N] and produces the necessary graphs. This procedure can then be reused in the next section.

Briefly describe the results and how the parameter c influences the solution and the error estimate.

III - Variable Step Size Implementation

Task 1: Implement a procedure that takes as arguments the parameters a, b, c ∈ R and z : R → R, the initial condition t₀, x₀, y₀ ∈ R, the final time t_e ∈ R, the initial step size h0, the error tolerance  and the maximum number of time steps N_max and calculates an approximation to the solution of the Fitz Hugh-Nagumo model (1.1)-(1.2) using the variable step size technique just described. The output of your procedure should be the number of time steps N and five arrays t[0..N], x[0..N], y[0..N], x^{^}[0..N], y^{^}[0..N] containing the grid points t_n and the numerical approximations x_n, y_n, x^{^}n and y^{^}n. Note that as far as input values are concerned the only difference with the constant step size procedure is that instead of the input parameter N we have the parameters h0, ε and N_max. The number of time steps N is an extra output value for the variable step procedure.

Because we do not know beforehand how many steps the procedure will take, we cannot use Maple Array data structures directly. Instead we will use table data structures and convert them to Array data structures at the end.

Task 2: Repeat Task 2 in part II, but use the variable step procedure rather than the constant step procedure. Use the same parameter values as in Task 2, but use h₀ = 10^-3, ε = 10^-5 and N_max = 400 (rather than specifying N as was done for the constant step case). In addition to the graphs of x_n, y_n versus t_n, the error estimate L_n versus t_n and the phase plot with x_n versus y_n, also generate a graph of the step size h_n versus t_n. Use a logarithmic scale for the h_n axis. You can either implement new procedures for this task (possibly using procedures from the previous section) or you can adapt the procedures from the previous section so that they can use both the constant and variable step methods. Briefly describe your observations.

References-

[Fit61] R. FitzHugh. Impulses and physiological states in theoretical models of nerve membrane. Biophysical Journal, 1:445-466, 1961.

[NAY62] J. S. Nagumo, S. Arimoto, and S. Yoshizawa. An active pulse transmission line simulating nerve axon. Proceedings of the Institute of Radio Engineers, 50:2061-2071, 1962.