maths

December 26, 2023

[ ]: import matplotlib.pyplot as plt

x = [1,2,3,4,6,7,8]
y = [2,7,9,1,5,10,3]
plt.scatter(x,y)
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.title('Scatter points')
plt.show()

1
[ ]: import matplotlib.pyplot as plt

x = [1,2,3,4,6,7,8]
y = [2,7,9,1,5,10,3]
plt.plot(x,y,'g+--')
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.title('My first graph!')
plt.show()

[ ]: import numpy as np
import matplotlib.pyplot as plt

x=np.arange(-10,10,0.001)
y1=np.sin(x)
y2=np.cos(x)
plt.plot(x,y1,'r',x,y2,'g',)
plt.title("sine and cosine curve ")
plt.xlabel("values of x")
plt.ylabel("values of sin(x) and cos(x)")

2
 plt.grid()
plt.show()

[ ]: import numpy as np
import matplotlib.pyplot as plt

x =np.arange(-10,10,0.001)

y =np.exp(x)
plt.plot(x,y)
plt.title("Exponential Curve")
plt.grid()
plt.show()

3
[ ]: import numpy as np
import matplotlib.pyplot as plt
x=np.linspace(0,2,100)

plt.plot(x,x,'r',label='linear')
plt.plot(x,x**2,'g',label='quadratic')
plt.plot(x,x**3,'b',label='cubic')
plt.xlabel('x-axis')
plt.ylabel('y-axis')
plt.title("Simple Plot")
plt.legend()
plt.show()

4
[ ]: from sympy import plot_implicit,symbols,Eq
x,y=symbols('x y')
p1=plot_implicit(Eq(x**2 + y**2,4),(x,-4,4),(y,-4,4),title='Circle:$x^2+y^2=4$')

