# Sum of Digits of 87694

##########################

#

# External input the interger

x=int(input("Give any nuber"))

# finding the number of digitis of the input

x2=x

i=0

while x2>1 :

x2=x2/10;

i=i+1

#print(i)

sumd=0.0

#

while i>=0 :

sumd=sumd+ (x-x%10**i)/10**i

x=x%(10**i)

i=i-1

print('sum of digit of is ',sumd)

# Bubble Short of a List of Numbers

##########################

#

# Reading List Data

input1=input("Entre list of elements with space : ") #Read the data as sting

input2=input1.split() # convert the string to the list of string by replacing space by comma

data=list(map(int,input2)) # map string element to Integer

#print(data)

n=len(data)

# Main Body

for i1 in range(n,0,-1):

for i2 in range(1,i1) :

if data[i2]< data[i2-1]:

data[i2], data[i2-1] = data[i2-1], data[i2]

#data Print

print(data)

# Intersection Short of a List of Numbers

##########################

#

# Reading List Data

input1=input("Entre list of elements with space : ") #Read the data as sting

input2=input1.split() # convert the string to the list of string by replacing space by comma

data=list(map(int,input2)) # map string element to Integer

#print(data)

n=len(data)

# Main Body

for i1 in range(1,n):

if data[i1]<data[i1-1]:

for i2 in range(i1,0,-1):

if data[i2]<data[i2-1]:

data[i2], data[i2-1] = data[i2-1], data[i2]

#data Print

print(data)

# variance

##################

x=[-2,0,1,6,4]

sum1=0.0

sum2=0.0

N=len(x)

#print(N)

####################

for i in range(0,N) :

sum1=sum1+x[i]

sum2=sum2+(x[i])**2

print('Average = ',sum1/N, ' Variance = ', sum2/N - (sum1/N)**2)

#

# Trapizoidal rule for Integration

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5]; #cos(x)

y=[1.0,0.995,0.98007,0.95537,0.92106,0.87758,0.82534,0.76484,0.69671,0.62161,0.5403,0.4536,0.36235,0.2675,0.16996,0.070737];

# Initialization

n=len(x); # number of loop

#print(n)

summ=0.0;

h=x[2]-x[1]

###############################

#

# Main body of Program

#

for i in range(1,n-1) :

summ=summ+y[i]

# end of for Loop

#

intg=h*(y[0]+2*summ+y[n-1])/2;

#

###### OutPut ################

#

print("the value of integration is ", intg)

#

# Simpson 1/3 rule for Integration

# Note that this rule is effective for n is

# devisable by 2, i.e., y(0) to y(n:devisable by 2)

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4]; #cos(x)

y=[1.0,0.995,0.98007,0.95537,0.92106,0.87758,0.82534,0.76484,0.69671,0.62161,0.5403,0.4536,0.36235,0.2675,0.16997];

# Initialization

n=len(x); # number of loop

sum2O=0.0;

sum2E=0.0;

h=x[2]-x[1]

###############################

#

# Main body of Program

#

for i in range(1,n-1) : # number element statrs i=0 and ends i=n-1

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sum2E=sum2E+y[i]

else : # odd for 2 divisable condition,

sum2O=sum2O+y[i]

# end of for Loop

#

intg=(h/3)*(y[0]+4*sum2O+2*sum2E+y[n-1]); # 1st element 1=0 and last element i=n-1

#

###### OutPut ################

#

print("the value of integration is ", intg)

#

# Simpson 1/3 rule for Integration

# Note that this rule is effective for n is

# devisable by 2, i.e., y(0) to y(n:devisable by 2)

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4]; #cos(x)

y=[1.0,0.995,0.98007,0.95537,0.92106,0.87758,0.82534,0.76484,0.69671,0.62161,0.5403,0.4536,0.36235,0.2675,0.16997];

# Initialization

n=len(x); # number of loop

sum2O=0.0;

sum2E=0.0;

h=x[2]-x[1]

###############################

#

# Main body of Program

#

for i in range(1,n-1) : # number element statrs i=0 and ends i=n-1

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sum2E=sum2E+y[i]

else : # odd for 2 divisable condition,

sum2O=sum2O+y[i]

# end of for Loop

#

intg=(h/3)*(y[0]+4*sum2O+2*sum2E+y[n-1]); # 1st element 1=0 and last element i=n-1

#

###### OutPut ################

#

print("the value of integration is ", intg)

#

# Simpson 1/3 rule for Integration

# Note that this rule is effective for n is

# devisable by 2, i.e., y(0) to y(n:devisable by 2)

#

############ Input ############

# Function

def f(x) :

import math

return math.log(x)

x0=4

xn=5.2

h=0.0001

n=int((xn-x0)/h)

x=[x0]

y=[f(x0)]

##### Creating aray of data

