## Taylor method with automatic differentiation

I teach numerical analysis course and in the one of the chapters, I teach my students about the numerical solution of ordinary differential equations (ODEs). One of numerical solution methods I teach is the so-called Taylor method. This method works by applying the Taylor series to the solution of the differential equation. Suppose we have one-dimensional differential equation as follow:

and .

Let’s assume that we only consider an autonomous ODE first. To find a solution at the next step, , we need to apply Taylor expansion, as follow:

.

Therefore, depending on how accurate we want our solution is, all we need is to compute the derivatives of . Unfortunately, computing all the derivatives is not an easy task for a computer. We have to derive all the derivatives by ourselves and then input them to the computer.

In this note, I will introduce an automatic differentiation as a tool to compute the derivatives of so that we can compute as accurately as possible. Let us consider the series expansion of as above, but I will use constant to replace the th derivative of ,

.

Let us derive the above series with respect to .

.

Notice that the derivative at the left hand side of the above equation is actually equal to since we know that satisfies the ordinary differential equation. Then we have,

.

Since is a series of then is also a series of . Then both sides of the above equation are series of , thus we can equate each of the coefficient of to get , and so on. The most difficult task in this problem is to find the series of since the function can be anything.

Let’s take an example. Suppose we want to integrate the following differential equation,

and .

Let us choose and compute . As before, we assume that is a series of as follows,

.

Differentiating both sides of the above equation, we get

.

The right hand side of the above is equal to

.

Therefore, we get,

, ,, and so on.

The real question is how to find , and so on. If the function only involves operations such as addition, subtraction and multiplication, it would be easy to find the ‘s. But if the function has operations such as division, exponential, logarithm and trigonometric functions then we have to think a bit harder on how to find ‘s. We go back to the previous example. Since the operations in this example are only multiplication and subtraction, we can easily compute , and so on. Using algebraic modification, we get the following,

, thus

, thus

, thus

and so on.

The more ‘s we compute, the more accurate our solution is. Using Taylor method this way, we can easily apply this method in a program. It should be clear why this method is called automatic differentiation because we don’t need to manually derive all the derivatives.

Next, I tried to numerically integrate the previous example. Below is the graphic of the error of Taylor method with automatic differentiation versus time. As we can see from the figure, that the error of our numerical integrator is so small and it is almost exact. The number of I used in this graph is 10.

I use python to numerically integrate this example. You can find the code to produce the above figure in the end of this note. One of the drawbacks of this method is when you change your ODE, you have to change your code as well. In other numerical integrators such as Runge Kutta, we don’t need to change the whole code, we just need to change the ODE and we should get the solution. Together with my student, I am trying to tackle this problem and I hope it will be posted in here.

import math as m import numpy as np import pylab as plt #preliminary constants order = 10 h = 0.1 t = 0.0 t_end = 5 y0 = 0.5 #defining solution sol = [] sol.append(y0) time = [] time.append(t) #start of the program while t<t_end: a = [] a.append(sol[-1]) for k in range(1,order): if k == 1: a.append(a[k-1]*(1-a[k-1])) else: sum = 0.0 for j in range(0,k-1): sum = sum + a[j]*(-a[k-1-j]) a.append((sum+a[k-1]*(1.0-a[0]))/k) sum = 0.0 for i in range(0,order): sum = sum + a[i]*h**i sol.append(sum) t = t + h time.append(t) #end of the program #real solution computation real_sol = [] for t in time: temp = m.exp(t)/(m.exp(t)+1.0) real_sol.append(temp) #plotting the error versus time sol = np.array(sol) real_sol = np.array(real_sol) time = np.array(time) error = sol-real_sol plt.plot(time[:],error[:]) plt.show()

## Files output in C++

All the programs that we normally use inputs from keyboard and give outputs on the screen. However, in a large scale computation, this situation is not effective. Thus, I have learnt how to produce outputs to a file, so that it can handle a large amount of data. I learned this back in 2007 and I have already forgotten about this code. Now, if I forget about this I know where to look.

Below I wrote a bit of source code that has to be compiled in C++. This program evaluates the exponential of a number from 0 to 3 using four different ways. The first three are using the Taylor polynomial of degree 1, 2 and 3 and the fourth way is using the built in exponential function from C++ itself. The output of this program is a file called “data.dat” that consists of five columns of numbers that can be used for plotting or other programs.

`#include <iostream>`

#include <fstream>

#include <cmath>

using namespace std;

double taylor(double x, int n);

int factorial(int n);

int main()

{

double number;

ofstream out_stream;

out_stream.open("data.dat");

for (number = 0 ; number <= 3.01 ; number = number + 0.01) {

out_stream << number << " ";

out_stream << taylor(number,1) << "\t";

out_stream << taylor(number,2) << "\t";

out_stream << taylor(number,3) << "\t";

out_stream << exp(number) << "\n";

}

out_stream.close();

return 0;

}

double taylor(double x, int n)

{

double y=1;

for(int k=1; k<=n; k++)

{

y = y + pow(x,k)/factorial(k);

}

return y;

}

int factorial(int n)

{

if (n <= 1)

return 1;

else

return n*factorial(n-1);

}

## Your says