5
[ ]: from sympy import plot_implicit,symbols,Eq
x,y=symbols('x y')
p3=␣
↪plot_implicit(Eq((y**2)*(2-x),(x**2)*(2+x)),(x,-5,5),(y,-5,5),title='Strophoid:

↪$y^2(a-x)=x^2 (a+x), a> 0$')

6
[ ]: from sympy import plot_implicit,symbols,Eq
x,y=symbols('x y')
p4=plot_implicit(Eq((y**2)*(3-x),x**3),(x,-2,5),(y,-5,5),title='Cissiod:
↪$y^2(a-x)=x^2,a>0$')

7
[ ]: from sympy import plot_implicit,symbols,Eq
x,y=symbols('x y')
p5=plot_implicit(Eq(4*(y**2),(x**2)*(4-x**2)),(x,-5,5),(y,-5,5),title='Lemniscate:
↪$a^2y^2=x^2(a^2-x^2)$')

8
[ ]: from sympy import plot_implicit,symbols,Eq
x,y=symbols('x y')
p6=plot_implicit(Eq(x**3+y**3,3*2*x*y),(x,-5,5),(y,-5,5),title='Folium of␣
↪De-Cartes:$x^3+y^3=3axy$')

9
[ ]: import numpy as np
import matplotlib.pyplot as plt

plt.axes(projection = 'polar')
r=3
rads = np.arange(0,(2 * np.pi),0.01)
for i in rads:
plt.polar(i,r,'g.')
plt.show()

10
[ ]: from pylab import *
theta=linspace(0,2*np.pi,1000)
r1=5+5*cos(theta)

polar(theta,r1,'r')
show()

11
[ ]: from pylab import *
theta=linspace(0,2*np.pi,1000)
r=2*abs(cos(2*theta))
polar(theta,r,'r')
show()

12
[ ]: import numpy as np
import matplotlib.pyplot as plt
import math

plt.axes(projection = 'polar')
a=3

rad = np.arange(0,(2*np.pi),0.01)
for i in rad:
r = a + (a*np.cos(i))
plt.polar(i,r,'g.')
r1 = a - (a*np.cos(i))
plt.polar(i,r1,'r.')
plt.show()

13
[ ]: import numpy as np
import matplotlib.pyplot as plt
def circle(r):
x = []
y = []

for theta in np.linspace(-2*np.pi,2*np.pi,100):

x.append(r*np.cos(theta))
y.append(r*np.sin(theta))
plt.plot(x,y)
plt.show()
circle(5)

14
[ ]: def cycloid(r):
x = []
y = []
for theta in np.linspace(-2*np.pi,2*np.pi,100):
x.append(r*(theta-np.sin(theta)))
y.append(r*(1-np.cos(theta)))

plt.plot(x,y)
plt.show()
cycloid(5)

15
[ ]: from sympy import *
r,t = symbols('r,t')

r1=4*(1+cos(t));
r2=5*(1-cos(t));
dr1=diff(r1,t)
dr2=diff(r2,t)
t1=r1/dr1
t2=r2/dr2
q=solve(r1-r2,t)

w1=t1.subs({t:float(q[1])})
w2=t2.subs({t:float(q[1])})
y1=atan(w1)
y2=atan(w2)
w=abs(y1-y2)
print('Angle between curves in radius is %0.3f'%(w))

Angle between curves in radius is 1.571

[ ]: from sympy import *

r,t = symbols('r,t')

16
 r1=4*(cos(t));
r2=5*(sin(t));
dr1=diff(r1,t)
dr2=diff(r2,t)
t1=r1/dr1
t2=r2/dr2
q=solve(r1-r2,t)

w1=t1.subs({t:float(q[0])})
w2=t2.subs({t:float(q[0])})
y1=atan(w1)
y2=atan(w2)
w=abs(y1-y2)
print('Angle between curves in radius is %0.4f'%(w))

Angle between curves in radius is 1.5708

[ ]: from sympy import *

t=Symbol('t')
r=Symbol('r')
r=4*(1+cos(t))
r1=Derivative(r,t).doit()
r2=Derivative(r1,t).doit()
rho=(r**2 + r1**2)**(1.5)/(r**2 + 2*r1**2 - r*r2)
rho1=rho.subs(t,pi/2)
print('The radius of curvature is %3.4f units'%rho1)

The radius of curvature is 3.7712 units

[ ]: from sympy import *

t,r,a,n=symbols('t r a n')
r=a*sin(n*t)
r1=Derivative(r,t).doit()
r2=Derivative(r1,t).doit()
rho=(r**2 + r1**2)**(1.5)/(r**2 + 2*r1**2 - r*r2)
rho1=rho.subs(t,pi/2)
rho1=rho.subs(n,1)
print('the radius of curvature is')
display(simplify(rho1))

the radius of curvature is

0.5
(𝑎2 )
2
[ ]: from sympy import *
from sympy.abc import rho,x,y,r,K,t,a,b,c,alpha

17
 y=(sqrt(x)-4)**2
y=a*sin(t)
x=a*cos(t)
dydx=simplify(Derivative(y,t).doit())/simplify(Derivative(x,t).doit())
rho=simplify((1+dydx**2)**1.5/(Derivative(dydx,t).doit()/(Derivative(x,t).
↪doit())))

print('Radius of curvature is')

display(ratsimp(rho))
t1=pi/2
r1=5
rho1=rho.subs(t,t1);
rho2=rho1.subs(a,r1);
print('\n\nRadius of curvature at r=5 and t= pi/2 is',simplify(rho2));
curvature=1/rho2;
print('\n\n Curvature at (5,pi/2) is',float(curvature))

Radius of curvature is
0.5
1
−𝑎 ( ) sin (𝑡)
sin2 (𝑡)

Radius of curvature at r=5 and t= pi/2 is -5

Curvature at (5,pi/2) is -0.2

[ ]: from sympy import *

from sympy.abc import rho,x,y,r,K,t,a,b,c,alpha
y=(a*sin(t))**(3/2)
x=(a*cos(t))**(3/2)
dydx=simplify(Derivative(y,t).doit())/simplify(Derivative(x,t).doit())
rho=simplify((1+dydx**2)**1.5/(Derivative(dydx,t).doit()/(Derivative(x,t).
↪doit())))

print('Radius of curvature is ')

display(ratsimp(rho))
t1=pi/4
r1=1;
rho1=rho.subs(t,t1);
rho2=rho1.subs(a,r1);
display('Radius of curvature at r=1 and t= pi/4 is',simplify(rho2));
curvature=1/rho2;
print('\n\n Curvature at (1,pi/4) is',float(curvature))

Radius of curvature is

18
 3.0 1.5
3.0
3.0 (𝑎 cos (𝑡)) ( (𝑎 sin (𝑡))
3.0 + 1) sin2 (𝑡) tan2 (𝑡)
(𝑎 cos (𝑡)) tan4 (𝑡)
− 1.5
(𝑎 sin (𝑡))
'Radius of curvature at r=1 and t= pi/4 is'
−2.52268924576114

Curvature at (1,pi/4) is -0.39640237166757364

[ ]: from sympy import *

x,y=symbols('x y')

u=exp(x)*(x*cos(y)-y*sin(y))
duxy=diff(u,x,y)
duyx=diff(u,y,x)
if duxy==duyx:
print('Mixed partial derivatives are equal')
else:
print('Mixed partial derivatives are not equal')

Mixed partial derivatives are equal

[ ]: from sympy import *

x,y=symbols('x y')
u=exp(x)*(x*cos(y)-y*sin(y))
display(u)
dux=diff(u,x)
duy=diff(u,y)
uxx=diff(dux,x)
uyy=diff(duy,y)
w=uxx+uyy
w1=simplify(w)
print('Ans:',float(w1))

(𝑥 cos (𝑦) − 𝑦 sin (𝑦)) 𝑒𝑥

Ans: 0.0

[ ]: from sympy import *

x,y,z=symbols('x,y,z')
u=x*y/z
v=y*z/x
w=z*x/y
dux=diff(u,x)
duy=diff(u,y)
duz=diff(u,z)
dvx=diff(v,x)

19
 dvy=diff(v,y)
dvz=diff(v,z)
dwx=diff(w,x)
dwy=diff(w,y)
dwz=diff(w,z)
J=Matrix([[dux,duy,duz],[dvx,dvy,dvz],[dwx,dwy,dwz]]);
print("The Jacobian matrix is\n")
display(J)
Jac=det(J).doit()
print('\n\nJ=',Jac)

The Jacobian matrix is

𝑦 𝑥
− 𝑥𝑦
𝑧2
⎡−𝑧𝑦𝑧 𝑧
𝑧 𝑦 ⎤
⎢ 𝑥2 𝑥 𝑥 ⎥
𝑧
⎣ 𝑦 − 𝑥𝑧
𝑦2
𝑥
𝑦 ⎦

J= 4

[ ]: from sympy import *

x,y,z=symbols('x,y,z')
u=x+3*y**2-z**3
v=4*x**2*y*z
w=2*z*z**2-x*y
dux=diff(u,x)
duy=diff(u,y)
duz=diff(u,z)
dvx=diff(v,x)
dvy=diff(v,y)
dvz=diff(v,z)
dwx=diff(w,x)
dwy=diff(w,y)
dwz=diff(w,z)
J=Matrix([[dux,duy,duz],[dvx,dvy,dvz],[dwx,dwy,dwz]]);
print("The Jacobian matrix is\n")
display(J)
Jac=det(J).doit()
print('\n\nJ=\n')

display(Jac)
J1=J.subs([(x,1),(y,-1),(z,0)])
print('\n\nJ at (1,-1,0):\n')
Jac1=Determinant(J1).doit()
display(Jac1)

20
 The Jacobian matrix is

1 6𝑦 −3𝑧 2
⎡8𝑥𝑦𝑧 4𝑥2 𝑧 4𝑥2 𝑦 ⎤
⎢ ⎥
⎣ −𝑦 −𝑥 6𝑧 2 ⎦

J=

4𝑥3 𝑦 − 24𝑥2 𝑦3 + 12𝑥2 𝑦𝑧 3 + 24𝑥2 𝑧3 − 288𝑥𝑦2 𝑧3

J at (1,-1,0):

20
[ ]: from sympy import *
from sympy.abc import rho,phi,theta
X=rho*cos(phi)*sin(theta);
Y=rho*cos(phi)*cos(theta);
Z=rho*sin(phi);
dx=Derivative(X,rho).doit()
dy=Derivative(Y,rho).doit()
dz=Derivative(Z,rho).doit()
dx1=Derivative(X,phi).doit()
dy1=Derivative(Y,phi).doit()
dz1=Derivative(Z,phi).doit()
dx2=Derivative(X,theta).doit()
dy2=Derivative(Y,theta).doit()
dz2=Derivative(Z,theta).doit()
J=Matrix([[dx,dy,dz],[dx1,dy1,dz1],[dx2,dy2,dz2]]);
print("The Jacobian matrix is")
display(J)
print('\n\nJ=\n')
display(simplify(Determinant(J).doit()))

The Jacobian matrix is

sin (𝜃) cos (𝜙) cos (𝜙) cos (𝜃) sin (𝜙)
⎡−𝜌 sin (𝜙) sin (𝜃) −𝜌 sin (𝜙) cos (𝜃) 𝜌 cos (𝜙)⎤
⎢ ⎥
⎣ 𝜌 cos (𝜙) cos (𝜃) −𝜌 sin (𝜃) cos (𝜙) 0 ⎦

J=

21
 𝜌2 cos (𝜙)

[ ]: from sympy import *

x,y=symbols('x y')
f=x**2+x*y+y**2+3*x-3*y+4
d1=diff(f,x)
d2=diff(f,y)
cp=solve((d1,d2),(x,y))
print('critical points are:',cp)
s1=diff(f,x,2)
s2=diff(f,y,2)
s3=diff(diff(f,y),x)
q1=s1.subs(cp)
q2=s2.subs(cp)
q3=s3.subs(cp)
delta=s1*s2-s3**2
print('funtion value is',delta)
if(delta>0 and s1<0):
print("f takes maximum")
elif(delta>0 and s1>0):
print("f takes minimum")
if(delta<0):
print("the point is a saddle point")
if(delta==0):
print("further tests required")

critical points are: {x: -3, y: 3}

funtion value is 3
f takes minimum

[ ]: import numpy as np
from matplotlib import pyplot as plt
from sympy import*
x=Symbol('x')
y=sin(1*x)
format
xo=float(pi/2)
dy=diff(y,x)
d2y=(diff(y,x,2))
d3y=(diff(y,x,3))
yat=lambdify(x,y)
dyat=lambdify(x,dy)
d2yat=lambdify(x,d2y)
d3yat=lambdify(x,d3y)
y=yat(xo)+((x-xo)/2)*dyat(xo)+((x-xo)**2/6)*d2yat(xo)+((x-xo)**3/24)*d3yat(xo)
print(simplify(y))
yat=lambdify(x,y)
print("%.3f"%yat(pi/2+10*(pi/180)))

22
 def f(x):
return np.sin(1*x)
x=np.linspace(-10,10)
plt.plot(x,yat(x),color='red')
plt.plot(x,f(x),color='green')
plt.ylim([-3,3])
plt.grid()
plt.show()

-2.55134749822365e-18*x**3 - 0.166666666666667*x**2 + 0.523598775598299*x +

0.588766483287943
0.995

[ ]: import numpy as np
A=np.matrix([[1,2,-1],[2,1,4],[3,3,4]])
B=np.matrix([[0],[0],[0]])
r=np.lining.matrix_rank(A)
n=A.shape[1]
if(r==n):
print("System has trivial solution")
else:
print("system has",n-r,"non tr")

23
 System has trivial solution

[ ]: f1=lambda x,y,z:(17-y+2*z)/20
f2=lambda x,y,z:(-18-3*x+z)/20
f3=lambda x,y,z:(25-2*x+3*y)/20
x0=0
y0=0
z0=0
count=1
e=float(input('Enter tolerable error:'))
print('\ncount\tx\ty\tz\n')
condition=True
while condition:
x1=f1(x0,y0,z0)
y1=f2(x1,y0,z0)
z1=f3(x1,y1,z0)
print('%d\t%0.4f\t%0.4f\t%0.4f\n'%(count,x1,y1,z1))
e1=abs(x0-x1);
e2=abs(y0-y1);
e3=abs(z0-z1);
count+=1
x0=x1
y0=y1
z0=z1
condition=e1>e and e2>e and e3>e
print('\nSolution: x=%0.3f,y=%0.3f and z=%0.3f\n'%(x1,y1,z1))

Enter tolerable error:0.0001

count x y z

1 0.8500 -1.0275 1.0109

2 1.0025 -0.9998 0.9998

3 1.0000 -1.0000 1.0000

4 1.0000 -1.0000 1.0000

Solution: x=1.000,y=-1.000 and z=1.000

[ ]: f1=lambda x,y,z:(17-y+2*z)/20
f2=lambda x,y,z:(-18-3*x+z)/20
f3=lambda x,y,z:(25-2*x+3*y)/20
x0=0

24
 y0=0
z0=0
count=1
e=float(input('Enter tolerable error:'))
print('\ncount\tx\ty\tz\n')
condition=True
while condition:
x1=f1(x0,y0,z0)
y1=f2(x1,y0,z0)
z1=f3(x1,y1,z0)
print('%d\t%0.4f\t%0.4f\t%0.4f\n'%(count,x1,y1,z1))
e1=abs(x0-x1);
e2=abs(y0-y1);
e3=abs(z0-z1);
count+=1
x0=x1
y0=y1
z0=z1
condition=e1>e and e2>e and e3>e
print('\nSolution: x=%0.3f,y=%0.3f and z=%0.3f\n'%(x1,y1,z1))

[ ]: import numpy as np
I=np.array([[4,3,2],[1,4,1],[3,10,4]])
print("\n Given matrix:\n",I)
w,v=np.linalg.eig(I)
print("\n Eigen values:\n",w)
print("\n Eigen vectors:\n",v)
print("Eigen value:\n",w[0])
print("\n Corresponding Eigen vector:",v[:,0])

Given matrix:
[[ 4 3 2]
[ 1 4 1]
[ 3 10 4]]

Eigen values:
[8.98205672 2.12891771 0.88902557]

Eigen vectors:
[[-0.49247712 -0.82039552 -0.42973429]
[-0.26523242 0.14250681 -0.14817858]
[-0.82892584 0.55375355 0.89071407]]
Eigen value:
8.982056720677654

Corresponding Eigen vector: [-0.49247712 -0.26523242 -0.82892584]

25
[ ]: import numpy as np
I=np.array([[1,-3,3],[3,-5,3],[6,-6,4]])
print("\n Given matrix:\n",I)
w,v=np.linalg.eig(I)
print("\n Eigen values:\n",w)
print("\n Eigen vectors:\n",v)

Given matrix:
[[ 1 -3 3]
[ 3 -5 3]
[ 6 -6 4]]

Eigen values:
[ 4.+0.00000000e+00j -2.+1.10465796e-15j -2.-1.10465796e-15j]

Eigen vectors:
[[-0.40824829+0.j 0.24400118-0.40702229j 0.24400118+0.40702229j]
[-0.40824829+0.j -0.41621909-0.40702229j -0.41621909+0.40702229j]
[-0.81649658+0.j -0.66022027+0.j -0.66022027-0.j ]]

[ ]: from sympy import *

x,y=symbols('x,y')
y=Function("y")(x)
y1=Derivative(y,x)
z1=dsolve(Eq(y1+y*tan(x)-y**3*sec(x),0),y)
display(z1)

[Eq(y(x), -sqrt(1/(C1 - 2*sin(x)))*cos(x)),

Eq(y(x), sqrt(1/(C1 - 2*sin(x)))*cos(x))]

[ ]: from sympy import *

x,y=symbols('x,y')
y=Function("y")(x)
y1=Derivative(y,x)
z1=dsolve(Eq(x**3*y1-x**2*y+y**4*cos(x),0),y)
display(z1)

[Eq(y(x), (x**3/(C1 + 3*sin(x)))**(1/3)),

Eq(y(x), (x**3/(C1 + 3*sin(x)))**(1/3)*(-1 - sqrt(3)*I)/2),
Eq(y(x), (x**3/(C1 + 3*sin(x)))**(1/3)*(-1 + sqrt(3)*I)/2)]

[ ]: f1=lambda x,y,z:(17-y+2*z)/20
f2=lambda x,y,z:(-18-3*x+z)/20
f3=lambda x,y,z:(25-2*x+3*y)/20
x0=0
y0=0
z0=0

26
 count=1
e=float(input("Enter tolerate error"))
print('count\tx\ty\tz\n')
condition=True
while condition:
x1=f1(x0,y0,z0)
y1=f2(x1,y0,z0)
z1=f3(x1,y1,z0)
print('%d\t%0.4f\t%0.4f\t%0.4f\n'%(count ,x1,y1,z1))
e1=abs(x0-x1);
e2=abs(y0-y1);
e3=abs(z0-z1);
count+=1
x0=x1
y0=y1
z0=z1
condition=e1>e and e2>e and e3>e
print('\n solution : x=% 0.3f,y=% 0.3f and z=% 0.3f\n' %(x1,y1,z1))

Enter tolerate error0.001

count x y z

1 0.8500 -1.0275 1.0109

2 1.0025 -0.9998 0.9998

3 1.0000 -1.0000 1.0000

solution : x= 1.000,y=-1.000 and z= 1.000

[ ]: from sympy import*

r,t=symbols('r,t')
r1=4*(1+cos(t))
r2=5*(1-cos(t))
dr1=diff(r1,t)
dr2=diff(r2,t)
t1=r1/dr1
t2=r2/dr2
q=solve(r1-r2,t)
w1=t1.subs({t:float(q[1])})
w2=t2.subs({t:float(q[1])})
y1=atan(w1)
y2=atan(w2)
w=abs(y1-y2)
print("Angle between curves in radius is%0.3f"%(w))

27
Angle between curves in radius is1.571

28