for i in range(1,n+1) :# We need run the loop 1 to n, i takes value form 1 to n

xi=x0+i*h

yi=f(xi)

x.append(xi)

y.append(yi)

# Initialization

sum2O=0.0;

sum2E=0.0;

###############################

#

# Main body of Program

#

for i in range(1,n) : # number element statrs i=0 and ends i=n

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sum2E=sum2E+y[i]

else : # odd for 2 divisable condition,

sum2O=sum2O+y[i]

# end of for Loop

#

intg=(h/3)*(y[0]+4*sum2O+2*sum2E+y[n]); # 1st element 1=0 and last element i=n

#

###### OutPut ################

#

print("the value of integration is ", intg)

#

# Simpson 3/8 rule for Integration

# Note that this rule is effective for n is

# divisable by 3, i.e., y(0) to y(n:devisable by 3)

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5]; #cos(x)

y=[1.0,0.995,0.98007,0.95537,0.92106,0.87758,0.82534,0.76484,0.69671,0.62161,0.5403,0.4536,0.36235,0.2675,0.16997,0.070737];

# Initialization

n=len(x); # number of loop

sum3O=0.0;

sum3E=0.0;

h=x[2]-x[1]

###############################

#

# Main body of Program

#

for i in range(1,n-1) : # number element statrs i=0 and ends i=n-1

if i%3==0 : # even for 3 divisable condition, % sign for modulus (vagses debe)

sum3E=sum3E+y[i]

else : # odd for 3 divisable condition,

sum3O=sum3O+y[i]

# end of for Loop

#

intg=(3*h/8)*(y[0]+3*sum3O+2*sum3E+y[n-1]); # 1st element 1=0 and last element i=n-1

#

###### OutPut ################

#

print("the value of integration is ", intg)

#

# Simpson 1/3 rule for Integration

# Note that this rule is effective for n is

# devisable by 2, i.e., y(0) to y(n:devisable by 2)

#

# Simpson 3/8 rule for Integration

# Note that this rule is effective for n is

# divisable by 3, i.e., y(0) to y(n:devisable by 3)

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

x=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7]; #cos(x)

y=[1.0,0.995,0.98007,0.95537,0.92106,0.87758,0.82534,0.76484,0.69671,0.62161,0.5403,0.4536,0.36235,0.2675,0.16997,0.070737,-0.0292,-0.128844];

# Initialization

n=len(x); # number of loop

sum3O=0.0;

sum3E=0.0;

sum2O=0.0;

sum2E=0.0;

h=x[2]-x[1]

###############################

# range devidation x=0.0 to 1.4 and x=1.4 to 1.7

# 1/3 rule applied in the range x=0.0 to x=1.4, n=0 to n=14

# which is devisable by 2

N3=14;

# 3/8 rule is applied in the range x=1.4 to 1.7, n=0 to n=3

# which is devisable by 3

#

###############################

# Main body of Program

###############################

# 1/3 rule

for i in range(1,N3) : # number element statrs i=0 and ends i=N3

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sum2E=sum2E+y[i]

else : # odd for 2 divisable condition,

sum2O=sum2O+y[i]

# end of for Loop

#

intg3=(h/3)*(y[0]+4*sum2O+2*sum2E+y[N3]);

# 1st element 1=0 and last element i=N3

#

###############################

# 3/8 rule

j=1;

for i in range(N3+1,n-1) : # number element statrs i=N3 and ends i=n-1

if j%3==0 : # even for 3 divisable condition, % sign for modulus (vagses debe)

sum3E=sum3E+y[i]

j=j+1

else : # odd for 3 divisable condition,

sum3O=sum3O+y[i]

j=j+1

# end of for Loop

#

intg8=(3*h/8)*(y[N3]+3*sum3O+2*sum3E+y[n-1]);

# 1st element i=N3 and last element i=n-1

#

###### OutPut ################

#

print("the value of integration is ", intg3+intg8)

#

#Python Code for Secant Method

# to find root of f(y)=y**3/3-2*y**2+y-4

#

def f(y):

return y**3/3-2*y**2+y-4

#

# INPUT

#def secant(f):

x0= 2 #put the initial starting point x0

xm1= 2.1 #put the second starting point x1

er=10**(-10) #error bar

n=1;

x=[xm1,x0];

#

# CALCULATION FOR ROOT OF F(Y)=0

#

while abs(f(x[n]))>10**(-er):

n=n+1;

x.append((x[n-1]*f(x[n-2])-x[n-2]*f(x[n-1]))/(f(x[n-2])-f(x[n-1])));

#

# OUTPUT

#

print("solution for f(y)=0 : y=",x[n])

# Python Program for Solution of Ordinary Differtial Equantion

# Runge-Kutta 2nd order

#

# Function

def f(x,y):

import math

return (2*math.sin(3*x) - x**2*y**2)/math.exp(y)

