MATHEMATICA TUTORIAL for the First Course. Part III.

 

---

Return to computing page
Return to computing page for the second course APMA0330
Return to Mathematica tutorial page for the first course
Return to Mathematica tutorial page for the second course
Return to the main page for the course APMA0330
Return to the main page for the course APMA0340

Contents

Part III Numerical Methods and Applications

a) Euler methods
b) Polynomial approximations
c) Runge-Kutta methods
d) Multistep methods
4) Numerov's method

Applications

 

 

3.1. Numerical Solution using DSolve and NDSolve

---

A complete analysis of a differential equation is almost impossible without utilizing computers and corresponding graphical presentations. We are not going to pursue a complete analysis of numerical methods, which usually requires:

an understanding of computational tools;
an understanding of the problem to be solved;
construction of an algorithm which will solve the given mathematical problem.

In this chapter, we concentrate our attention on presenting some simple numerical algorithms for finding solutions of the initial value problems for the first order differential equations in the normal form:

assuming that the given problem has a unique solution.

A general approach to find numerically solutions of the initial value problems and then plot them using Mathematica's capability is as follows:

The above script was used to solve (and then plot) the non-linear differential equstion \( y' = (1-x)/(1+y^5) \) subject to the initial condition \( y(0) =1.\) To see a numerical value at, say x=3.0, we type:

Since its solution is known in implicit form, there is an alternative way to plot the graph:

Example 1.3.1:

Mathematica represents the solution through the special function, called Fresnel function.
?FresnelS
FresnelS[z] gives the Fresnel integral S(z).  >>

This is a plot of the numerical solution. Now we calculate its approximate value at x=6.5 as follows:

This indicates that y(6.5) approximately equal to 0.639918. We can improve accuracy by choosing the option PrecisionGoal

There are two main approaches to find a numerical value for the solution to the initial value problem y' = f(x,y), y( x0 )= y0 when its solution can be obtained using either DSolve (for analytical solution, if it is possible) or NDSolve (for numerical solution) command. The first one is based on definition of the slope function f(x,y) as an expression. The second one (which is recommended) is to define the slope as a function. To see how it works, let us consider an example.

Example 1.3.2: Consider the slope function f(x,y) = 1/(3*x - 2*y + 1). Then we find a numerical value treating f(x,y) as an expression.

Now we use the function DSolve to find explicit solution:

Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >>
Now we define solution as a function of x:

We present two other options to find a numerical value of a solution to the initial value problem using NDSolve.

or

 

3.2. Euler's Methods

---

I. Euler Method

In 1768, Leonhard Euler (St. Petersburg, Russia) introduced a numerical method that is now called the Euler method or the tangent line method for solving numerically the initial value problem:

There are several ways to implement the Euler numerical method for solving initial value problems. We demonstrate its implementation in a series of codes. To start, define the initial point and then the slope function:

Next, define the step size:

Now we define the Euler method itself:
Create the table of approximations using Euler's rule:

Plot with some options:

or
or

Another way is to make a loop.

Explanation. First of all, it is always important to clear all previous assignment to all the variables that we are going to use, so we have: Clear[y]

The basic structure of the loop is:
Do[some expression with n, {n, starting number, end number}]
or with option increment, denoted by k:
Do[some expression with n, {n, starting number, end number, k}]

The function of this Do loop is to repeat the expression, with n taking values from “starting
number” to “ending number,” and therefore repeat the expression (1+ (ending number) - (starting
number))  many times. For our example, we want to iterate 10 steps, so n will go from 0, 1, 2, ...,
until 9, which is 10 steps in total. There is just one technical issue: we must have n as an
integer. Therefore we let y[n] denote the actual value of y(x) when x=n*h. This way everything is an
integer, and we have our nice do loop.
Finally put y[10] , which is actualy y(2), then shift+return, and we have our nice answer.

First, we start with output of our program, which is, perhaps, the most important part of our program design. You'll notice that in the program description, we are told that the "solution that it produces will be returned to the user in the form of a list of points. The user can then do whatever one likes with this output, such as create a graph, or utilize the point estimates for other purposes.

Next, we must decide upon the information that must be passed to the program for it to be able to achieve this. In programming jargon, these values are often referred to as parameters. So what parameters would an Euler program need to know in order to start solving the problem numerically? Quite a few should come to mind:

The slope function f(x, y)
The initial x-value, x₀
The initial y-value, y₀
The final x-value,
The number of subdivisions to be used when chopping the solution interval into individual jumps or the step size.

To code these parameters into the program, we need to decide on actual variable names for each of them. Let's choose variable names as follows:

f, for the slope function f(x, y)
x₀,      for the initial x-value, x₀
y₀,      for the initial y-value, y₀
x_n,      for the final x-value
Steps, for the number of subdivisions
h,        for step size

There are many ways to achieve this goal. However, it is natural to have the new subroutine that looks similar to build-in Mathematica's commands. So we want our program might look something like this:

In order to use Euler's method to generate a numerical solution to an initial value problem of the form:

We have to decide upon what interval, starting at the initial point x0, we desire to find the solution. We chop this interval into small subdivisions of length h, called step size. Then, using the initial condition \( (x_0,y_0) \) as our starting point, we generate the rest of the solution by using the iterative formulas:

to find the coordinates of the points in our numerical solution. The algorithm is terminated when the right end of the desired interval has been reached.

If you would like a built-in Euler’s method, unluckily there is none. However, we can define it ourself. Simply copy the following command line by line:

Now we have our euler function:   euler[f(x,y), {x,x0,x1},{y,y0},steps]
Then this script will solve the differential equation y’=f(x,y), subject to the initial condition y(x0)=y0, and generate all values
between x0 and x1. The number of steps for the Euler’s method is specified with steps.
To solve the same problem as above, we simply need to input:

Another example:

As you can see, the output is really a big table of values, and you can really just read the last one off to get y(1). If the only one final value is needed, then Return command should be replaced with

We can plot the output as follows

Or we can plot it as

The following code, which uses a slightly different programming paradigm, implements Euler's method for a system of differential equations:

And the usage message is:

euler::usage = "euler[F, t0, Y0, b, n] gives the numerical solution to {Y' == F[t, Y], Y[t0] == Y0} over the interval\n
[t0, b] by the n-step Euler's method. The result is in the form of a table of {t, Y} pairs."

Note that this function uses an exact increment h rather than converting it explicitly to numeric form using N. Of course you can readily hange
that or accomplish the corresponding thing by giving an approximate number for one of the endpoints.

Next we plot the points by joining them with a curve:

Another way without writing a subroutine:

The results can also be visualized by connecting the points:

Next, we demonstrate application of Function command:

Here is another steamline of the Euler method (based on Mathematica function).

Now we can plot our results, making sure we only refer to values of x and y that we have defined:

Now we are going to repeat the problem, but using Mathematica lists format instead. The idea is still the same---we set a new x value, compute the corresponding y value, and save the two quantities in two lists. In particular, in a list we are simply storing values; In previously presented code, we were actually storing rules, which are more complicated to store and evaluate. We start out with just one value in each list (the initial conditions). Then we use the Append command to add our next pair of values to the list.

Now we can plot our results:

Another option:

Euler method:

Euler's method is first-order accurate, so that errors occur one order higher starting at powers of \( h^2 . \)
Euler's method is implemented in NDSolve as ExplicitEuler:

Euler's method is first-order accurate, so that errors occur one order
higher starting at powers of \( h^2 . \)
Euler's method is implemented in NDSolve as ExplicitEuler:

 

 

II. Backward Euler Method

---

 

Backward Euler formula:

where h is the step size. Consider the following initial value problem:

Here is the Mathematica code that solve this problem:

Then we plot solution:

NSolve has to spend time to compute all roots to the equation (which can be computationally expensive). FindRoot does a pretty fast search looking for only a single root, so it is quick for complex equations.
If you care about all possible roots, or if you have no clue where the roots of the equation may be, FindRoot is a terrible choice. If you only care about a single root and have a rough idea of where it might be, though, FindRoot will find it quickly.

 

 

III.Trapezoid Rule

---

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 above 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 Mathematica syntax will be as follows:

 

 

IV.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 trapezoid rule and gives an explicit formula for computing y(n+1) in terms of the values of x. 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 or it is commonly referred to as the Heun method (although discovered by C. Runge):

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 1.3.1:

Solve the IVP: \( y'= 1/(3x-2y+1), 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:

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 {\em 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.

 

V. 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%

Notice how each improved approximation approaches the exact value, which results in the modified Euler’s method being the most accurate method of approximating the value of the differential equation y’=1/(2*x-3*y+5). Note that the modified Euler’s method isn’t always the most accurate approximation method, as the improved Euler’s method can be more accurate than the modified Euler’s method depending on the differential equation in question.

 

 

1.3.3. Polynomial Approximations

---

II. Polynomial Approximations

Since the first order Taylor series approximation is identical with Euler’s method, we start with the second order one. Now let us
take example where y’=1/(2x-3y+5), y(0)=1 and let us try to find out y(1).

Approximations using Taylor series expansions is a much easier task to perform using Mathematica because there are many iterations of the same calculations performed in order to achieve the desired answer.

Approximations with using Taylor series expansions in the first-order is actually the Euler algorithm: y_n+1=y_n+hf(x_n,y_n)=y_n+h/(1-2x_n+y_n), n=0,1,2,…

They are similar in that they both take in the current term and add the approximate deviation to the next term in order to yield the value of the next term. If you look at the Mathematica file, the approximation of the sequence is the last y term of the sequence, which in this case is y₁₀=1.443523310.

For each step in the sequence, the expansion requires 3 addition, 1 subtraction, 3 multiplication, and 1 division operations. Since there are 8 steps to this sequence, in order to obtain the desired approximation, the code performs 80 operations, which would be hectic work to do by hand.

 

 

II. Second Order Polynomial Approximation

---

The second order Taylor approximation is just adding the second order differential deviation to the next term in the same equation used for the first order Taylor expansion. The approximation is y10=1.435290196. This algorithm requires at least 22 operations per step, which means that the entire sequence requires 220 steps. However, this is not significant since Mathematica does all the work.

Since the first order Taylor series approximation is identical with Euler’s method,
we start with the second order one. Now let us
take example where y’=1/(2x-3y+5), y(0)=1 and let us try to find out y(1).

Using Mathematica command NDSolve, we check the obtained value:

Another example.

III. Third Order Polynomial Approximation

---

The third order Taylor approximation is adding a third order differential deviation to the equation for the 2nd order expansion. The approximation for this is y₁₀=1.435283295.
Here is the syntax for the third order Taylor approximation (which takes too long to execute):

The total number of numerical operations here is 420. It is obvious at this point why using mathematical programs such as Mathematica is a desired approach for such problems.

Another example:

IV. Regular Plotting vs Log-Plotting

---

Plotting the absolute value of the errors for Taylor approximations, such as the ones from the previous example, is fairly easy to do in normal plotting.

The normal plotting procedure involves taking all of the information regarding the Taylor approximations and putting them into one file, then calculate the exact solution, then find the difference between each approximation and the exact solution, and plotting the solution.

Looking at the code, the syntax is pretty much self-explanatory. The only part that may confuse students is the syntax for the plot command. The plot command need not be in the same syntax as the file has it. If you wish to, you can alter the syntax to form a curve instead of dots.

Log-plotting is not much different than plotting normally. In other programs such as Mathematica, there are a couple of fixes to make, but in Maple only one line of code is required, as seen in the file. Make sure that before you use the logplot command to open the plots package.

 

1.3.4. Runge-Kutta Methods

---

Both previously discussed rules (improved Euler and modified Euler) are particular cases of the family of numerical methods known as the Runge-Kutta methods.

The idea of Runge--Kutta methods is to take successive (weighted) Euler steps to approximate a Taylor series. In this way function evaluations (and not derivatives) are used. For example, consider the one-step formulation of the midpoint method.

Extending the approach in (1), repeated function evaluation can be used to obtain higher-order methods. Denote the Runge--Kutta method for the approximate solution to an initial value problem at \( t_{n+1} = t_n +h \) by

It is generally assumed that the row-sum conditions hold:

These conditions effectively determine the points in time at which the function is sampled and are a particularly useful device in the derivation of high-order Runge--Kutta methods.

The coefficients of the method are free parameters that are chosen to satisfy a Taylor series expansion through some order in the time step \( h . \) In practice other conditions such as stability can also constrain the coefficients.

Explicit Runge--Kutta methods are a special case where the matrix A is strictly lower triangular:

It has become customary to denote the method coefficients \( {\bf c} = [c_i ]^T ,\ {\bf b} = [b_i]^T \) , and \( {\bf A} = [a_{i,j}] \) using a Butcher table, which has the following form for explicit Runge--Kutta methods.

The row-sum conditions can be visualized as summing across the rows of the table.

Notice that a consequence of explicitness is \( c_1 =0 ,\) so that the function is
sampled at the beginning of the current integration step.

Example.The Butcher table for the explicit midpoint method is given by:

See
http://reference.wolfram.com/language/tutorial/NDSolveExplicitRungeKutta.html

 

 

The Euler’s method is sometimes called the first order Runge-Kutta Method, and the Heun’s
method the second order one. In this section, we consider a family of Runge-Kutta methods. To make our exposition uniform, we choose the step size h=0.1, and test each algorithm on the following initial value problem:

\[ y' = x^2 - y^2 , \qquad y(0)=0. \]

Correspondingly, we define the slope function f[x_,y_]=x*x-y*y;

I. Second Order Runge-Kutta Approximations

---

Second Order Optimal Runge-Kutta Approximation when the step size is specified: \begin{equation} \label{E35.opt} y_{n+1} =y_n + \frac{h}{4}\,f(x_n , y_n ) + \frac{3h}{4} \,f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation}

The same algorithm when the number of steps is specified:

If only one final value is needed, then Return command should be replaced with

 

Second Order Ralson Approximation:

\begin{equation*} %\label{E35.ralson} y_{n+1} =y_n + \frac{h}{3}\,f(x_n , y_n ) + \frac{2h}{3} \,f\left( x_n + \frac{3h}{4} , y_n + \frac{3h}{4} \,f(x_n , y_n ) \right) , \quad n=0,1,2,\ldots . \end{equation*}
The same algorithm with fixed step size (which is denoted by h) is specified:

Implicit Runge-Kutta of order 2:

\begin{equation} \label{E36.implicit} y_{n+1} = y_n + \frac{h}{2}\left( k_1 + k_2 \right) , \end{equation} where \begin{align*} k_1 &= f\left( x_n + \frac{\sqrt{3} -1}{2\sqrt{3}}\,h , \ y_n + \frac{h}{4}\,k_1 + \frac{\sqrt{3} -2}{4\sqrt{3}}\,hk_2\right) , \\ k_2 &= f\left( x_n + \frac{\sqrt{3} +1}{2\sqrt{3}}\,h , \ y_n + \frac{\sqrt{3} +2}{4\sqrt{3}}\,hk_1 + \frac{h}{4}\,k_2\right) , \end{align*}

 

First we present code when the step size (h) is given:

Then the same algorithm when the number of steps is specified:

Using Mathematica command NDSolve, we check the obtained value:

 

II. Third Order Runge-Kutta Approximation

---

We start with the following third order algorithm credited to Kutta:

\begin{equation} \label{E35.9} y_{n+1} = y_n + \frac{h}{6}\, ( k_1 + 4 k_2 + k_3 ), \end{equation} where \[ \begin{array}{l} k_1 = f( x_n , y_n ) , \\ k_2 = f( x_n + h/2 , y_n + k_1 h/2 ) , \\ k_3 = f( x_n +h , y_n + 2h k_2 - h k_1 ) . \end{array} \]

If you want to see only the last value (or any other particular value) use its modification:

 

 

III. The Nystrom Algorithm

---

The Nystrom algorithm:

\begin{equation} \label{E35.nystrom} y_{n+1} =y_n + \frac{h}{8}\,(2k_1 + 3k_2 + 3k_3 ), \qquad n=0,1,2,\ldots , \end{equation} where \[ k_1 = f(x_n , y_n ), \ k_2 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_1 \right) , \ k_3 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_2 \right) ; \]

Example.

 

IV. The Nearly Optimal Algorithm

---

The nearly optimal algorithm (used in MATLAB ode23)

\begin{equation} \label{E35.optimal} y_{n+1} =y_n + \frac{h}{9} \,(2k_1 + 3k_2 + 4k_3 ), \qquad n=0,1,2,\ldots , \end{equation} in which the stages are computed according to \[ k_1 = f(x_n , y_n ), \quad k_2 = f \left( x_n + \frac{h}{2} , y_n + \frac{h}{2}\,k_1 \right) , \quad k_3 = f \left( x_n + \frac{3h}{4} , y_n + \frac{3h}{4}\,k_2 \right) ; \]

The same algorithm when the step size is specified:

 

 

V. The Heun Algorithm

---

Heun's algorithm:

\begin{equation} \label{E35.heun} y_{n+1} =y_n + \frac{h}{4} \,(k_1 + 3k_3 ), \qquad n=0,1,2,\ldots , \end{equation} where \[ k_1 = f(x_n , y_n ), \ k_2 = f\left( x_n + \frac{h}{3} , y_n + \frac{h}{3}\,k_1 \right) , \ k_3 = f\left( x_n + \frac{2h}{3} , y_n + \frac{2h}{3}\,k_2 \right) . \]

 

 

VI. Bogacki and Shampine Approximation

---

 

 

VII. Fourth Order Runge-Kutta Approximation

---

Let us introduce the classical fourth order Runge-Kutta algorithm:

\begin{equation} \label{E35.10} y_{n+1} = y_n + \frac{h}{6}\, ( k_1 + 2 k_2 + 2 k_3 + k_4 ), \end{equation} where \[ \begin{array}{l} k_1 = f_n = f(x_n , y_n ), \\ k_2 = f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2} k_1 \right) , \\ k_3 =f\left( x_n + \frac{h}{2}, y_n + \frac{h}{2} k_2 \right) , \\ k_4 = f( x_n + h, y_n + h k_3 ). \end{array} \]

