-
Notifications
You must be signed in to change notification settings - Fork 0
/
newton_raphson_matrix.py
133 lines (97 loc) · 3.01 KB
/
newton_raphson_matrix.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
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
from __future__ import division
import copy
from math import exp
from matrices import Matrix
E = 220e-3
R = 500
I_SA = 0.6e-6
I_SB = 1.2e-6
kT_q = 25e-3
EPSILON = 1e-9
def newton_raphson_matrix_solve():
"""
Solves for the nodal voltages of the nonlinear diode in Q3 by Newton-Raphson.
:return: the solved voltages, as well stepwise values for the error, voltage and f vector
"""
error_values = []
norm_values = []
vA_values = []
vB_values = []
fA_values = []
fB_values = []
iteration = 1
v_n = Matrix.empty(2, 1)
f = Matrix.empty(2, 1)
F = Matrix.empty(2, 2)
update_f(f, v_n)
f_0 = copy.deepcopy(f)
update_jacobian(F, v_n)
error_values.append('{:.3e}'.format(f.two_norm / f_0.two_norm))
norm_values.append('{:.3e}'.format(f.two_norm))
vA_values.append(v_n.scaled_values(1000)[0])
vB_values.append(v_n.scaled_values(1000)[1])
fA_values.append('{:.3e}'.format(f.values[0]))
fB_values.append('{:.3e}'.format(f.values[1]))
while f.two_norm / f_0.two_norm >= EPSILON:
v_n -= inverse_2x2(F) * f
update_f(f, v_n)
update_jacobian(F, v_n)
iteration += 1
error_values.append('{:.3e}'.format(f.two_norm / f_0.two_norm))
norm_values.append('{:.3e}'.format(f.two_norm))
vA_values.append(v_n.scaled_values(1000)[0])
vB_values.append(v_n.scaled_values(1000)[1])
fA_values.append('{:.3e}'.format(f.values[0]))
fB_values.append('{:.3e}'.format(f.values[1]))
return v_n, error_values, norm_values, vA_values, vB_values, fA_values, fB_values
def update_f(f, v_n):
"""
Updates the f vector.
:param f: the f vector to update
:param v_n: the nodal voltages
"""
v_a, v_b = v_n.values
f[0][0] = f_a(v_a, v_b)
f[1][0] = f_b(v_a, v_b)
def update_jacobian(F, v_n):
"""
Updates the Jacobian matrix.
:param F: the Jacobian matrix to update
:param v_n: the nodal voltages
"""
v_a, v_b = v_n.values
F[0][0] = dfa_dva(v_a, v_b)
F[0][1] = dfa_dvb(v_a, v_b)
F[1][0] = dfb_dva(v_a, v_b)
F[1][1] = dfb_dvb(v_a, v_b)
def f_a(v_a, v_b):
return v_a + R * I_SA * exp_f_term(v_a, v_b) - E
def f_b(v_a, v_b):
return I_SA * exp_f_term(v_a, v_b) - I_SB * exp_f_term(0, -v_b)
def dfa_dva(v_a, v_b):
return 1 + R * I_SA * exp_df_term(v_a, v_b)
def dfa_dvb(v_a, v_b):
return - R * I_SA * exp_df_term(v_a, v_b)
def dfb_dva(v_a, v_b):
return I_SA * exp_df_term(v_a, v_b)
def dfb_dvb(v_a, v_b):
return - I_SA * exp_df_term(v_a, v_b) - I_SB * exp_df_term(0, -v_b)
def exp_f_term(v_a, v_b):
return exp((v_a - v_b) / kT_q) - 1
def exp_df_term(v_a, v_b):
return exp((v_a - v_b) / kT_q) / kT_q
def inverse_2x2(A):
"""
Inverts a 2x2 matrix and returns a copy.
:param A: the matrix to invert
:return: the inverted matrix
"""
a = A[0][0]
b = A[0][1]
c = A[1][0]
d = A[1][1]
inverse = Matrix([
[d, -b],
[-c, a]
])
return inverse.scalar_divide(a * d - b * c)