#

# Initialization

y=[5.0]

x=[0.0]

h=0.005 # step size

xn = 0.5

n = round(xn/h) # total number (interger) of step requird to reach xn

#

# Main body of Program

#

i=1

while i <= n :

xim=x[i-1]

yim=y[i-1]

xi=xim+h

k1=f(xim,yim)

k2=f(xim+h, yim + k1*h)

yi=yim + (k1/2 + k2/2)*h # y obtained by Runge-Kutta 2nd order method

x.append(xi) # append add xi in the list of x

y.append(yi)

i=i+1

#

# Output

#

print('The value of function at x=',x[n],'is ',y[n])

print('Actual value of y at x=',x[n],'is',4.99715389095)

#

# Python Program for Solution of Ordinary Differtial Equantion

# Runge-Kutta 4th order

#

# Function

def f(x,y):

import math

return (2*math.sin(3*x) - x**2*y**2)/math.exp(y)

#

# Initialization

y=[5.0]

x=[0.0]

h=0.005 # step size

xn = 0.5

n = round(xn/h) # total number (interger) of step requird to reach xn

#

# Main body of Program

#

i=1

while i <= n :

xim=x[i-1]

yim=y[i-1]

xi=xim+h

#

k1=f(xim,yim)

k2=f(xim + h/2, yim + k1*h/2)

k3=f(xim + h/2, yim + k2*h/2)

k4=f(xim + h, yim + k3*h)

#

yi=yim + (k1 + 2*k2 + 2*k3 + k4)*h/6 # y obtained by Runge-Kutta 4th order method

x.append(xi) # append add xi in the list of x

y.append(yi)

i=i+1

#

# Output

#

print('The value of function at x=',x[n],'is ',y[n])

print('Actual value of y at x=',x[n],'is',4.99715389095)

#

# Python Program for Solution of Ordinary Differtial Equantion

# RC Circuit (Charging) Runge-Kutta 4th order

#

# RC circuit charging equation : R Q'[t] + Q[t]/C = V0

# Values : C=5*10**(-6) F, R=5 Ohm, V0=10 volt, Q[t=0]=0

#

# Analytic Solution : Q[t] = C V0(1- Exp[-t/(RC)])

#

# Function

def qp(t,q):

# q[t=0]=0.0

# C=5*10**(-3)

# R=5.0

# V0=10.

return 10/5 - q/(25*10**(-6)) #Q'[t] = V0/R - Q[t]/(C*R)

#

# Analytic Solution

def Aq(t) :

import math

return 50*10**(-6)*(1-math.exp(-t/(25*10**(-6)))) #C V0 (1-math.exp(-t/(C*R))

#

# Initialization

t=[0.0] # time aray for numerical solution

tA=[0.0] # time aray for analytic solution

q=[0]

qA=[0]

#C=5*10**(-6)

#R=5.0

h= 10**(-6) # step size : 1 micro second

tn = 2*10**(-4) # we are looking 0 to 100 micro second

n = round(tn/h) # total number (interger) of step requird to reach xn

hA=10**(-8) # Step size for Analytic Solution : 10 ns

nA=round(tn/hA) # Number of data for analytic solution

#

# Main body of Program

#

i=1

while i <= n :

tim=t[i-1]

qim=q[i-1]

ti=tim+h

#

qi=qim + qp(tim,qim)*h # y obtained by Runge-Kutta 4th order method

#

t.append(ti) # append add xi in the list of x

q.append(qi)

i=i+1

j=1

while j<=nA :

tAj=tA[0] + j*hA

qAj=Aq(tAj) # correpsonding analytic values

tA.append(tAj)

qA.append(qAj)

j=j+1

#

# Output plot

#

tdata=[t[0],t[5],t[10],t[15],t[20],t[30],t[40],t[50],t[60],t[80],t[100],t[120],t[140],t[160],t[180],t[200]] # some particular numerical data

vdata=[q[0],q[5],q[10],q[15],q[20],q[30],q[40],q[50],q[60],q[80],q[100],q[120],q[140],q[160],q[180],q[200]]

#import matplotlib

#from matplotlitb import pyplot as plt

#

import matplotlib.pyplot as plt

#import numpy as np # numpy is used for array

#X = np.linspace(0, tn, 256, endpoint=True) # one dimentional array for t

#AP1=np.Av(X) # creating one dimentional array for v(t)

#

plt.figure(figsize=(8,5))

plt.plot(tA, qA, color="green", linewidth=1.0, linestyle="-")

plt.plot(t, q,color="red",linewidth=1.0, linestyle="--")

plt.plot(tdata, vdata,"ro")

plt.axis([0, tn, 0, q[0]])

plt.ylim(0, 0.00006) # Set y-axis limit

plt.ylabel('Charge')

plt.xlabel('time')

