### Archive

Posts Tagged ‘C++’

## Taylor method with automatic differentiation

June 14, 2016 1 comment

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:

$\frac{dy}{dt} = f(y)$ and $y(t) = y_0$.

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

$y(t+h) = y(t) + y'(t) h + y''(t)\frac{h^2}{2} +y'''(t)\frac{h^3}{6}+\dots$.

Therefore, depending on how accurate we want our solution is, all we need is to compute the derivatives of $y(t)$. 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 $y(t+h)$ so that we can compute $y(t+h)$ as accurately as possible. Let us consider the series expansion of $y(t+h)$ as above, but I will use constant $a_i$ to replace the $i-$th derivative of $y(t)$,

$y(t+h) = a_0 + a_1 h + a_2 h^2 + \dots$.

Let us derive the above series with respect to $h$.

$y'(t+h) = a_1 + 2 a_2 h + 3a_3 h^2 + \dots$.

Notice that the derivative at the left hand side of the above equation is actually equal to $f(y(t+h))$ since we know that $y(t+h)$ satisfies the ordinary differential equation. Then we have,

$f(y(t+h)) = a_1 + 2 a_2 h + 3a_3 h^2 + \dots$.

Since $y(t+h)$ is a series of $h$ then $f(y(t+h))$ is also a series of $h$. Then both sides of the above equation are series of $h$, thus we can equate each of the coefficient of $h^j$ to get $a_1$$a_2$ and so on. The most difficult task in this problem is to find the series of $f(y(t+h))$ since the function $f(y)$ can be anything.

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

$\frac{dy}{dt} = y(1-y)$ and $y(0)=0.5$.

Let us choose $h=0.1$ and compute $y(0.1)$. As before, we assume that $y(0.1)$ is a series of $h$ as follows,

$y(0.1) = a_0 + a_1 h + a_2 h^2 + \dots$.

Differentiating both sides of the above equation, we get

$y'(0.1) = a_1 + 2 a_2 h + 3 a_3 h^2 + \dots$.

The right hand side of the above is equal to

$y'(0.1) = y(0.1)(1-y(0.1)) = (a_0 + a_1 h + \dots)(1 - (a_0 + a_1 h + \dots)) = b_0 + b_1 h + b_2 h^2 + \dots$.

Therefore, we get,

$a_1 = b_0$$a_2 = b_1/2$,$a_3 = b_2/3$,  and so on.

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

$b_0 = a_0 \times (1-a_0)$, thus $a_1 =a_0 \times (1-a_0)$
$b_1 = (a_0 \times -a_1 ) + (a_1 \times (1-a_0))$, thus $a_2 =((a_0 \times -a_1 ) + (a_1 \times (1-a_0)))/2$
$b_2 = (a_0 \times -a_2 ) + (a_1 \times (-a_1))+(a_2 \times (1-a_0))$, thus $a_3 =((a_0 \times -a_2 ) + (a_1 \times (-a_1))+(a_2 \times (1-a_0)))/3$

and so on.

The more $a$‘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 $a$ 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()

#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); }