We revisit our example, and solve it in the following way:

Other codes of classical Runge-Kutta method (the first one is when the step size h is specified, and another one when the number of steps is given):

Another version when the number of steps is specified.

Here is a slightly different version of the above code:

We apply this numerical method to find an approximate values for the solution to the IVP with slope function \( f(x,y) = x^2 - y^2 \)

Kutta's algorithm of order four:

\begin{equation} \label{E35.11} y_{n+1} = y_n + \frac{h}{8}\, ( k_1 + 3 k_2 + 3 k_3 + k_4 ), %\eqno{(3.12)} \end{equation} where \[ \begin{array}{l} k_1 = f_n = f(x_n , y_n ), \\ k_2 = f\left( x_n + \frac{h}{3}\, , y_n + \frac{h}{3}\, k_1 \right) , \\ k_3 =f\left( x_n + \frac{2h}{3}\, , y_n - \frac{h}{3}\, k_1 + h k_2 \right) , \\ k_4 = f( x_n + h, y_n + h k_1 - h k_2 + h k_3 ). \end{array} \] If one wants to see the whole list of intermediate values calculated by this algorithm, then replace the last line with the following: Example.

Here is a slightly different version of Kutta's algorithm:

Now the same algorithm, but with the step size specified: If one wants to see the whole list of intermediate values calculated by this algorithm, then replace the last line with the following: Example.