plt.title('Plot for Variation of Load Volage over Time (RC-Charging)')

plt.show()

#

# Python Program for Solution of Ordinary Differtial Equantion

# RC Circuit (Discharging) Runge-Kutta 4th order

#

# RC circuit equation : C V'[t] + V[t]/R = 0

# Values : C=5*10**(-6) F, R=5 Ohm, V[t=0]=3 V

#

# Analytic Solution : V[t] = V[t=0] Exp[-t/(RC)]

#

# Function

def vp(t,v):

# v0=3.0

# C=5*10**(-6)

# R=5.0

return -v/(25*10**(-6)) #v0/(C*R)

#

# Analytic Solution

def Av(t) :

import math

# v0=3.0

# C=5*10**(-6)

# R=5.0

return 3*math.exp(-t/(25*10**(-6))) #v0*math.exp(-t/(C*R)

#

# Initialization

t=[0.0] # time aray for numerical solution

tA=[0.0] # time aray for analytic solution

v=[3.0]

vA=[3.0]

#C=5*10**(-6)

#R=5.0

h= 10**(-6) # step size : 1 micro second

tn = 10**(-4) # we are looking 0 to 100 micro second

n = round(tn/h) # total number (interger) of step requird to reach xn

hA=10**(-8) # Step size for Analytic Solution : 10 ns

nA=round(tn/hA) # Number of data for analytic solution

#

# Main body of Program

#

i=1

while i <= n :

tim=t[i-1]

vim=v[i-1]

ti=tim+h

#

k1=vp(tim,vim)

k2=vp(tim + h/2, vim + k1*h/2)

k3=vp(tim + h/2, vim + k2*h/2)

k4=vp(tim + h, vim + k3*h)

#

vi=vim + (k1 + 2*k2 + 2*k3 + k4)*h/6 # y obtained by Runge-Kutta 4th order method

t.append(ti) # append add xi in the list of x

v.append(vi)

i=i+1

j=1

while j<=nA :

tAj=tA[0] + j*hA

vAj=Av(tAj) # correpsonding analytic values

tA.append(tAj)

vA.append(vAj)

j=j+1

#

# Output plot

#

tdata=[t[0],t[20],t[40],t[50],t[60],t[80],t[100]] # some particular numerical data

vdata=[v[0],v[20],v[40],v[50],v[60],v[80],v[100]]

#import matplotlib

#from matplotlitb import pyplot as plt

#

import matplotlib.pyplot as plt

#import numpy as np # numpy is used for array

#X = np.linspace(0, tn, 256, endpoint=True) # one dimentional array for t

#AP1=np.Av(X) # creating one dimentional array for v(t)

#

plt.figure(figsize=(8,5))

plt.plot(tA, vA, color="green", linewidth=1.0, linestyle="-")

plt.plot(t, v,color="red",linewidth=1.0, linestyle="--")

plt.plot(tdata, vdata,"ro")

plt.axis([-0.0000005, tn+0.0000009, 0, v[0]+0.1])

plt.ylabel('Voltage')

plt.xlabel('time')

plt.title('Plot for Variation of Load Volage over Time (RC-Discharging)')

plt.show()

#

# Python Program for Solution of Ordinary Differtial Equantion

# Radioactive Decay : Runge-Kutta 4th order

#

# Radioactive deacy equation : dN/dt = - N/T

# Values : T decay constant = 1

#

# Analytic Solution : N[t] = N[t=0] Exp[-t/T]

#

# Function

def NM(t,N):

T=1

return -N/T

#

# Analytic Solution

def ANM(t) :

import math

T=1

N0=100

return N0*math.exp(-t/T) #N0*math.exp(-t/T)

#

# Initialization

t=[0.0] # time aray for numerical solution

tA=[0.0] # time aray for analytic solution

N=[100.0]

NA=[100.0]

#C=5*10**(-6)

#R=5.0

h= 10**(-2) # step size for analytical solution

tn = 5 # we are looking 0 to 5 second

m = round(tn/h) # total number (interger) of step requird to reach xn

hA=10**(-5) # Step size for Analytic Solution

mA=round(tn/hA) # Number of data for analytic solution

#

# Main body of Program

#

i=1

while i <= m :

tim=t[i-1]

Nim=N[i-1]

ti=tim+h

#

k1=NM(tim,Nim)

k2=NM(tim + h/2, Nim + k1*h/2)

k3=NM(tim + h/2, Nim + k2*h/2)

k4=NM(tim + h, Nim + k3*h)

#

Ni=Nim + (k1 + 2*k2 + 2*k3 + k4)*h/6 # y obtained by Runge-Kutta 4th order method

t.append(ti) # append add xi in the list of x

N.append(Ni)

i=i+1

j=1

while j<=mA :

tAj=tA[0] + j*hA

NAj=ANM(tAj) # correpsonding analytic values

