Lagrange Interpolation

#Lagrange interpolation
n=input(" No of data points = ")
xd=[]
yd=[]
print "Enter data points"
for i in range(n):
    xi,yi=input()
    xd.append(xi)
    yd.append(yi)

x=input(" Enter the x value for which y is to be found ")
l=0
for i in range(n):
    lg=1
    for j in range(n):
        if i!= j:
            lg=lg*(x-xd[j])/float(xd[i]-xd[j])
           
    l=l+yd[i]*lg
   

print " Y value = ",l

Newton forward interpolation

#Newton forward interpolation
from numpy import*

def fact(n):
    if n<=1:
        return(1)
    else:
        return(n*fact(n-1))
       
n=input("Enter no of elements ")
y=zeros([n,n])
x0,h=input("Enter the value of x[0] and finite difference ")
x=[x0+i*h for i in range(n)]

print "Enter y values"
for i in range(n):
    print "Enter y value for x[%d]"%(i)
    y[i,0]=input()
for i in range(1,n):
    k=n-i
    for j in range(k):
        y[j,i]=y[j+1,i-1]-y[j,i-1]
   
print y

xp=input(" Enter the x value for which y is to be found")
u=(xp-x[0])/float(h)
#print "u=",u
yp=y[0,0]
for i in range(1,n):
    t=1
    for j in range(i):
        t=t*(u-j)
    yp=yp+t*y[0,i]/fact(i)
   
print "The value of y is  %f "%(yp)

Gamma function value x>0 - Purushothaman sir

http://www.scribd.com/doc/59295714/gamma-python-plot

Pie Chart

from pylab import *

figure(1, figsize=(8,8))
#axes([0.1, 0.1, 0.8, 0.8])

mycolors=['red', 'blue', 'green', 'brown']
mylabels=['Red', 'Blue', 'Green', 'Brown']
fracs=[40, 30, 20, 10]
pie(fracs,labels=mylabels,colors=mycolors)

show()

Parametric curves

from pylab import*
a=5
b=3
wx=4
wy=5
dx=0
dy=0
t=linspace(0,10,200)
x=a*cos(wx*t-dx)
y=b*sin(wy*t-dy)
figure(1,figsize=(4,4))
plot(x,y)
show()

Polar Graphs- Butterfly plot

from pylab import*
theta=linspace(0,8*pi,1000)
r=exp(sin(theta))-2*cos(4*theta)+sin((2*theta-pi)/24.0)
#r=cos(5*theta)*cos(5*theta)+sin(3*theta)+0.3
figure(1,figsize=(10,10))
polar(theta,r,'red')
show()

Bessel function plot

#Bessel function plot
from pylab import*
def fact(n):
    if n<=1:
        return(1)
    else:
        return(n*fact(n-1))
xdat=linspace(0,8*pi,200)

for n in range(5):
    jnx=[]
    for x in xdat:
        sum=0
        for k in range(0,50):
            sum+=(-1)**k*(x/2)**(n+2*k)/(fact(k)*fact(n+k))
        jnx.append(sum)
   
    plot(xdat,jnx)
   
show()

Power series- sine series

#sine series
from pylab import*
def fact(n):
    if n<=1:
        return(1)
    else:
        return(n*fact(n-1))
x=linspace(0.0,8*pi,100)
y=sin(x)
sinp=[]

for k in x:
    sum=0
    for i in range(0,31):
        sum+=(-1)**i*k**(2*i+1)/fact(2*i+1)
    sinp.append(sum)
   
plot(x,y,x,sinp)
show()

Fourier series- Square wave - Sawtooth wave

#Fourier Series-Square Wave
#f(t)=a/2+(2a)/pi[sin(wt)+sin(3wt)/3+sin(5wt)/5+........]

from pylab import*
t=linspace(0.0,2*pi,100)
a=2.0
w=2.0
y=a/2
for i in range(1,200,2):
    y=y+2*a*sin(i*w*t)/(pi*i)
plot(t,y)
show()

#Fourier Series-Saw-tooth wave

#f(t)=a/2+a/pi[sin(wt)+sin(2wt)/2+sin(3wt)/3+sin(4wt)/4+sin(5wt)/5+........]

from pylab import*

t=linspace(0.0,2*pi,100)

a=2.0

w=5.0

y=a/2

for i in range(1,100):

    y=y+a*sin(i*w*t)/(pi*i)

plot(t,y)

show()

Damped Harmonic oscillator

#Damped harmonic oscillator
#dx/dt=v
#dv/dt=-w*w*x-b*v
#Feynman Newton algorithm
from pylab import*
from math import*
x0=10
v0=20
w=2
b=0.05
t=0 #initial time
h=0.01 #time step size

xdat=[]
vdat=[]#velocity data store here
time=[]#time stored here