Now we apply the same code with modified last line

Another version:

Gill's method:
\begin{equation} \label{E35.12} y_{n+1} = y_n + \frac{h}{6} \left[ k_1 + 2 \left( 1 - \frac{1}{\sqrt{2}} \right) k_2 + 2 \left( 1 + \frac{1}{\sqrt{2}} \right) k_3 + k_4 \right] , \end{equation} where \[ \begin{array}{l} k_1 = f(x_n , y_n ), \\ k_2 = f\left( x_n + \frac{h}{2}\, , y_n + \frac{h}{2}\, k_1 \right) , \\ k_3 =f\left( x_n + \frac{h}{2}\, , y_n - \left( \frac{1}{2} - \frac{1}{\sqrt{2}} \right) h k_1 + \left( 1 - \frac{1}{\sqrt{2}} \right) h k_2 \right) , \\ k_4 = f( x_n + h, y_n - \frac{1}{\sqrt{2}} h k_2 + \left( 1 + \frac{1}{\sqrt{2}} \right) h k_3 ). \end{array} \]

Another version:

Then we exercute the code:

 

Application of the method:

Another version:

Fifth order Butcher's method:
\begin{equation} \label{E36.butcher} y_{k+1} = y_k + \frac{h}{90} \left( 7k_1 + 32 k_3 + 12 k_4 + 32 k_5 + 7 k_6 \right) , \qquad k=0,1,\ldots , \end{equation} where \[ \begin{align*} k_1 &= f(x_k , y_k ), \qquad &k_2 = f \left( x_k + h/4 , y_k + hk_1 /4 \right) , \\ k_3 &= f \left( x_k + h/4 , y_k + k_1 h/8 + k_2 h/8 \right), \qquad &k_4 = f \left( x_k + h/2 , y_k - k_2 h/2 + k_3 h \right) , \\ k_5 &= f \left( x_k + 3h/4 , y_k + 3k_1 h/16 +9 k_4 h/16 \right) , \qquad &k_6 = f \left( x_k + h , y_k - 3k_1 h/7 + 2k_2 h/7 + 12k_3 h/7 -12k_4 h/7 + 8k_5 h/7 \right) . \end{align*} \]