tA.append(tAj)

NA.append(NAj)

j=j+1

#

# Output plot

#

tdata=[t[0],t[50],t[100],t[200],t[300],t[400],t[500]] # some particular numerical data

Ndata=[N[0],N[50],N[100],N[200],N[300],N[400],N[500]]

#import matplotlib

#from matplotlitb import pyplot as plt

#

import matplotlib.pyplot as plt

#import numpy as np # numpy is used for array

#X = np.linspace(0, tn, 256, endpoint=True) # one dimentional array for t

#AP1=np.Av(X) # creating one dimentional array for v(t)

#

plt.plot(tA, NA, color="green", linewidth=1.0, linestyle="-")

plt.plot(t, N,color="red",linewidth=1.0, linestyle="--")

plt.plot(tdata, Ndata,"ro")

plt.axis([0, tn, 0, N[0]+3])

plt.ylabel('N(t)')

plt.xlabel('time')

plt.title('Plot for Radioactive Decay with decay constant =1')

plt.show()

#Python code for Newton method

#our function f(y) = y3/3 − 2y2 + y − 4

#

#

def f(y):

return y**3/3-2*y**2+y-4

# derivative of the above function

def df(y):

return y**2-4*y+1

#

#

#def newton(f,df):

x0=float(input("entre the initial value for the solution with decimal point\n"))

er=int(input("entre value of N for error bar 10**(-N), say, 7\n"))

n=0;

x=[x0];

while abs(f(x[n]))>10**(-er):

n=n+1;

x.append(x[n-1]-f(x[n-1])/df(x[n-1]));

print("solution for f(y)=0 : y=",x[n])

# Sum of odd infinte series

##############################

sumO=0.0

for i in range(0,2000000) :

sumO=sumO + (-1)**i*(1/(1+2*i))

print('value of pi=',4*sumO)

#

# Python Program for Solution of Ordinary Differtial Equantion

# Modified Eular's method

#

# Newton Law of Motion : Simple Harmonic Motion

#

# F = m (d^2 x)/dt^2 = -k x, can be written as 1st order equation

# dv/dt = - x , k=1, m=1

# dx/dt = v

#

# Analytical Soln : x = cos(t) and v = - sin(t),

#

# Function

def fv(x,v,t):

return -x

#

def fx(x,v, t) :

return v

# Initialization

v=[0.0] # intial velocity zero

x=[1.0] # initial position 0.25 meter

t=[0.0] # initial time 0

h=0.0001 # step size

tn = 5 # final time

n = round(tn/h) # total number (interger) of step requird to reach xn

#

# Main body of Program

#

i=1

while i <= n :

tim=t[i-1]

vim=v[i-1]

xim=x[i-1]

#

ti=tim+h

vi = vim + fv(xim,vim,tim)*h # v obtained by Eular's method

xi= xim + fx(xim,vim,tim)*h

#

vmi=vim + (fv(xim,vim,tim)+fv(xi,vi,ti))*h/2 # v obtained by modified Eular's method

xmi=xim + (fx(xim,vim,tim)+fx(xi,vi,ti))*h/2

#

t.append(ti) # append add ti in the list of t

v.append(vmi)

x.append(xmi)

i=i+1

#

# Output

#

print('at time t ', tn, 'velocity v=',v[n],'and position x= ',x[n])

# Analytical Soln : x = 0.1 cos(w t) and v = 0.1 w sin(w t), w = sqrt(k/m)

import math

print('Actual value of v=',-math.sin(tn),'x',math.cos(tn),)

#

# Newton Forward Difference Method for Interpolation

#

############ Input ############

#x=[0,0.001,0.002,0.003,0.004,0.005];

#y=[1.121,1.123,1.1255,1.127,1.128,1.1285];

x=[0,0.05, 0.1,0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5];

y=[0.0, 0.05004, 0.10034, 0.15114, 0.20271, 0.25534, 0.30934, 0.36503, 0.42279, 0.48305, 0.5463];

z0=0.07; # find y at z0

# Initialization

n=len(x); # number of loop

summ=y[0];

coef=1;

u=(z0-x[0])/(x[1]-x[0]);

###############################

#

# Main body of Program

#

for k in range(1,n) :

i=1

while i<n :

y[i-1]=y[i]-y[i-1]; #forward difference

i=i+1

# End of inner while loop

coef=coef*(u-k+1)/k # coefficient

summ=summ+coef*y[0]

# end of outer for loop

#

###### OutPut ################

#

print("the value of function at ",z0," is ", summ)

#

# Python Program for Solution of Ordinary Differtial Equantion

# Newton's Law of Cooling : Runge-Kutta 4th order

#

# Suppose

#

# Cooling dy/dt = -(1/27)*(y-65) and initial condition y(t=0)=200

#

