Archive

Posts Tagged ‘ictp’

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_1a_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_0a_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. error_plot

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()

Puzzle on the road – how many combinations of food that I can have ?

July 7, 2011 Leave a comment

Let me start with the introduction of why I come up with this puzzle. I was invited to go to Trieste, Italy for a summer school in ICTP (click here for information of ICTP). I was lucky cause ICTP sponsored me so I did not have to pay for my living cost while I was there. Everyday, I ate in a cafeteria and I did not need to pay, cause ICTP gave me coupons that I can use to buy my lunch and my dinner. Cafeteria usually provides a set of first course meal, which is usually a pasta, a set of second course meal which is usually either beef, chicken, fish or vegetarian meal and a set of cooked vegetable.

One full meal that I can buy with a coupon,  consists of one of the first courses, one of the second courses and two of the cooked vegetables. In the first course meal, there are three choices of pasta and I am only allowed to take one of these. In the second course meal, there are also three choices  whether you want, say beef or fish or vegetarian meal and again I am allowed to take one of these. In the set of vegetable, there are three choices and I am allowed to take two of these three choices.

The question is how many combinations of meal that I can have.

To make it more interesting, let us assume that these following situations can occur.  I am sometimes so hungry so I take the full meal. But there is another time when I do not feel like eating so I do not eat. There is also another time when I want to eat but actually I am not that hungry so I only choose the first course and one of the cooked vegetable.

Anyway, it was a nice experience for me to be able to visit ICTP. I made a lot of friends for different countries there. I gained a lot of information that is really really useful for my research.

Here is a picture taken from my cell phone while I was there.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ps. I am not sure what is the correct answer. I have had already a number in my mind and let us confirm whether my guess is correct.

Categories: math puzzle Tags: , , , ,