v0=v0+h*(-w*w*x0-b*v0)/2

while(t<=250):
   
    x1=x0+h*v0
    v1=v0+h*(-w*w*x1-b*v0)
   
    xdat.append(x1)
    vdat.append(v1)
    time.append(t)
    x0=x1
    v0=v1
    t=t+h
   
figure(1)
title("Damped harmonic oscillator motion")
xlabel(" Time")
ylabel("x")
plot(time,xdat)
grid(True)
figure(2)
title("Damped harmonic oscillator motion")
xlabel("Time")
ylabel("Velocity")
plot(time,vdat)
grid(True)
figure(3)
title("Damped harmonic oscillator motion")
xlabel("x")
ylabel("velocity")
plot(xdat,vdat)
grid(True)

show()

Projectile motion

#projectile motion
from pylab import*
from math import*
ax=0
ay=-9.8
theta=30
v=20
vx0=v*cos(theta*pi/180)
vy0=v*sin(theta*pi/180)
x0=0
y0=0
t=0 #initial time
h=0.001 #time step size

xdat=[x0]
ydat=[y0]#position data stored here
vdat=[v]#velocity data store here
time=[t]#time stored here

while(y0>=0):
    x1=x0+h*vx0
    y1=y0+h*vy0
    vx1=vx0+h*ax
    vy1=vy0+h*ay
    xdat.append(x1)
    ydat.append(y1)
    vdat.append(sqrt(vx1*vx1+vy1*vy1))
    time.append(t)

    x0=x1
    y0=y1
    vx0=vx1
    vy0=vy1
    t=t+h

figure(1)
title("Motion of Projectile")
xlabel(" x")
ylabel("y")
plot(xdat,ydat,"r-")


figure(2)
title("Motion of projectile")
xlabel(" Time")
ylabel("Velocity")
plot(time,vdat,"k-")
show()


   

Harmonic oscillator

#Harmonic oscillator
#dx/dt=v
#dv/dt=-w*w*x
#Feynman-Newton method
#x(t+h)=x(t)+h*v(t+h/2)
#v(t+h/2)=v(t-h/2)+h*a(t)
#v(1/2)=v0+(h/2)a0
from pylab import*
from math import*
x0=5
v0=5
w=1
t=0 #initial time
h=0.001 #time step size

xdat=[]
vdat=[]#velocity data store here
time=[]#time stored here
v0=v0+h*(-w*w*x0)/2
while(t<=100):
    x1=x0+h*v0
    v1=v0+h*(-w*w*x1)
   
   
    xdat.append(x1)
    vdat.append(v1)
    time.append(t)
    x0=x1
    v0=v1
    t=t+h
   
figure(1)
title("Harmonic oscillator motion")
xlabel(" Time")
ylabel("x")
plot(time,xdat)
grid(True)

figure(2)
title("Harmonic oscillator motion")
xlabel("Time")
ylabel("Velocity")
plot(time,vdat)
grid(True)

figure(3)
title("Harmonic oscillator motion")
xlabel("x")
ylabel("Velocity")
plot(xdat,vdat)
grid(True)

show()

Advanced Python programming notes by Purushothaman sir

Advanced Python programming

Simulation falling through viscous medium

#Falling of a body through viscous medium
#Loss of weight   V*d1*g-V*d2*g=V*d1*g(1-d2/d1)=m*a where a=g(1-d2/d1)
#Drag in the case of ball 6*pi*eta*a*v=k*v
#dy/dt=v
#dv/dt=a-kv
from pylab import*
from math import*
d1=2000.0
d2=750.0
g=9.8

height=2000
y0=0
v0=0
a=(1.0- d2/d1)*g #including buoyant force
k=0.2 #drag coefficient
h=0.001
t=0
y=[y0]
v=[v0]
time=[t]
 
while height-y0>0:
    v1=v0+h*(a-k*v0)
    y1=y0+h*v0
 
    t=t+h
    y.append(height-y1)
    v.append(fabs(v1))
    time.append(t)
    y0=y1
    v0=v1

 
figure(1)
title("Fall through viscous medium")
xlabel(" Time")
ylabel("Position")
plot(time,y)


figure(2)
title("Fall through viscous medium")
xlabel(" Time")
ylabel("Velocity")
plot(time,v)
show()


 

Simulation free fall of a body

#Free fall of a body
#We are solving two first order differential equations by Euler method
#1. dy/dt=v initial condition y0
#2. dv/dt=a initial condition v0


from pylab import*
from math import*

a=-9.8 #acceleration 
v0=0.0 #initial velocity
y0=100 #initial height
t=0 #initial time
h=0.01 #time step size
ydat=[y0]#position data stored here
vdat=[v0]#velocity data stored here
time=[t]#time stored here