# Analytic Solution : y[t] = 65 + 135*Exp[-t/27]

# From this solution one can find value of A from the consition y(1)=35

#

# Function

def yt(t,y):

return -(1/27)*(y-65)

#

# Analytic Solution

def Ay(t) :

import math

return 65+135*math.exp(-t/27)

#

# Initialization

t=[0.0] # time aray for numerical solution

tA=[0.0] # time aray for analytic solution

yd=[200.0] # initial condition

yA=[200.0] # initial condition

#C=5*10**(-6)

#R=5.0

h = 1 # step size : 1 micro second

tn = 10 # we are looking 0 to 3 hour

n = round(tn/h) # total number (interger) of step requird to reach xn

hA=10**(-5) # Step size for Analytic Solution : 10 ns

nA=round(tn/hA) # Number of data for analytic solution

#

# Main body of Program

#

i=1

while i <= n :

tim=t[i-1]

ydim=yd[i-1]

ti=tim+h

#

k1=yt(tim,ydim)

k2=yt(tim + h/2, ydim + k1*h/2)

k3=yt(tim + h/2, ydim + k2*h/2)

k4=yt(tim + h, ydim + k3*h)

#

ydi=ydim + (k1 + 2*k2 + 2*k3 + k4)*h/6 # y obtained by Runge-Kutta 4th order method

t.append(ti) # append add xi in the list of x

yd.append(ydi)

i=i+1

j=1

while j<=nA :

tAj=tA[0] + j*hA

yAj=Ay(tAj) # correpsonding analytic values

tA.append(tAj)

yA.append(yAj)

j=j+200

#

# Output plot

#

tdata=[t[0],t[1],t[2],t[3],t[4],t[5],t[6],t[7],t[8],t[9],t[10]] # some particular numerical data

ydata=[yd[0],yd[1],yd[2],yd[3],yd[4],yd[5],yd[6],yd[7],yd[8],yd[9],yd[10]]

#

import matplotlib.pyplot as plt

import numpy as np # numpy is used for array

X = np.linspace(0, tn, 256, endpoint=True) # one dimentional array for t

#

plt.plot(tA, yA, color="green", linewidth=1.0, linestyle="-")

#plt.plot(t, yd,color="red",linewidth=1.0, linestyle="--")

plt.plot(tdata, ydata,"ro")

plt.ylabel('Voltage')

plt.xlabel('time')

plt.show()

print('The value of function at x=',t[n],'is ',yd[n],'Actual :', Ay(t[n]))

#print('Actual value of y at x=',x[tn],'is',4.99715389095)

#

#Python Code for Secant Method

# to find root of f(y)=y**3/3-2*y**2+y-4

#

def f(y):

return y**3/3-2*y**2+y-4

#

# INPUT

#def secant(f):

x0= 2 #put the initial starting point x0

xm1= 2.1 #put the second starting point x1

er=10**(-10) #error bar

n=1;

x=[xm1,x0];

#

# CALCULATION FOR ROOT OF F(Y)=0

#

while abs(f(x[n]))>10**(-er):

n=n+1;

x.append((x[n-1]*f(x[n-2])-x[n-2]*f(x[n-1]))/(f(x[n-2])-f(x[n-1])));

#

# OUTPUT

#

print("solution for f(y)=0 : y=",x[n])

#

# Monte Carlo Methods : Hit-or-Miss Method

#

# Integration of sin^2(1/x) over -pi to pi

#

# calling mathematical funtion

#

import math

#

# Function deffination

def f(x) :

return (math.sin(1/x))**2

#

# Range of integration

#

a=-math.pi

b=math.pi

#

# total number of sample point

N=10**(5);

Nfx=0; # number of sample point lies with in the curve

i=0; # loop continuation

#

# Main body of program

#

while i <=N :

#

# generation of random variable for f(x) = y and and x=z

#

import random #calling random function

x=a + random.random()*(b-a); # creat floating random number within [0.0,1.0]

y=random.random()

#

if y <= f(x) :

Nfx=Nfx+1 # number of point within the curve

i=i+1

# Area of sample space

A= (b-a)*1 # here 1 is the maximum value

#

# Output

#

print("the value of integration : ",A*Nfx/N)

#

# Python Program for Solution of Ordinary Differtial Equantion

# Modified Eular's method

#

# Function

def f(x,y):

import math

return (2*math.sin(3*x) - x**2*y**2)/math.exp(y)

#

# Initialization

y=[5.0]

x=[0.0]

h=0.005 # step size

xn = 0.5

n = round(xn/h) # total number (interger) of step requird to reach xn

#

# Main body of Program

#

i=1

while i <= n :

xim=x[i-1]

yim=y[i-1]

xi=xim+h

yi = yim + f(xim,yim)*h # y obtained by Eular's method

ymi=yim + (f(xim,yim)+f(xi,yi))*h/2 # y obtained by modified Eular's method

