Basic exercises of calculations and plotting.
Exercise 4
Exercise 9
Reads datacolumns YYYY-MM-DD XXXXXX from a textfile kwh.dat
Class for testing
Must be in the same directory
Calculates the integrals using trapezoidals
from scipy import *
from pylab import *
from math import *
import sys
import matplotlib.pyplot as plt
#task 6
L9=list(range(101)) #skapar lista L9 med 100 värden mellan 0-100
L10 = [k/100 for k in L9] #skapar lsita L10 där varje värde från L9 / 100
print(L10)
s = 0 #task 7
for i in range(0, 500): #räkne-loop går igenom i 0-500 (slutar vid 500)
s = s + i #adderar i's stegrande värde med alla föregående
print(s)
ss = [0] #task 7
for i in range(1, 500): #gör som ovanstående
ss.append(ss[i-1] + i) #men sparar värje iterering i lista
print(ss)
xplot=list(range(100)) #task 8
for i in range(0, 100): #räknar variabel i mellan 0-100
xplot[i] = i/100 #lägger in variabel i/100 i listan xplot
print(xplot)
yplot=list(range(100)) #task 9
for i in range(0, 100): #räknar variabel i mellan 0-100
yplot[i]=math.atan(i/100) #skapar arctan-värden mellan 0-1
print(yplot)
plot (xplot, yplot) #task 10, plottar x, y värden i listorna
print(plot)
p = 0.0
L=list(range(200)) #task 11
for i in range(1, 200):
L[i] = 1/math.sqrt(i)
p = p + L[i]
print(p)
h = 1/1000 #task 12 extra
a = -0.5
n = 0
u=list(range(1000))
u[0]=math.exp(0)
u[1]=math.exp(h*a)
u[2]=math.exp(2*h*a)
for n in range(0, 997):
u[n+3] = u[n+2] + h*a*((23/12)*u[n+2] - (4/3)*u[n+1] + (5/12)*u[n])
td=list(range(1000))
for n in range(0, 1000):
td[n] = h*n
plot (u, td)
print(plot)
td2=list(range(1000))
for n in range(0, 1000):
td2[n] = abs(math.exp(a*td[n])-u[n])
plot (u, td2)
plt.xlabel("u")
plt.ylabel("td2")
print(plot)Exercise 2from scipy import *
from pylab import *
from math import *
import sys
import matplotlib.pyplot as plt
x=0.5
a=8
c=0.
for i in range(201): # Task 1
c=x #saves x-value for previous iteration to c
x=math.sin(x)-a*x+30
if abs(c-x) <= 1.e-8: #if the difference between previous and new x-value
print('The result after {num} iterations is {res}'.format(num=i,res=x))
break
if i == 200:
print('The condition was not met within 200 iteration steps.')
break
x=0.5 #task 2
a=0.5
xplot=linspace(5, 30, num=26) #x-värden
yplot= []
#yplot=list(range(26))
for i in range(0, 26):
yplot.append(math.sin(xplot[i])-a*xplot[i]+30) #y-värden
plot (xplot, yplot) #plottar x, y värden i listorna
plot (xplot, xplot)
print(plot)
elements = []
y = 1.0
n = 0 # Task 3
while abs(y) >= 1.e-9: # Loopa så länge elementen är större än 1E-9
n = n + 1 # Stega upp 1 steg för varje iteration
y = sin(n)**2/n # N:te elementets värde
elements.append(y) # Elementet lagras i array elements
result = 'The loop had to do {} iteration steps .'
print(result.format(n))
negative = False
a = 0.25
c = 0.0
x = 1.0
i = 0 # Task 4
for i in range(1000):
c = x
x = 0.2*x-a*(x**2-5)
if x < 0:
negative = True
if abs(c-x) <= 1.e-9:
print('Sequence converged to x={num} after {res} iterations.'.format(num=x,res=i))
break
if i == 999:
print('No convergence detected.')
break
if negative == True:
print('There were negative elements in the sequence.')
a = 0.5
c = 0.0
x = 1.0
i = 0 # Task 5
minuslist = []
pluslist = []
for i in range(1000):
c = x
x = 0.2*x-a*(x**2-5)
if x < 0:
minuslist.append(x)
if x > 0:
pluslist.append(x)
if abs(c-x) <= 1.e-9:
print('Sequence converged to x={num} after {res} iterations.'.format(num=x,res=i))
break
if i == 999:
print('No convergence detected.')
break
def my_func(a):
c = 0.0
x = 1.0
i = 0 # Task 6
for i in range(31):
c = x
x = 0.2*x-a*(x**2-5)
if abs(c-x) <= 1.e-9:
return True
if i == 30:
return False
#if my_func(1) == True
print('Returned True')
def function(a, x):
c = 0.0
neg = []
pos = []
i = 0 # Task 8
for i in range(31):
c = x
x = 0.2*x-a*(x**2-5)
if x < 0:
neg.append(x)
if x > 0:
pos.append(x)
if abs(c-x) <= 1.e-9:
return neg and pos and True
if i == 30:
return neg and pos and False
#if my_func(1) == True
#print('Returned True')Exercise 3from scipy import *
from pylab import *
from math import *
import sys
import matplotlib.pyplot as plt
def f(m): #Task 3
L = [n-m/2 for n in range(m)] #L[0] = -m/2, L[-1] = m/2 -1
return 1 + L[0] + L[-1] #1 + -m/2 + m/2 -1 = 0 Always!!
print(f(100)) #Calls function with argument 100
L0 = [0, 20, 30, 40] #Task 4
L1 = [20, 0, 50, 60]
L2 = [30, 50, 0, 70]
L3 = [40, 60, 70, 0]
distance = [L0, L1, L2, L3]
reddistance=[] #Defines reddistance as a list
for i in range(1, len(L0)): #Counts from 1 to number of columns
temp = distance[i] #For each cycle each row stored in temp
reddistance.append(temp[0:i]) #Saves only i elements from each row
temp1=[]
temp2=[]
reddistance2=[]
for i in range(1, len(L0)): #Nr of iterations from 1 to nr of coulombs
temp1 = distance[i] #Stores first row in temporary list 1
for n in range(0, i): #Sorts out wanted elements in temp list 1
temp2.append(temp1[n]) #Saves the rows wanted elements in temp 2
reddistance2.append(temp2)
temp2=[]
reddistance3=[]
l=1
while l <= len(L0) - 1: #While loop
reddistance3.append(distance[l][0:l]) #Slicing elements of rows in dist
l = l + 1
A = {'banana','apple','coconut', 'orange'} #Task 5
B = {'ananas', 'orange', 'apple', 'mango'}
C = {}
def task5(A, B):
C = A.union(B) - A.intersection(B) #Symmetric difference of sets
return C
print(task5(A, B))
print(A.symmetric_difference(B)) #Does same thing but not same order
A.intersection(B) #Leaves A unchanged but presents intersection
A.intersection_update(B) #Changes A to only contain intersection
def func(x): #Task 8 and 9
y = (3*x**2 - 5) #Defines function to valuate by bisec
#y = math.atan(x) #Defines function to valuate by bisec
return y
a = -1.5 #Interval startpoint
b = -0.4 #Interval endpoint
t = 1e-2 #Precision of solution interval
def bisec(a, b, t):
while abs(a-b) > t: #If interval is bigger than tolerance
if func(a)*func((a+b)/2) <= 0: #If signchange in first interval,
b = (a+b)/2 #b gets new position
elif func((a+b)/2)*func(b) <= 0: #If signchange in second interval,
a = (a+b)/2 #a gets new position
elif func(a) == 0: #The root is exactly on position a
return (a, a, a)
elif func(b) == 0: #The root is exactly on position b
return (b, b, b)
else:
print('No sign-change in interval endpoints.')
return (0, 0, 0)
return (a, b, (a+b)/2) #Third parameter is presumed root
print(bisec(a, b, t)) #For checkingExercise 4
from scipy import *
from pylab import *
from math import *
from numpy import *
import sys
import matplotlib.pyplot as plt
def ff(r,fi): #Task 1
return r*exp(1j*fi) #Define complex function
phi=linspace(0,2*pi, 100) #Creates array-list of values betw 0 and 2pi
rad=linspace(0.0, 1.0, 10)
z=ff(rad,phi) #Defines variable z as result of ff
plt.plot(z.real,z.imag)
for r in rad: #For every radie (r) gör en cirkel
z=ff(r,phi) #Calls complex function ff
plt.plot(z.real,z.imag, linewidth=2) #Plottar cirkel
def f(x): #Task 2
y = x**2 + 5*x - 10 #Arbitrary function
return y
def fp(x): #Derivative of function above
y = 2*x + 5
return y
x0 = pi/2. #Starting point
tol = 1e-10 #Tolerance level
def newton(f, fp, x0, tol):
x = x0
conv = False
for n in range(400):
x0 = x
x = x - f(x)/fp(x)
if abs(x-x0) < tol: #Criterion for convergence
conv = True
break
return (x, conv, n)
x,conv,n = newton(f, fp, x0, tol)
print(x0, f(x0), x, f(x), conv)Exercise 5from scipy import *
from scipy.optimize import fsolve
from pylab import *
from numpy import *
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
def f(x):
y = sin(omega*x)
return y
a = 0
b = pi/2.
omega = 2*pi
yomega = []
print(quad(f,a,b)) #Task 1, integral from a to b of function f
def inte(w): #Task 2
a=0 #Function we calculate its integral
b=pi/2 #Omega (w) depends on input
def f(x): #Calls on quad with its own defined f
return sin(w*x) #Only returns first value [0]
return quad(f,a,b)[0]
print(inte(2*pi)) #Omega constant 2pi same as Task 1
xomega = linspace(0,2*pi,1000) #X-values for plot
yomega = [inte(w) for w in xomega]
plt.plot(xomega, yomega) #Y-values recalling inte-function
def p(x): #Task 3
y = x**2 + x - 3
return y
x0= 5
root = fsolve(p, x0)[0]
print('The root of function p(x):', root)
x1=0 #Task 4, a range from 1 to 5
def q(a): #Takes varying a as input
def r(x): #Function with varying a to solve roots foor
return a*x**2 + x - 3
return fsolve(r,x1) #Return root with wanted a
A = linspace(1,5) #Array with all our x-values
y_axis = [q(a) for a in A] #Calculate y-values for the x-values from q(a)
plt.plot(A, y_axis) #PlotsExercise 6from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
M = array([[1.,2.],[3.,4.]])
def func(M): #Task1
if (M==M.T).all(): #If matrixes are symmetric
a=1
elif (M==-M.T).all(): #If matrixes are skew-symmetric
a=-1
else:
a = 0
return a
v=array([1,0])
u=array([0,1])
"""this function takes two vectors as input and returns true if they
are ortoghonal, else it returns False"""
def h(u,v): #Task2
a=v@u
if a == 0: #If vectors are othogonal =0
return True
else:
return False
w=array([2,5,6])
def g(w): #Task3
b = w/norm(w)
return b
def k(w): #Task3
m=0.
n = len(w)
for i in range (n):
m = (w[i])**2 + m
m = sqrt(m)
return w/m
l = pi/3 #Task 4
N = array([[cos(l),sin(l)],[-sin(l),cos(l)]])
print(dot(N,N.T))
print(dot(N,inv(N)))
K=eye(20) #Task 5, creates unit 20x20 matrix
K = 4*K #Gives 4 on diagonal values
i, j = np.indices(K.shape)
K[i == j-1] = 1 #Superdiagonal with 1
K[i == j+1] = 1 #Subdiagonal with 1
L=sl.eig(K) #L is eigenvalues of matrix K
M=eye(20) #Task 6, creates unit 20x20 matrix
M = 4*M #Gives 4 on diagonal values
i, j = np.indices(M.shape)
M[i == j-1] = 1 #Superdiagonal with 1
M[i == j+1] = -1 #Subdiagonal with -1
O=sl.eig(M) #O is eigenvalues of matrix M
#More complicated matrix for the eigenvalues?Exercise 7from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
U = array([1,2,3]) # Task 1
V = array([4.5,5,6]) # U, V = arbitrary vectors
def vecdot(U, V): # Takes two vectors as input
s = 0 # Summarasition begins with 0
if len(U) != len(V): # Raises error message, diffrent dimensions
print("Vectors not of same dimensions.")
return
for i in range(len(U)): # Nr of iterations = nr of dimensions
s = U[i]*V[i] + s # Multiplicates each elements
return s # Returns summation
W = array([[4.5,5,6],[4,2,0],[3,4,8]]) # Task 2, arbitrary matrix
def matvec(W, V): # Matrix W and vector V as input
T = zeros((len(V)))
if len(W) != len(V):
print("Matrix and vector not of same dimension.")
return
for i in range(len(W)):
T[i] = vecdot(W[i], V) # Use each row on matrix as vector
return T
X = array([[7,0,2.2],[9,5,7],[5,4,1]]) # Task 3, another arbitrary matrix
def matmat(W, X): # Matrix W and X as input
l,m = W.shape # l = rows, m = elements in each row
R = zeros((l,m)) # Creates empty matrix of size l x m
temp = zeros(l) # temporary storage of separate columns as vector
if shape(W) != shape(X):
print("The matrices is not of the same dimensions.")
return
for n in range(m): # Takes column as a vector times matrix
for i in range(l): # Conversion column to vector
temp[i] = X[i,n] # Storing each column as temporary vector
R[n] = matvec(W, temp) # Calling function matvec to multiply
i = 0 # Resets counter for inner loop
R = transpose(R)
return R
L = array([[4.5,0,0],[4,2,0],[3,4,8]]) # Lower triangular matrix
def triang(L, V): # Task 4
M = zeros((len(V))) # Creates empty vector
for i in range(len(V)):
k = 0
j = 0
for j in range(i+1): # Remember the hotel nights
k = L[i,j]*V[j] + k # skips the zero positions in triang-matrix
M[i] = k
return MExercise 8from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
A = diag(10*[10.]) + diag(9*[1.],-1) + diag(9*[1.],+1) + diag(6*[2.],-4) + diag(6*[2.],4)
def task1(matrix): # Task 1
L = np.zeros(shape(matrix))
for i in range ((shape(matrix))[0]):
for j in range (i):
L[i,j] = matrix[i,j]
R = np.zeros(shape(matrix))
for i in reversed(range(shape(matrix)[0])):
for j in reversed(range (i)):
R[j,i] = matrix[j,i]
D = np.zeros(shape(matrix))
for i in range ((shape(matrix))[0]):
D[i,i] = matrix[i,i]
return(L,R,D)
L,R,D = task1(A)
b = dot(A,range(10))
def gaussjacobi(A,b): # Task 2
L,R,D = task1(A)
x=ones(len(A))
for i in range(100): # Solves eq system A@x = b
x = np.linalg.solve(D,-(L+R)@x+b)
defect = sqrt(norm(A@x-b)**2/len(A)) # Task 5
return x, defect
gaussjacobi(A,b)
def task3(D): #Task 3
diagvec = np.zeros(len(D))
for i in range(len(D)):
diagvec[i] = D[i,i]
return diagvec
d = task3(D)
def diagdom(matrix): #Task 4
L,R,D = task1(matrix)
d = task3(D)
for row in range(len(matrix)):
rowsum = 0
for i in range(len(matrix)):
rowsum += L[i,row] + R[i,row]
if d[i] < rowsum:
raise Exception("Not strictly diagonally row dominant matrix.")
return False
return True
def systemmatrix(d, k=-1000): # Task 6
zero = np.zeros((2,2))
I = np.diag((1,1))
K = array([[k, 0.5],[0.5, k]])
D = array([[-d, 1.0],[1.0, -d]])
matrix = np.bmat([[zero, I], [K, D]])
return matrix
dval = np.arange(0, 100, 0.2) # Task 7
for d in dval:
z = sl.eig(systemmatrix(d))[0] #[0] since we only want eigenvalues return
plt.scatter(z.real, z.imag) # alternativ plt.plot(z.real,z.imag, ',')
plt.title('Eigenvalues of parameter depending system matrix')
plt.xlabel('eigenvalue real part')
plt.ylabel('eigenvalue imaginary part')
plt.grid(which="both", axis="both", ls="-")Exercise 9
Reads datacolumns YYYY-MM-DD XXXXXX from a textfile kwh.dat
from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
newfile=open('newfile.dat','w') # Task 1
newfile.write('1 st line \n')
newfile.write('2 st line \n')
newfile.close()
file=open('newfile.dat','r')
text1=file.read()
text2=file.read()
def testfile():
file=open('newfile.dat','r')
return file
text1=list(testfile()) # Task 2
text2=list(testfile())
def kwhdat(): # Task 3
file=open('kwh.dat','r') # öppnar för läsning av fil YYYY-MM-DD XXXXXX
return file # returnerar rad för rad
fildata = list(kwhdat()) # Läser in filen som en lista "fildata"
filarray = np.asarray(fildata) # Konverterar listan till en array
yearmonthday = []
kwh = []
for line in filarray: # Går igenom raderna i filen kwh.dat
temp1, temp2 = line.split(" ") # Sorterar ut datumet och kwh-värdet
temp2 = int(temp2) # kwh-stäng görst om till int
kwh.extend([temp2]) # kwh-integer läggs till i kwh-array
yearmonthday.append(temp1) # datumet läggs till i listan yearmonthday
yearmonthday.reverse() # Task 4
kwh.reverse()
kwh = np.asarray(kwh) # gör om kwh-lista till array
monthly = np.diff(kwh) # Task 5
del yearmonthday[0] # tar bort startdatumet
plt.title('Elförbrukning per månad')
plt.xlabel('Datum') # Task 6
plt.ylabel('KWh')
xval = np.arange(1, 155, 1)
plt.plot(xval, monthly, '-')
ma = max(monthly) # Task 7, högsta värdet ur monthly
for i in range(len(monthly)): # letar efter positionen för det värdet
if monthly[i] == ma: # om värdet är det högsta
maxi = i # spara positionen
mi = min(monthly) # Lägsta värdet ur monthly
for i in range(len(monthly)):
if monthly[i] == mi:
mini = i
difference = np.zeros(len(monthly)-1) # Task 8
for i in range(len(difference)): # letar efter positionen för det värdet
difference[i] = monthly[i+1] - monthly[i]
diffmax = max(difference) # Task 7, högsta värdet ur monthly
for i in range(len(difference)): # letar efter positionen för det värdet
if difference[i] == diffmax: # om värdet är det högsta
maxpos = i # spara positionen
diffmin = min(difference) # Lägsta värdet ur monthly
for i in range(len(difference)):
if difference[i] == diffmin:
minpos = i
average = np.mean(monthly) #Task 9
print("Lägsta förbrukning:", yearmonthday[mini], monthly[mini], "kWh.") # Task 10
print("Högsta förbrukning:", yearmonthday[maxi], monthly[maxi], "kWh.")
print("Högsta förbrukningsökning mellan", yearmonthday[maxpos],"och", yearmonthday[maxpos+1], difference[maxpos], "kWh.")
print("Högsta förbrukningsminskning mellan", yearmonthday[minpos],"och", yearmonthday[minpos+1], difference[minpos], "kWh.")
print("Genomsnittlig månadsförbrukning:", average, "kWh.")Exercise 10from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
class RationalNumber:
'''
This is a docstring for the class RationalNumber.
@Properties:
gcd:
int: returns greatest common divisor
shorten:
list: returns shortened numerator and denominator
factor_for_shortening:
int: returns factor used for highest shortening
can_be_shortened:
bool: returns True if fraction can be shortened, else False
convert2float:
float: returns calculated real number of the fraction
del_value:
exception: numerator and denominator cant deletes seperately
'''
def __init__(self, numerator, denominator):
if not isinstance(numerator, int):
raise TypeError('Numerator should be of type int.')
if not isinstance(denominator ,int):
raise TypeError('Denominator should be of type int.')
self._numerator = numerator
self._denominator = denominator
def gcd(self, a, b):
if b==0: # Task 2, computes the greatest common divisor
return a # performs an inplace modification of the object
else:
return self.gcd(b, a%b)
def shorten(self):
shnumerat = self._numerator # shortened numerator
shdenomi = self._denominator # shortened denominator
factor = self.gcd(shnumerat , shdenomi)
shnumerat = shnumerat // factor
shdenomi = shdenomi // factor
return shnumerat, shdenomi
def factor_for_shortening(self): # Task 3
shnumerat = self._numerator # shortened numerator
shdenomi = self._denominator # shortened denominator
factor = self.gcd(shnumerat , shdenomi)
return factor
def can_be_shortened(self): # Task 3
numshort, denshort = self.shorten()
if numshort != self._numerator and denshort != self._denominator:
return True
else:
return False
def convert2float(self):
return float(self._numerator)/float(self._denominator)
def __add__(self, other): # __add__ is + sign
p1, q1 = self._numerator, self._denominator
if isinstance(other, RationalNumber):
p2, q2 = other.numerator, other.denominator
elif isinstance(other, int):
p2, q2 = other, 1
else:
raise TypeError('Should be rational or of type int.')
return RationalNumber(p1*q2 + p2*q1, q1*q2)
def __repr__(self):
return '{}/{}'.format(self._numerator, self._denominator)
def __radd__(self, other):
return self + other # reverses the role of self and other
def fget_n(self):
return self._numerator
def fset_n(self, numerator):
self._numerator = numerator
return
def fget_d(self):
return self._denominator
def fset_d(self, denominator):
self._denominator = denominator
return
def del_value(self): # Task 5
raise Exception('Numerator and denominator belongs together.')
numerator = property(fget = fget_n, fset = fset_n , fdel = del_value)
denominator = property(fget = fget_d, fset = fset_d , fdel = del_value)
q=RationalNumber(10,20) # command q.__init__(10,20) is executed
"""
q.numerator # returns 10
q.denominator # returns 20
q.convert2float() # returns 0.5
RationalNumber.convert2float(q) # same as RationalNumber.convert2float(q)
q.shorten() # returns (1, 2)
q.can_be_shortened() # returns True
q.factor_for_shortening() # returns 10, (10/20)/10 = 1/2
q.numerator = 5 # sets q.numerator to value 5
q.numerator # returns 5
q.factor_for_shortening() # returns 5, (5/20)/5 = 1/4
q.shorten() # returns (1, 4)
"""Exercise 11from scipy import *
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
import scipy.linalg as sl
import scipy.sparse as sp
class PolyGon:
'''
Takes corner point coordinates as input and computes the edges.
Input is a list of arrays with shape (2,0) for a points x,y-coordinate.
Parent of subclasses Rectangle and SpecialRectangle.
'''
def __init__(self, corners):
self.corners = corners
n = len(corners)
x = np.zeros(n+1)
y = np.zeros(n+1)
x[-1] = (self.corners[0])[0]
y[-1] = (self.corners[0])[1]
for i in range(n):
x[i] = (self.corners[i])[0]
y[i] = (self.corners[i])[1]
edges = np.zeros((n,2))
for i in range(n):
edges[i] = (x[i+1] - x[i]), (y[i+1] - y[i])
self._x, self._y = x, y
self.edges = edges
def plot(self):
plt.plot(self._x, self._y)
pass
class Rectangle(PolyGon):
'''
Inherits plot and edges from parent PolyGon.
'''
def area(self):
width = (max(self._x)-min(self._x))
length = (max(self._y)-min(self._y))
return width*length
class SpecialRectangle(Rectangle):
'''
Inherits area from parent Rectangle. Plot and edges is inherited from
grandparent PolyGon.
'''
def __contains__(self, other):
if max(other._x) <= max(self._x) and min(other._x) >= min(self._x):
if max(other._y) <= max(self._y) and min(other._y) >= min(self._y):
return True
else:
return False
kvadrat = [np.array([0,0]), np.array([1,0]), np.array([1,1]), np.array([0,1])]
kvadrat2 = [np.array([0.1,0.1]), np.array([0.9,0.1]), np.array([0.9,0.9]), np.array([0.1,0.9])]
A = SpecialRectangle(kvadrat)
B = SpecialRectangle(kvadrat2)
A.plot()
B.plot()
plt.gca().set_aspect('equal') # Samma skala för x,y-axel
B in A # True
A in B # FalseExercise 12Class for testing
from scipy import *
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
import scipy.linalg as sl
import scipy.sparse as sp
def crazzy_plus(a, b):
if not isinstance(a, int) or not isinstance(b, int):
raise ValueError('Endast heltal får användas.')
pass
first = str(a-b)
last = str(a+b)
svar = int(first+last)
print(a, '+', b, '=', svar)
return svarUse this code to test the class crazzy_plus from exercise12.pyMust be in the same directory
from exercise12 import crazzy_plus
import unittest
class TestIdentity(unittest.TestCase):
def test(self):
result = crazzy_plus(5,1)
expected = 46
self.assertEqual(result, expected)
def Test2(self):
self.assertRaises(ValueError, crazzy_plus(2.2,8))
if __name__=='__main__':
unittest.main()Calculates the integrals using trapezoidals
from scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
""" Task 1. Arbitrary mathematical function."""
def f(x):
y = sin(x)+x
return y
""" Primitive function of the mathematical function above."""
def F(x):
y = -cos(x)+(x**2)/2
return y
""" Function ctrapezoidal calculates the integral I of f(x) defined above using
trapezoids which has the width h (delta-x). a is lower limit and b upper limit
of the integral. n is the amount of trapezoids created within the interval.
It also returns h, the delta-x which is needed for task 3."""
a = 0
b = 10
n = 20
def ctrapezoidal(f, a, b, n):
h = (b-a)/n
I = (f(a)+f(b))/2
for i in range(n):
I = f((a+i*h)) + I
I = I*h
return(I,h)
""" Task 2. Required accuracy of the integral as input tolerance.
Begins with making two trapezoidals then doubles the amount of every
iteration. The value for previous iteration is compared with the actual if
the absolute value of the difference which is the error estimation.
The loop has a maxi amount of iteration so it never does more than 20 even
if the required tolerance not is reached."""
def call(tolerance= 1.e-3):
n = 2
maxi = 20
current = ctrapezoidal(f, a, b, n)[0]
for i in range(maxi+1):
previous = current
n = n*2
current = ctrapezoidal(f, a, b, n)[0]
if abs(previous-current) <= tolerance:
break
print("Required accuracy: ",tolerance)
print("Reached accuracy: ",abs(previous-current))
print("Nr of trapezoidals used:",n-1)
return current
""" Exact mathematical solution of the integral of the function f(x) with
primitive function F(x) using the fundamental theorem of calculus."""
def integral(F, a, b):
F = F(b) - F(a)
return F
""" Task 3. Makes a graph for comparison between the relationship between
increased accuracy with increasing number of trapezoidals used for calculating
the integral of f(x). Taking input of number of points wanted in the graph
for every point a new calculation of the integral has to be done. To know
the true error the primitive function F(x) of the function has to be known
for calculating how much each calculated intregral deviates from the exact."""
def errorplot(points=20):
xplot = zeros(points)
yplot = zeros(points)
n = 2
for i in range(points):
yplot[i], xplot[i] = ctrapezoidal(f, a, b, n)
yplot[i] = abs(yplot[i]-integral(F, a, b))
n = n*2
plt.loglog(xplot, yplot)
plt.grid(which="both", axis="both", ls="-")
plt.title('Trapezoidal approximation accuracy')
plt.xlabel('Width of trapezoidals')
plt.ylabel('Error')
return
""" Task 4. Takes two input vectors which one is for x-coordinates and the
other is for the y-coordinates of discrete points of an unknown function.
The vandermonde(v) returns the Vandermonde matrix which later is used
for creating an interpolation polynomial to fit this discrete points."""
xarray = array([0.0, 0.5, 1.0, 1.5, 2.0, 2.5])
yarray = array([-2.0, 0.5, -2.0, 1.0, -0.5, 1.0])
def vandermonde(v):
m = (len(v))
A = zeros((m,m))
for i in range(m):
A[:,m-i-1] = v**i
return A
""" Task 5. Takes the vandermonde matrix and y-values array and solves the
values for the coefficient vector c. A*c = y. Outside after the function
the Alexandre-Théophile Vandermonde matrix is stored as X and the coefficient
vector as c."""
def interpoly(A,y):
c = solve(A,y)
return(c)
X = vandermonde(xarray)
c = (interpoly(X,yarray))
""" Task 6. The polyval-function takes only the coefficient vector c as input
and calculates the interpolation polynomial value for the input value z.
In the for-loop the last coeff-element is first for lowest order so it goes
backwards why there is a minus sign in front of i."""
def polyval(c, z):
m = (len(c))
p = 0.
for i in range(m):
p += c[m-i-1]*z**i
return p
""" Task 7. Makes a plot of the interpolated polynomial and compares it with
the given point-values. x-values is in range from -0.02 to 2.57 in steps of
0.01. The discrete points is directly taken from xarray and yarray and for
plotting the interpolated polynomial polyval(c, z) is called to get the
polynomial values in the wanted range and steps."""
def vanderplot(xarray, yarray):
plt.title('Interpolated polynomial of discrete stars')
plt.xlabel('x-values')
plt.ylabel('y-values')
plt.plot(xarray, yarray, '*')
xval = np.arange(-0.02, 2.57, 0.01)
yval = polyval(c, xval)
plt.plot(xval, yval)
plt.grid(which="both", axis="both", ls="-")
returnInterval analysis using classesfrom scipy import *
from scipy.optimize import fsolve
from pylab import *
import numpy as np
import sys
import matplotlib.pyplot as plt
from scipy.integrate import quad
import scipy.linalg as sl
class Interval:
'''
The class Interval represents the left and right endpoints of types
integer or float in an closed interval.
'''
def __init__(self, left, right):
self.left, self.right = left, right
def __repr__(self):
return '[{}, {}]'.format(self.left, self.right)
def __add__(self, other):
a, b = self.left, self.right
if isinstance(other, Interval):
c, d = other.left, other.right
elif isinstance(other, int) or isinstance(other, float):
c, d = other, other
return Interval(a + c, b + d)
def __radd__(self, other):
return self + other
def __sub__(self, other):
a, b = self.left, self.right
if isinstance(other, Interval):
c, d = other.left, other.right
elif isinstance(other, int) or isinstance(other, float):
c, d = other, other
return Interval(a - d, b - c)
def __rsub__(self, other):
self.left, self.right = -self.right, -self.left
return self + other
def __mul__(self, other):
a, b = self.left, self.right
if isinstance(other, Interval):
c, d = other.left, other.right
elif isinstance(other, int) or isinstance(other, float):
c, d = other, other
return Interval(min(a*c, a*d, b*c, b*d), max(a*c, a*d, b*c, b*d))
def __rmul__(self, other):
return self * other
def __truediv__(self, other):
a, b = self.left, self.right
if isinstance(other, Interval):
c, d = other.left, other.right
elif isinstance(other, int) or isinstance(other, float):
c, d = other, other
if (c or d) == 0:
raise ZeroDivisionError
pass
else:
return Interval(min(a/c, a/d, b/c, b/d), max(a/c, a/d, b/c, b/d))
def __contains__(self, point):
if point >= self.left and point <= self.right:
return True
else:
return False
def __pow__(self, n):
a, b = self.left, self.right
if not isinstance(n, int):
raise TypeError('Power should be of type int.')
pass
if n%2 == 0: # n is even
if a >= 0:
return Interval(a**n, b**n)
elif b < 0:
return Interval(b**n, a**n)
else:
return Interval(0, max(a**n, b**n))
else: # n is odd
return Interval(a**n, b**n)
xl = linspace(0., 1, 1000) # 1000 intervals between
xu = linspace(0., 1, 1000)+0.5 # [0, 0.5] --> [1, 1.5]
lower = np.zeros(1000)
upper = np.zeros(1000)
for i in range(1000):
I = Interval(xl[i], xu[i])
P = 3*I**3 - 2*I**2 - 5*I - 1
lower[i] = P.left
upper[i] = P.right
plt.title('Interval analysis of a 3:rd degree polynomial')
plt.xlabel('I(x) = [x, x + ½]')
plt.ylabel('P(I) = 3I³ -2I² -5I -1')
plt.grid(which="both", axis="both", ls="-")
"""plt.plot(xl, lower, color='red')
plt.plot(xl, upper, color='red')"""
plt.fill_between(xl, xl, xu, color='grey', alpha=0.5)
plt.fill_between(xl, lower, upper, color='red', alpha=0.7)
"""plt.plot(xl, xl, color='grey')
plt.plot(xl, xu, color='grey')"""
I1 = Interval(1,4)
I2 = Interval(-2,-1)
I3 = Interval(0,3)