Archive

Posts Tagged ‘analisis numerik’

Taylor method with automatic differentiation – 2

July 20, 2016 Leave a comment

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:

  1. Addition between a constant and a polynomial
  2. Multiplication between  a constant and a polynomial
  3. Addition between polynomials
  4. Multiplication between polynomials
  5. Dividing a constant by a polynomial
  6. Dividing a polynomial by a polynomial
  7. Logarithm of a polynomial
  8. Exponentiation of a polynomial
  9. Power of a polynomial

Operation 1 and 2 are trivial. Suppose x(h) is a polynomial, then x(h)=\sum x_i h^i for some real constant x_i. Suppose \alpha is a real constant, then the addition between \alpha and x(h) is a polynomial y(h) where y_0 = x_0 + \alpha and y_i = x_i for all i \geq 1. The multiplication between \alpha and x(h) is a polynomial z(h) where z_i = \alpha x_i for all i.

Operation 3 and 4 are not difficult. Suppose x(h), y(h) are polynomials. Then the addition between x(h) and y(h) is a polynomial z(h) where z_i = x_i + y_i for all i. The multiplication between x(h) and y(h) is a polynomial w(h) where w_i = \sum_{j=0}^{i} x_j y_{i-j}.

Operation 5 and 6 are a little bit tricky. Suppose x(h) is a polynomial and \alpha is a constant. We want to compute y(h) where y(h) = \alpha / x(h). Notice that x(h)y(h)=\alpha. Therefore, we are going to take advantage of this equality to find the coefficient of y(h). Thus we shall have y_0 = \alpha / x_0, and y_i = - (\sum_{j=1}^{i} x_j y_{i-j})/x_0. For operation 6, let’s suppose z(h) = y(h) / x(h). As before, notice that x(h)z(h) = y(h), thus we can compute the coefficient of z(h) using this equality. Thus, z_0 = y_0 / x_0 and z_i = (y_i - \sum_{j=1}^{i} x_j z_{i-j} )/x_0.

Now, let’s consider the next operation, which is a logarithm of a polynomial. Suppose x(h) is a polynomial and we are going to find y(h)=\log(x(h)). Let’s differentiate both side with respect to h, then we have y'(h) = x'(h) / x(h). Thus by exploiting a division operation above, we have y_0 = \log x_0 and y_i = (ix_i - (\sum_{j=1}^{i-1} x_j ((i-j)y_{i-j})))/(ix_0).

Now, let’s consider the exponential operator, which is an exponentiation of a polynomial. Suppose x(h) is a polynomial and we are going to find y(h)=\exp(x(h)). Let’s take logarithm on both side and then differentiate it, we then have y'(h)/y(h) = x'(h), or y'(h)=y(h)x'(h). Therefore, by exploiting a multiplication operation, we shall have y_0  = \exp x_0 and y_i = (\sum_{j=0}^{i-1} y_j ((i-j)x_{i-j}))/i.

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

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 y' = y(1-y) with initial condition y(0)=0.5. 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, x+y cannot be written as just x+y but addpoly(x,y). 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.

My first Coursera course

December 3, 2014 Leave a comment

In the middle of this year, I received a Fulbright Visiting Scholar that allowed me to go to Athens, OH, US. I spent four months in Ohio University starting in September, which means that I took a leave in the university I am working. One of my activities while I was in Ohio University was I took a course from Coursera. At that time I was eager to learn and to know python programming language therefore I decided to take the course, An Introduction to Interactive Programming in Python, offered by Rice University taught by the following professors: Joe Warren, Scott Rixner, John Griener, and Stephen Wong.

My first coursera was nothing but really, really good. The class could not definitely be better. The intructors have put so much effort in preparing the class and at the same time they were really enjoyable and fun. I really learned a lot from this class and I recommend everyone to take this class. After nine weeks, I finished the course and proudly said that I got 96.9% with distinction mark.

But what I wanted to tell here is the content of the course. From knowing nothing of Python I now know that the word ‘python’ is originally from the Month Python’s Flying Circus not the reptilia. I also know how useful this new programming language is. Python can be used to build an interactive game, to conduct some scientific computations, to analyse data science and many more. There are already many communities that use Python and they are (still) growing. One more good thing is the fact that it is free. Python works on either Windows, Linux or Mac and the installation on each device is not that difficult. I am using both Windows and Ubuntu. For Windows, I use mainly Python(x,y) and for Ubuntu, I don’t have to do anything as it is already inside. However, in this course, I did not have to use all of these as one of the instructors has built a nice Python web application that is called Codeskulptor. It will allows us to make a python scripts on your web browser so actually you don’t have to install python in your computer.

From the course, I learned to make interactive games. Even though the games I made were not perfect but I think it is playable and fun enough. I learned to make games like rock-paper-scissors-spock-lizzard, number guessing game, stopwatch game, pong, memory, and asteroids. The following links are the python scripts for all of the above games I made during the course. The link will bring you to a codeskulptor website, in which I had already written a python script. You just need to press the play button on the top left of your screen.

1. RPSLS  This game is called rock-paper-scissors-lizzard-spock game. It is an extention of rock-paper-scissors and first appeared in The Big Bang Theory series.

2. Number guessing game This game, as you might guess, is a number guessing game. The computer hide a number between a small number and a big number and you have finite tries to guess what number the computer is hiding.

3. Stopwatch This game will train your reflexes. It will show a stopwatch and you have to press the pause button every whole second.

4. Pong Oh, please tell me you know this game.

5. Memory There are 16 cards facing down, and all of them come in pair. Figure out the pair with the smallest number of tries.

6. Blackjack The blackjack game, yes, it is the blackjack game that you always know.

7. Asteroids Asteroids is an arcade space shooter game and I think I played this game when I was a kid. I played this game using an Atari device.

Other than those games, I made two another python scripts just for fun. The first is a fan and the second is the illustration of Hilbert’s curve. Hope you like it.

Result of Numerical Analysis mid-semester exam announced

April 7, 2009 8 comments

After long and long examining almost a hundred exam, I finally finish the marking.

Please go to my teaching page and see the result of the mid-semester exam result.

Enjoy guys.

Any comments are welcome.