x.append(xi) # append add xi in the list of x

y.append(ymi)

i=i+1

#

# Output

#

print('The value of function at x=',x[n],'is ',y[n])

print('Actual value of y at x=',x[n],'is',4.99715389095)

########################################################

# Program for Least Square Fitting ###################

########################################################

#

#

# Input x-values and y-values

#

import numpy as np

x=np.array([2.23, 4.78, 7.21, 9.37, 11.64, 14.23, 16.55, 18.70, 21.05])

y=np.array([139, 123, 115, 96, 62, 54, 10, -3, -13])

#

# creating list of {x1*y1, x2*y2, x3*y3, x4*y4, ...., xn*yn}

n=len(x)

xy=np.zeros(n) # creating array of length n with zeros

x2=np.zeros(n)

for i in range(0,n) :

xy[i]=x[i]*y[i]

x2[i]=x[i]*x[i]

# Note that above operation can be done without loop as

# xy=x*y and x2=x*x

#

# Calculation of average of lists x, y and xy

xav = np.mean(x)

yav=np.mean(y)

xysum=sum(xy)

x2sum=sum(x2)

#

# Equantion of line : y = a + b x and corresponding data file

b = (xysum - n*xav*yav)/(x2sum - n*xav*xav)

a=yav - b*xav

xt=np.linspace(x[0],x[n-1],50)

yf= a + b*xt

#

# Output as Figure

import matplotlib.pyplot as plt

plt.plot(x, y, 'ro', label='data')

plt.plot(xt, yf, 'b-', label='fitted curve')

plt.xlabel('t, time')

plt.ylabel('v, velocity')

plt.legend(loc='upper right')

plt.show()

#

# Save figure in the directory path : /Users/Books/Python Prog

#plt.savefig("Path/leat2fit.png", dpi = 200)

"""

Program for Guass Seidel Method

"""

import numpy as np

#A=np.array([[20,15,10,45],[-3,-2.249,7,1.751],[5,1,3,9]],float)

# diagonally dominant form for matrix with cofficient of x1, x2, x3

A=np.array([[5,1,3,9],[20,15,10,45],[-3,-2.249,7,1.751]],float)

#A=np.array([[5,-2,3,-1],[-3,9,1,2],[2,-1,-7,3]],float)

#A=np.array([[4,2,-1,-1],[3,-5,1,3],[1,0,2,-4]],float)

# Initialization values of x1, x2, x3

xn=np.shape(A)[0] # total number of variable

x=np.zeros(xn,float)

# itaration number

l=int(input("entre the number of iteration, say 50 \n"))

for k in range(l+1):

x[0]=(A[0,3] - (A[0,2]*x[2] + A[0,1]*x[1]))/A[0,0]

x[1]=(A[1,3] - (A[1,2]*x[2] + A[1,0]*x[0]))/A[1,1]

x[2]=(A[2,3] - (A[2,1]*x[1] + A[2,0]*x[0]))/A[2,2]

print('the solutions are','x1=',x[0],',','x2=',x[1],',','x3=',x[2])

#print(A[0,0]*x[0] + A[0,1]*x[1] + A[0,2]*x[2],',',A[1,0]*x[0] + A[1,1]*x[1] + A[1,2]*x[2],',', A[2,0]*x[0] + A[2,1]*x[1] + A[2,2]*x[2])

# Fibonacci series

#

#

# Initialization

x=[0,1]

# number of series element

n=20

for i in range(2,n+1) :

xi=x[i-1]+x[i-2]

x.append(xi)

#

# Print

print("Fibonacci Series \n",x)

#Python code for Bisection method

#our function f(y) = x**2-2

#

#

def f(x):

return x**2-2

#

# INPUT

#

a=0 # 1ST POINT

b=2 # 2ND POINT

n=30 #number of loop

er=10**(-7) #error bar

#

# CALCULATION FOR FINDING ROOT OF F(Y)=0

#

for i in range (n) : # here rang(N) defines how many times loop will continue, here 50 times

c=(a+b)/2

if abs(f(c))<er : # condition for error bar of solution (1ST IF LOOP)

print("found zero at x=",c)

break #break terminate the progam when the above conditon is satisfied

# first if ends here

if f(a)*f(c)<0 : #2nd if loop

a=a

b=c

else: # this else for 2nd if loop

a=c

b=b # 2nd if ends here

else : # this else for "for" loop

# OUTPUT

print("solution with the given error bar is not found")

# for loop ends here

#

# Trapizoidal rule for Integration

#

############ Input ############

#x=[1.0,1.2,1.4,1.6,1.8,2.0,2.2];

#y=[2.7183,3.3201,4.0552,4.9530,6.0496,7.3891,9.0250];

# Data for upper curve x-> H and y->B function : 2 Tanh[(HU + 500)/800]

