## Taylor method with automatic differentiation – 2

In the previous post, we have talked about Taylor method to solve ordinary differential equation (ODE) numerically. All the derivatives in the Taylor method are computed using automatic differentiation. However, the drawback of this method is when we change our ODE, we have to write another code. In this post, we are going to present a way such that when we want to compute another ODE, we don’t need to change the majority of our code.

Continuing our analysis in the last post, we are going to need operations on polynomials. We are going to define operations such as:

- Addition between a constant and a polynomial
- Multiplication between a constant and a polynomial
- Addition between polynomials
- Multiplication between polynomials
- Dividing a constant by a polynomial
- Dividing a polynomial by a polynomial
- Logarithm of a polynomial
- Exponentiation of a polynomial
- Power of a polynomial

Operation 1 and 2 are trivial. Suppose is a polynomial, then for some real constant . Suppose is a real constant, then the addition between and is a polynomial where and for all . The multiplication between and is a polynomial where for all .

Operation 3 and 4 are not difficult. Suppose are polynomials. Then the addition between and is a polynomial where for all . The multiplication between and is a polynomial where .

Operation 5 and 6 are a little bit tricky. Suppose is a polynomial and is a constant. We want to compute where . Notice that . Therefore, we are going to take advantage of this equality to find the coefficient of . Thus we shall have , and . For operation 6, let’s suppose . As before, notice that , thus we can compute the coefficient of using this equality. Thus, and .

Now, let’s consider the next operation, which is a logarithm of a polynomial. Suppose is a polynomial and we are going to find . Let’s differentiate both side with respect to , then we have . Thus by exploiting a division operation above, we have and .

Now, let’s consider the exponential operator, which is an exponentiation of a polynomial. Suppose is a polynomial and we are going to find . Let’s take logarithm on both side and then differentiate it, we then have , or . Therefore, by exploiting a multiplication operation, we shall have and .

For the last operation of this post, we shall consider the power of a polynomial. Suppose is a polynomial, we are going to find for any real number . First, we need to take logarithm on both sides and then differentiate the result with respect to , or . By exploiting multiplication and division, we will get the coefficient of .

To facilitate these operation, we build a Python module that will handle these operations. We will also write another code to integrate an ODE. This integration code will work in general without any dependencies on the system or on the order of the system.

def integrate(F,y_init,t_init = 0.0 ,tStop = 1.0 ,h = 0.05 ): dim = len(y_init) y_sol = [] for i in range(0,dim): y_sol.append([]) y_sol[i].append(y_init[i]) t_sol = [] t_sol.append(t_init) t = t_init while t < tStop: y = [] for i in range(0,dim): y.append([]) y[i].append(y_sol[i][-1]) for k in range(0,20): dy = F(y) for i in range(0,dim): y[i].append(dy[i][k]/float(k+1)) n = len(y[0]) for i in range(0,dim): temp = 0.0 for j in range(0,n): temp = temp + y[i][j]*h**j y_sol[i].append(temp) t_sol.append(t+h) t = t + h return np.array(t_sol) ,np.transpose(np.array(y_sol))

All the operations above and the last code will be written in a python file so that if we need to use this Taylor method, we can just import it. Let’s take an example to integrate one ODE. We will create another code that specifies the ODE and also imports all the operations above and then numerically integrate it.

from mytaylor import * import matplotlib.pyplot as plt import numpy as np import math as m #define your ODE here def fun(y): z = subtractpoli(y[0],multiplypoli(y[0],y[0])) return [z] #define your computation parameter here y_init = [0.5] t_init = 0 h = 0.1 tStop = 5.0 #integrating ODE T, Y = integrate(fun,y_init,t_init,tStop,h) #real solution computation real_sol = [] for t in T: temp = m.exp(t)/(m.exp(t)+1.0) real_sol.append(temp) #plotting the error versus time real_sol = np.array(real_sol) error = Y[:,0]-real_sol plt.plot(T[:],error[:]) plt.show()

In the above code, we try to numerically integrate with initial condition . This is the same system we try to integrate in the previous post. I hope you can see the difference between the last code and the code we have shown in the previous post. We really don’t need to consider about Taylor method at all, we just need to write our system in the above code if we want to compute the numerical solution of our code. Notice that we have written all the operations and the integrating code in a file named mytaylor.py which is imported in the first line of the above code. In the fifth line, we can write any ODE we want to integrate, however, the way we write our ODE will be slightly different. Due to polynomial operations that we have defined in this post, cannot be written as just but addpoly(). Once we define our ODE, we can just integrate it like in line 15.

Another thing we would like to do is to extend the operation of polynomial. We have not included trigonometry operations on polynomials and et cetera. We would like also to test this Taylor method to a stiff system. Please find the following two files (mytaylor.py and test_case.py) for interested readers if you would like to test yourselves.

## Your says