-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathLabExam.py
More file actions
220 lines (174 loc) Β· 5.84 KB
/
LabExam.py
File metadata and controls
220 lines (174 loc) Β· 5.84 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
"""Converted from LabExam.ipynb"""
# 1. To find gradient of π = π₯2π¦ + 2π₯π§ β 4
#To find gradient of scalar point function.
from sympy.vector import *
from sympy import symbols
N=CoordSys3D('N')
x,y,z=symbols('x y z')
A=N.x**2*N.y+2*N.x*N.z-4
delop=Del()
display(delop(A))
gradA=gradient(A)
print(f"\n Gradient of {A} is \n")
display(gradA)
# 2. To find curl of πΉβ = π₯2π¦π§πΜ+ π¦2π§π₯πΜ+ π§2π₯π¦οΏ½
#To find curl of a vector point function
from sympy . vector import *
from sympy import symbols
N= CoordSys3D ('N')
x,y,z= symbols ('x y z')
A=N.x ** 2*N.y*N.z*N.i+N.y ** 2*N.z*N.x*N.j+N.z ** 2*N.x*N.y*N.k
delop =Del ()
curlA = delop . cross (A)
display ( curlA )
print (f"\n Curl of {A} is \n")
display ( curl (A))
# 3. To find gradient of πΉβ = π₯2π¦π§πΜ+ π¦2π§π₯πΜ+ π§2π₯π¦πΜ
#To find divergence of F=x^2yi+yz^2j+x^2zk
from sympy . physics . vector import *
from sympy import var
var ('x,y,z')
v= ReferenceFrame ('v')
F=v[0] ** 2*v[1]*v.x+v[1]*v[2] ** 2*v.y+v[0] ** 2*v[2]*v.z
G= divergence (F,v)
F=F. subs ([(v[0],x) ,(v[1],y) ,(v[2],z)])
print (" Given vector point function is ")
display (F)
G=G. subs ([(v[0],x) ,(v[1],y) ,(v[2],z)])
print (" Divergence of F=")
display (G)
# 4. Obtain a root of the equation π₯3 β 2π₯ β 5 = 0 between 2 and 3 by regulafalsi method. Perform 5 iterations.
# Regula Falsi method
from sympy import *
x= Symbol ('x')
g = input ('Enter the function ') #%x^3-2*x-5; % function
f= lambdify (x,g)
a= float ( input ('Enter a valus :')) #2
b= float ( input ('Enter b valus :')) # 3
N=int( input ('Enter number of iterations :')) #5
for i in range (1,N+1):
c=(a*f(b)-b*f(a))/(f(b)-f(a))
if ((f(a)*f(c)<0)):
b=c
else :
a=c
print ('itration %d \t the root %0.3f \t function value %0.3f \n'%(i,c,f(c)))
# 5. Find a root of the equation 3π₯ = cos π₯ + 1, near 1, by Newton Raphson method. Perform 5 iterations.
from sympy import *
x= Symbol ('x')
g = input ('Enter the function ') #%3x -cos(x)-1; % function
f= lambdify (x,g)
dg = diff (g);
df= lambdify (x,dg)
x0= float ( input ('Enter the intial approximation ')); # x0=1
n= int( input ('Enter the number of iterations ')); #n=5;
for i in range (1,n+1):
x1 =(x0 - (f(x0)/df(x0)))
print ('itration %d \t the root %0.3f \t function value %0.3f \n'%(i, x1 ,f(x1))); # print all iteration value
x0 = x1
# 6.1 Using Trapezoidal Rule Evaluate β«11+π₯2
# Definition of the function to integrate
def my_func(x):
return 1 / (1 + x ** 2)
# Function to implement trapezoidal method
def trapezoidal(x0, xn, n):
h = (xn - x0) / n # Calculating step size
# Finding sum
integration = my_func(x0) + my_func(xn) # Adding first and last terms
for i in range(1, n):
k = x0 + i * h # i-th step value
integration = integration + 2 * my_func(k) # Adding areas of the trapezoids
# Proportioning sum of trapezoid areas
integration = integration * h / 2
return integration
a = 0 # lower limit
b = 6 # upper limit
n = 6 # number of sub intervals
result = trapezoidal(a, b, n)
print('%3.5f' % result)
# 6.2 simpson's 1/3 Evaluate β«11+π₯2
# Function to implement the Simpson's one-third rule
def simpson13(x0, xn, n):
h = (xn - x0) / n # calculating step size
# Finding sum
integration = (my_func(x0) + my_func(xn))
k = x0
for i in range(1, n):
if i % 2 == 0:
integration = integration + 4 * my_func(k)
else:
integration = integration + 2 * my_func(k)
k += h
# Finding final integration value
integration = integration * h * (1 / 3)
return integration
# Input section
lower_limit = float(input("Enter lower limit of integration: "))
upper_limit = float(input("Enter upper limit of integration: "))
sub_interval = int(input("Enter number of sub intervals: "))
result = simpson13(lower_limit, upper_limit, sub_interval)
print("Integration result by Simpson's 1/3 method is: %0.6f" % result)
# 6.3 simpson's 3/8 Evaluate β«11+π₯2
def simpsons_3_8_rule(f, a, b, n):
h = (b - a) / n
s = f(a) + f(b)
for i in range(1, n, 3):
s += 3 * f(a + i * h)
for i in range(3, n - 1, 3):
s += 3 * f(a + i * h)
for i in range(2, n - 2, 3):
s += 2 * f(a + i * h)
return s * 3 * h / 8
def f(x):
return 1 / (1 + x ** 2) # function here
a = 0 # lower limit
b = 6 # upper limit
n = 6 # number of sub intervals
result = simpsons_3_8_rule(f, a, b, n)
print('%3.5f' % result)
# 7. Runge-Kutta methodApply the Runge Kutta method to find the solution of dy/dx = 1 + (y/x) at y(2) taking h = 0.2. Given that y(1) = 2.
from sympy import *
import numpy as np
def RungeKutta(g, x0, h, y0, xn):
x, y = symbols('x y')
f = lambdify([x, y], g)
xt = x0 + h
Y = [y0]
while xt <= xn:
k1 = h * f(x0, y0)
k2 = h * f(x0 + h / 2, y0 + k1 / 2)
k3 = h * f(x0 + h / 2, y0 + k2 / 2)
k4 = h * f(x0 + h, y0 + k3)
y1 = y0 + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
Y.append(y1)
x0 = xt
y0 = y1
xt = xt + h
return np.round(Y, 2)
RungeKutta('1+(y/x)', 1, 0.2, 2, 2)
# 8. Apply Milneβs predictor and corrector method to solve dy/dx = x2 + (y/2) at y(1.4).Given that y(1)=2, y(1.1)=2.2156, y(1.2)=2.4649, y(1.3)=27514. Use corrector formula thrice.
# Milne's method to solve first-order DE
# Use corrector formula thrice
x0 = 1
y0 = 2
y1 = 2.2156
y2 = 2.4649
y3 = 2.7514
h = 0.1
x1 = x0 + h
x2 = x1 + h
x3 = x2 + h
x4 = x3 + h
def f(x, y):
return x ** 2 + (y / 2)
y10 = f(x0, y0)
y11 = f(x1, y1)
y12 = f(x2, y2)
y13 = f(x3, y3)
y4p = y0 + (4 * h / 3) * (2 * y11 - y12 + 2 * y13)
print('Predicted value of y4 is %3.3f' % y4p)
y14 = f(x4, y4p)
for i in range(1, 4):
y4 = y2 + (h / 3) * (y14 + 4 * y13 + y12)
print('Corrected value of y4 after iteration %d is %3.5f' % (i, y4))
y14 = f(x4, y4)