HU = [-3000, -2700, -2400, -2100, -1800, -1500, -1200, -900, -600, -300, 0.0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]; #H1 (x)

BU = [-1.9923, -1.9837, -1.9657, -1.9281, -1.8507, -1.6966, -1.4078, -0.9242, -0.2487, 0.4898, 1.1092, 1.5232, 1.7597, 1.8828, 1.9438, 1.9732, 1.9873, 1.994, 1.9972, 1.9987, 1.9994]; # B1 (y)

# Data for Lower curve x-> H and y->B function : 2 Tanh[(HD - 500)/800]

HD = [-3000, -2700, -2400, -2100, -1800, -1500, -1200, -900, -600, -300, 0.0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]

BD = [-1.9994, -1.9986, -1.9972, -1.994, -1.9873, -1.9732, -1.9437, -1.8828, -1.7597, -1.5232, -1.1092, -0.4898, 0.2487, 0.9242, 1.4078, 1.6966, 1.8507, 1.9281, 1.9657, 1.9837, 1.9923] # B2

# Initialization

nU=len(HU); # number of Upper-Loop

nD=len(HD); # number of Lower-Loop

#print(nU, nD)

###############################

summU=0.0;

summD=0.0;

###############################

hU=HU[2]-HU[1]

hD=HD[2]-HD[1]

###############################

#

# Main body of Program

#

# Integration on upper curve

#

for i in range(1,nU-1) :

summU = summU + BU[i]

# end of for Loop

#

intgU=hU*(BU[0]+2*summU+BU[nU-1])/2;

#

#

# Integration on Lower curve

#

for i in range(1,nD-1) :

summD = summD + BD[i]

# end of for Loop

#

intgD=hD*(HD[0]+2*summD+HD[nD-1])/2;

#

###### OutPut ################

#

print("The area of BH-Loop ", intgU-intgD)

print("Theoretical Area of BH-Loop based on cosidered function", 3994.34)

#

# Simpson 1/3 rule for Integration

# Note that this rule is effective for n is

# devisable by 2, i.e., y(0) to y(n:devisable by 2)

#

############ Input ############

# Data for upper curve x-> H and y->B function : 2 Tanh[(HU + 500)/800]

HU = [-3000, -2700, -2400, -2100, -1800, -1500, -1200, -900, -600, -300, 0.0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]; #H1 (x)

BU = [-1.9923, -1.9837, -1.9657, -1.9281, -1.8507, -1.6966, -1.4078, -0.9242, -0.2487, 0.4898, 1.1092, 1.5232, 1.7597, 1.8828, 1.9438, 1.9732, 1.9873, 1.994, 1.9972, 1.9987, 1.9994]; # B1 (y)

# Data for Lower curve x-> H and y->B function : 2 Tanh[(HD - 500)/800]

HD = [-3000, -2700, -2400, -2100, -1800, -1500, -1200, -900, -600, -300, 0.0, 300, 600, 900, 1200, 1500, 1800, 2100, 2400, 2700, 3000]

BD = [-1.9994, -1.9986, -1.9972, -1.994, -1.9873, -1.9732, -1.9437, -1.8828, -1.7597, -1.5232, -1.1092, -0.4898, 0.2487, 0.9242, 1.4078, 1.6966, 1.8507, 1.9281, 1.9657, 1.9837, 1.9923] # B2

# Initialization

nU=len(HU); # number of Upper-Loop

nD=len(HD); # number of Lower-Loop

#print(nU, nD)

###############################

sumUE=0.0;

sumUO=0.0;

sumDE=0.0;

sumDO=0.0;

###############################

hU=HU[2]-HU[1]

hD=HD[2]-HD[1]

###############################

#

# Main body of Program

#

# Integration on upper curve

#

for i in range(1,nU-1) : # number element statrs i=0 and ends i=n-1

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sumUE=sumUE+BU[i]

else : # odd for 2 divisable condition,

sumUO=sumUO+BU[i]

# end of for Loop

#

intgU=(hU/3)*(BU[0]+4*sumUO+2*sumUE+BU[nU-1]); # 1st element 1=0 and last element i=n-1

#

#

# Integration on Lower curve

#

for i in range(1,nD-1) : # number element statrs i=0 and ends i=n-1

if i%2==0 : # even for 2 divisable condition, % sign for modulus (vagses debe)

sumDE=sumDE+BD[i]

else : # odd for 2 divisable condition,

sumDO=sumDO+BD[i]

# end of for Loop

#

intgD=(hD/3)*(BD[0]+4*sumDO+2*sumDE+BD[nD-1]); # 1st element 1=0 and last element i=n-1

#

###### OutPut ################

#

print("The area of BH-lopp ", intgU-intgD)

print("Theoretical Area of BH-Loop based on cosidered function", 3994.34)