Preface

---

This is a tutorial made solely for the purpose of education and it was designed for students taking Applied Math 0330. It is primarily for students who have very little experience or have never used Mathematica before and would like to learn more of the basics for this computer algebra system. As a friendly reminder, don't forget to clear variables in use and/or the kernel.

Finally, the commands in this tutorial are all written in bold black font, while Mathematica output is in normal font. This means that you can copy and paste all commands into Mathematica, change the parameters and run them. You, as the user, are free to use the scripts for your needs to learn the Mathematica program, and have the right to distribute this tutorial and refer to this tutorial as long as this tutorial is accredited appropriately.

Heun's Methods

---

Since for many problems the Euler rule requires a very small step size to produce sufficiently accurate results, many efforts have been devoted to the development of more efficient methods. Our next step in this direction includes Heun's method, which was named after a German mathematician Karl Heun (1859--1929), who made significant contributions to developing the Runge-Kutta method. Actually, the algorithm named after Heun was discovered by Carl Runge.

I. Trapezoid Rule

---

Suppose that we want to solve the initial value problem for the first order differential equation

which is equivalent to the integral equation. To find its numerical solution, we pick up the uniform set of grid points \( x_n = x_0 + h\, n , \quad n=0,1,2,\ldots ; \) where h is a step size chosen fixed for simplicity. As usual, we denote by y_n approximate values of unknown function at mesh points \( y_n \approx \phi (x_n ) , \quad n=0,1,2,\ldots ; \) where φ is the actual solution of the given initial value problem (IVP for short).

The concept of the Trapezoidal Rule in numerical methods is similar to the trapezoidal rule of Riemann sums. The Trapezoid Rule is generally more accurate than the Euler approximations, and it calculates approximations by taking the sum of the function of the current and next term and multiplying it by half the value of h. Therefore the syntax will be as follows:

This is an implicit method: the value y_n+1 appears on both sides of the equation, and to actually calculate it, we have to solve an equation which will usually be nonlinear. One possible method for solving this equation is Newton's method. We can use the Euler rule to get a fairly good estimate for the solution, which can be used as the initial guess of Newton's method.

 

II. Heun's Formula / Improved Euler Method

---

The Improved Euler’s method, also known as the Heun formula or the average slope method, gives a more accurate approximation than the Euler rule and gives an explicit formula for computing y_n+1. The basic idea is to correct some error of the original Euler's method. The syntax of the Improved Euler’s method is similar to that of the trapezoid rule, but the y value of the function in terms of y_n+1 consists of the sum of the y value and the product of h and the function in terms of x_n and y_n.

Improved Euler formula or the average slope method is commonly referred to as the Heun method (although discovered by C. Runge):

Since it is actually the simplest version of predictor-corrector method, the recurrence can be writen as Therefore, one of the versions on the Heun method in Mathematica is as follows
Clear[y]
f[x_, y_] := x^2 + y^2
y[0] = 1;
h:=0.1;
Do[k1 = h f[h n, y[n]];
k2 = h f[h (n + 1), y[n] + k1];
y[n + 1] = y[n] + .5 (k1 + k2),
{n, 0, 9}]
y[10]

The main difference here is that instead of using y[n+1]=y[n]+h*f(x[n],y[n]), we use

i.e. we are averaging the slope. For a build-in version of this, input the following:

Instead of calculating intermediate values k1 and k2 , one can define ynew directly:

The function heun is basically used in the same way as the function euler. Then to solve the above problem, do the following:

Then read the answer off the table.

Another option:

Then implement Heun's rule:

 

Example:

Solve the IVP: \( y'= 1/(3x-2y+1),\quad y(0)=0 \)

First, we find its solution using Mathematica's command DSolve:

To evaluate the exact solution at x=1.0, we type

Then we use Heun's numerical method with step size h=0.1:

Now we demonstrate the use of while command:

Example:

Consider the initial value problem for the integro-differential equation First, we recall y₁ found previously based on Euler rule that we use as a predictor: Now we correct this value based on trapezoid rule: We approximate again the double integral in the right hand-side with the trapezoid rule: Therefore, The general case of Heun's method yields where We break the integral into sum Using again the trapezoid rule, we get The integral in the right hand-side can be iterated: with \( J_1 = \frac{h}{2} \left( y_1 + 50 \right) . \) Now we ask Mathematica to do calculations.

 

Make a computational experiment: in the improved Euler algorithm, using the corrector value as the next prediction, and only after the second iteration use the trapezoid rule.

The ``improved'' predictor-corrector algorithm:

is realized using the following Mathematica script:


Then we apply this code to numerically solve the initial value problem \( y' = (1+x)\,\sqrt{y} , \quad y(0) =1 \) with step size h=0.1: we get approximate value 9.00778 while the regular Heun's method gives 8.99148.

 

III. Modified Euler Method / Midpoint Method

---

The Modified Euler’s method is also called the midpoint approximation. This method reevaluates the slope throughout the approximation. Instead of taking approximations with slopes provided in the function, this method attempts to calculate more accurate approximations by calculating slopes halfway through the line segment. The syntax of the Modified Euler’s method involves the sum of the current y term and the product of h with the function in terms of the sum of the current x and half of h (which defines the x value) and the sum of the current y and the product of the h value and the function in terms of the current x and y values (which defines the y value).

The midpoint method can be shown to have a local error of 2, so it is second-order accurate. The midpoint method is implemented in NDSolve as "ExplicitMidpoint":

Modified Euler formula or explicit midpoint rule or midpoint Euler algorithm:

Therefore, the Mathematica syntax is as follows:

Another option:

Then implement modified Euler method/ midpoint rule:

Let’s take all of the approximations and the exact values to compare the accuracy of the approximation methods:

x values	Exact	Euler	Backwards Euler	Trapezoid	Improved Euler	Modified Euler
0.1	1.049370088	1.050000000	1.057930422	1.0493676	1.049390244	1.049382716
0.2	1.097463112	1.098780488	1.118582822	1.0974587	1.097594738	1.097488615
0.3	1.144258727	1.146316720	1.182399701	1.1442530	1.144322927	1.144297185
0.4	1.189743990	1.192590526	1.249960866	1.1897374	1.189831648	1.189795330
0.5	1.233913263	1.237590400	1.322052861	1.2339064	1.234025039	1.233977276
0.6	1.276767932	1.281311430	1.399792164	1.2767613	1.276904264	1.276844291
0.7	1.318315972	1.323755068	1.484864962	1.3183102	1.318477088	1.318404257
0.8	1.358571394	1.364928769	1.580059507	1.3585670	1.358757326	1.358671110
0.9	1.397553600	1.404845524	1.690720431	1.3975511	1.397764204	1.397664201
1.0	1.435286691	1.443523310	1.830688225	1.4352865	1.435521666	1.435407596
Accuracy	N/A	99.4261%	72.4513%	99.9999%	99.9836%	99.9916%