#Implementation Euler method
while(y0>0):
    y1=y0+h*v0
    v1=v0+h*a
    t=t+h
    ydat.append(y1)
    vdat.append(fabs(v1))
    time.append(t)
    y0=y1
    v0=v1
   
figure(1)
title("Free fall of a body")
xlabel(" Time")
ylabel("Position")
plot(time,ydat)


figure(2)
title("Free fall of a body")
xlabel(" Time")
ylabel("Velocity")
plot(time,vdat)
show()

Python notes by Pramode- See the link

http://radiantbytes.com/books/python-latex/

Matrix operations using numpy

# This example shows how a matrix can be read from terminal and to do various operations
from numpy import *
print "Enter order of matrix"
m=input("Rows=")
n=input("Columns=")
xx=zeros([m,n],float)
print "enter matrix elements"
for i in range(m):
    print "Enter elements of row ",i+1
    for j in range(n):
        xx[i][j]=input()

print "xx \n",xx
print
print "xx+2 \n",xx+2
print
print "xx+xx \n",xx+xx
print
print "xx*xx \n",xx*xx
print
print "xx-xx \n",xx-xx
print
print "xx*2 \n",xx*2
print
print "xx/xx \n",xx/xx

Bisection method

# Bisection method

def f(x):
return (x**2-2)
e=0.00001
a=float(input('Enter limit a   '))
b=float(input('Enter limit b   '))
if f(a)*f(b)<=0.0:
while abs(a-b)>=e:

x=(a+b)/2  
if(f(x)*f(a)>0):
a=x
else:
b=x
print'The best approximation is %.5f'%x
else:
print' sorry !!! There is no root  between',a,'and ',b

Fixed Iteration method

# Fixed iteration method to find approximate root
from math import sqrt
#function definition
def f(x):
    f1=1.0/sqrt(1+x)
    return f1

e=0.0001
x0=float(input("Enter starting value   "))

while abs(f(x0)-x0)>=e:
    x0=f(x0)
  
  
print 'The best approximation to root is %.4f'%(x0)

Newton-Raphson method

#Newton-Raphson method to find approximate root

from math import *
e=0.00001
#function definition
def f(x):
    return x*sin(x)+cos(x)

#differential of the function
def fp(x):
    return x*cos(x)


x0=float(input("Enter the approximate root \t"))

while 1:
    x1=x0-f(x0)/fp(x0)

    if abs(x0-x1)<e:
        break
    else:
        x0=x1

print "The best approximation to root is   %.5f"%x1

Functions

#program to calculate sin(x) using two functions power(x,n) and fact(n)
from math import pi
def power(x,n):
    if n==1:
        return x
    return x*power(x,n-1)
def fact(n):
    if n==0:
        return 1
    return n*fact(n-1)

x=raw_input("Enter angle in degrees  ")
xr=float(x)*pi/180.0
term=xr
sum=0
n=3
sign=1
while term>0.00001:
    sum+=sign*term
    term=power(xr,n)/fact(n)
    n+=2
    sign=-sign
print "sin(",x,")=",sum

Sorting in multiple fields using lists of tuples

#list of tuples containing data of solids
#each tuple contains name of substance, density(10**3), Young's modulus(10**10), Rigidity modulus(10**10)
solids=[('Aluminium',2.7,6.9,2.4),('Brass',8.6,9,3.45),('Copper',8.89,11,3.4),\
        ('Glass',2.6,6,3.0),('Iron',7.5,10,3.5)]

print "Substance\t\t\tDensity\t\t\tYoung's modulus\t\t\tRigidity modulus"
print
for item in sorted(solids,key=lambda a:(-a[1],a[2]),reverse=True):
    print item[0],"\t\t\t",item[1],"\t\t\t\t",item[2],"\t\t\t\t",item[3]

Dictionary example

#program to determine resistance of a resistor using colour bands--uses dictionary
color1={'black':0,'brown':1,'red':2,'orange':3,\
       'yellow':4, 'green':5,'blue':6,'violet':7,\
       'gray':8,'white':9}
color2={'black':1e0,'brown':1e1,'red':1e2,'orange':1e3,\
       'yellow':1e4, 'green':1e5,'blue':1e6,'violet':1e7,\
       'gray':1e8,'white':1e9,'gold':1e-1,'silver':1e-2}
color3={'brown':1,'red':2,'green':0.5,'blue':0.25,\
        'violet':0.1,'gray':0.05,'gold':5,'silver':10,'none':20}
resist=0
while 1:
    b1=str(raw_input("Enter colour first band = ")).lower()
    if b1 in color1:
        resist=color1[b1]
        break
    else:
        print "Enter a valid color"
        continue
   