Another version of Butcher's method:

 

3.5. Multistep Methods

---

Predictor-corrector of order 2:

To apply yjis method, use the following commands:

If you want to see the whole table of values:

Example.

Predictor-corrector method of order 2:

or using Function command

Predictor-corrector of order 3:

Preprocessing includes:
Nsteps = number of steps
h = step size
f[u_,v_] = slope function
Initial values:
x[0] =
y[0] =
Generate first values using a single-step numerical method or actual solution:
y[1] =
y[2] =
Arrays of independent and dependent variables corresponding to uniform grid:
Array[x, Nsteps]; Array[y, Nsteps];
i=1;
While[i < Nsteps + 1, x[i] = x[i - 1] + h; i++];

 

 

Variable coefficient recurence, corresponding to explicit Adams' method of the second order where the first value is determined from the actual solution (however, it can be found using any single-step method, say obtained by classical Runge-Kutta algorithm) can be solved numerically:

 

1.3.6. Accuracy

---

Sometimes you may be interested to find out what methods are being used in NDSolve.
Here you can see the coefficients of the default 2(1) embedded pair.

You also may want to compare some of the different methods to see how they perform for a specific problem.

Implicit Runge--Kutta methods have a number of desirable properties.

The Gauss--Legendre methods, for example, are self-adjoint, meaning that they provide the same solution when integrating forward or backward in time.

http://reference.wolfram.com/language/tutorial/NDSolveImplicitRungeKutta.html

A generic framework for implicit Runge--Kutta methods has been implemented. The focus so far is on methods with interesting geometric properties and currently covers the following schemes:

"ImplicitRungeKuttaGaussCoefficients"
"ImplicitRungeKuttaLobattoIIIACoefficients"
"ImplicitRungeKuttaLobattoIIIBCoefficients"
"ImplicitRungeKuttaLobattoIIICCoefficients"
"ImplicitRungeKuttaRadauIACoefficients"
"ImplicitRungeKuttaRadauIIACoefficients"

The derivation of the method coefficients can be carried out to arbitrary order and arbitrary precision.

 

 

1.3.7. Numerov Method

---