-
Notifications
You must be signed in to change notification settings - Fork 0
/
Newton_Raphson_Method.py
95 lines (74 loc) · 1.96 KB
/
Newton_Raphson_Method.py
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
# -*- coding: utf-8 -*-
"""
Created on Sun Aug 30 11:55:30 2020
@author: Kari Ness
"""
"""
Newton-Raphson method (N-R)
"""
import sympy
import numpy as np
import matplotlib.pyplot as plt
#code terminates based on number of iterations
def iterations(num_it,f,D,P,d):
table = np.zeros((num_it,6))
#current iteration,i
i = 0;
while i < num_it:
#iteration in table
table[i,0] = i
deltaD,Kt,r_num = NewtonRaphson(f,D,P,d);
#displacement in table
table[i,1]=d
#deltaD in table
table[i,5] = deltaD
#Kt in table
table[i,4] = Kt
#residual in table
table[i,2] = r_num
#convergence criteria
table[i,3] = r_num/P
d = d+deltaD;
i += 1;
print(table)
return table
#code terminates based on tolerance
#def tolerances(tolerance):
def NewtonRaphson(f,D,P,d):
#finds residual
r = P-f;
#numerical values of r
r_num = r.subs({D:d});
#derivative of residual, drdd
drdd = sympy.diff(r,D);
#finds tangent stiffness
Kt = -drdd.subs({D:d});
#finds deltaD
deltaD = r_num/Kt;
return deltaD,Kt,r_num
def plot_eq_path(f,table,D):
drdd = sympy.diff(f,D);
d = np.linspace(table[0,1],table[-1,1]*12,100)
f_num = []
Kt = []
for i in range(0,100):
f_num.append(f.subs({D:d[i]}));
Kt.append(drdd.subs({D:d}))
plt.plot(d,f_num)
limit_points = sympy.solve(Kt, D)
print('Limit points: ')
for i in range(len(limit_points)):
print(str(limit_points[i][0]) + ','+ str(f.subs({D:limit_points[i][0]})))
zero_points = sympy.solve(f, D)
print('Zero points: ' + str(zero_points))
def main():
D = sympy.Symbol('D')
f = 1/2*(2*D-3*D**2+D**3);
#Target value, P
P = 0.125;
#dependent variable, d
d = 0;
num_it = 5;
table = iterations(num_it,f,D,P,d)
plot_eq_path(f,table,D)
main()