resist *=10
while 1:
    b2=str(raw_input("Enter colour second band = ")).lower()
    if b2 in color1:
        resist +=color1[b2]
        break
    else:
        print "Enter a valid color"
        continue

while 1:
    b3=str(raw_input("Enter colour third band = ")).lower()
    if b3 in color2:
        resist *=color2[b3]
        break
    else:
        print "Enter a valid color"
        continue

while 1:
    b4=str(raw_input("Enter colour fourth band = ")).lower()
    if b4 in color3:
        tol=color3[b4]
        break
    else:
        print "Enter a valid color"
        continue


print resist," +-",tol,"%"


      

Lists---Tuples---Dictionaries

Lists are what they seem - a list of values. Each one of them is numbered, starting from zero - the first one is numbered zero, the second 1, the third 2, etc. You can remove values from the list, and add new values to the end. Example: Your many cats' names.

"A list"

cats = ['Tom', 'Snappy', 'Kitty', 'Jessie', 'Chester']

"If you want to print for example the first name on the list you would do the following"

print cats[0]

"Also if you want to add a name to the list you would do this"

cats.append('Catherine')

"To remove a name you would do this"

del cats[0]

"A list" cats = ['Tom', 'Snappy', 'Kitty', 'Jessie', 'Chester'] "If you want to print for example the first name on the list you would do the following" print cats[0] "Also if you want to add a name to the list you would do this" cats.append('Catherine') "To remove a name you would do this" del cats[0]

Tuples are just like lists, but you can't change their values. The values that you give it first up, are the values that you are stuck with for the rest of the program. Again, each value is numbered starting from zero, for easy reference. Example: the names of the months of the year.

"Here is an example of a tuple"

months = ('January','February','March','April','May','June',\
'July','August','September','October','November','  December')

"So that’s the difference between lists and tuples. A list can be modified but a tuple cannot be modified in anyway(unless changed in the source code of the program)"

"Here is an example of a tuple" months = ('January','February','March','April','May','June',\ 'July','August','September','October','November',' December') "So that’s the difference between lists and tuples. A list can be modified but a tuple cannot be modified in anyway(unless changed in the source code of the program)"

There are also dictionaries

Dictionaries are similar to what their name suggests - a dictionary. In a dictionary, you have an 'index' of words, and for each of them a definition. In python, the word is called a 'key', and the definition a 'value'. The values in a dictionary aren't numbered - tare similar to what their name suggests - a dictionary. In a dictionary, you have an 'index' of words, and for each of them a definition. In python, the word is called a 'key', and the definition a 'value'. The values in a dictionary aren't numbered - they aren't in any specific order, either - the key does the same thing. You can add, remove, and modify the values in dictionaries. Example: telephone book.

"Here is a simple phonebook"

phonebook = {'Andrew Parson':8806336, \
'Emily Everett':6784346, 'Peter Power':7658344, \
'Lewis Lame':1122345}

"To add a name to the phonebook you would do the following"

phonebook['Gingerbread Man'] = 1234567

"If you want to remove a name you would do the same as you would do in a list"

del phonebook['Andrew Parson']

daniweb IT discussion community

Value of pi using Monte Carlo Method

#To determine value of pi using Monte Carlo Method
from random import *
j=0
for i in range(10000000):
    x=random()
    y=random()
    if (x**2+y**2)<=1:j+=1
print "Value of pi = ",4.0*j/i

by Prof. Purushothaman Vadakkedam

Trapezoidal and Simpson's rule -Integration

#This program calculates integral of 1/1+x in the interval 0 to 1
#using both trapezodal and simpsons rule using Lists without
#importing any modules
h=.001
N=1000
x=[h*i for i in range(N)]
y=[1/(1+i) for i in x]
trapezoidal=h/2*(y[0]+y[-1]+2*

sum(y[1:-2]))
print "tapezoidal rule gives%.6f"%trapezoidal
simpsons=(h/3)*(y[0]+y[-1]+4*sum(y[1:-2:2])+2*sum(y[2:-3:2]))
print "simpsons rule gives %.6f"%simpsons

 

by Prof. Ahamed Ashraf 

 

Simple version of orthogality problem

from math import *
dotproduct=[sin(pi*i/180)*cos(pi*i/180)for i in range(180)]

print sum(dotproduct)

by Prof. Ahamed Ashraf

Arrays-2

#To verify orthogonality of sine and cosine functions
#Principle: Create two vctors X and Y using the values of sine
#and cosine for the same angle within relevant range.
#Show that sum of products xi.yi=0'
import Numeric, math
sum=0
for i in range(180):
    theta=math.pi*i/180
    x=math.cos(theta)
    y=math.sin(theta)
    sum+=x*y
print "The sum of products xi.yi=%e is negligible. Hence orthogonal in the range (0-180)degrees"%sum


by Prof. Purushothaman Vadakkedam