源起

经过两次电路程序设计作业的洗礼——谁叫我这个憨憨不想用Matlab呢（逃），对python的掌握更进了一步，也对电路的求解有了更深的认识。索性就把第二次实验，即设计程序求解高阶电路的过程和思路搬到博客上来其实是因为老师说放在博客上会有额外的加分，也算是对这些天工作的一个小小的总结

框架

程序采用了python的tkinter库来实现可视化操作，虽然界面比较朴素，但也还算是勉强实现了人机交互的功能。然后在使用了matplotlib,numpy以及scipy实现微分方程的求解和作图

这次大作业算是对我们手下留情了，所有的高阶原件都在同一个支路中，所以可以直接使用戴维宁等效化简再求解即可，从而免去了设状态方程的过程（主要是python求解微分方程必须要化成标准方程的形式，这样的话会使得问题的难度成倍增大）。

这样一来，只要解决了等效电路的求解，问题就迎刃而解了

建立系数矩阵+高斯消元

将高阶原件所在的支路首尾端点连起来，建立一个含有假想的新元件的支路

首先将这条支路上的原件设为值为0的电流源，代入求解开路电压uoc

然后再将该原件变为值为1欧姆的电阻，代入求得电压u1

即可得到req=uoc/u1−1req=uoc/u_1-1

具体求解时，即用节点电压法，以节点n为接地端，构造系数矩阵，然后结合高斯消元，得到的r列向量里的每个值，即对应每个点的节点电压

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
def rec(): #构建系数矩阵
	for i in range(1,n+1):
		R[i]=0.0
		for j in range(1,n+1):
			L[i][j]=0.0
	for i in range(1,n):
		global flag
		flag=0
		for j in range(1,n+1):
			x=min(i,j)
			y=max(j,i)			
			if flag:
				break
			if link[i][j]==0:
				continue
			for h in range(1,num[i][j]+1):
				if K[h][i][j]==-1:
					flag=1
					cx=conx[h][i][j]
					cy=cony[h][i][j]
					if K[1][cx][cy]==2:
						for q in range(1,n+1):
							L[i][q]=0
						L[i][x]+=1
						L[i][y]-=1
						R[i]=A[h][i][j]*A[1][cx][cy]
						break
					if K[1][cx][cy]==3:
						for q in range(1,n):
							L[i][q]=0
						L[i][x]+=1
						L[i][y]-=1
						L[i][cx]-=1/A[1][cx][cy]*A[h][i][j]
						L[i][cy]+=1/A[1][cx][cy]*A[h][i][j]
						R[i]=0
						break

				if K[h][i][j]==-2:
					cx=conx[h][i][j]
					cy=cony[h][i][j]
					L[i][cx]+=A[h][i][j]
					L[i][cy]-=A[h][i][j]

				if K[h][i][j]==1:
					flag=1
					for p in range(1,n+1):
						L[i][p]=0
					L[i][x]=1
					L[i][y]=-1
					R[i]=A[h][i][j]
					break
				if K[h][i][j]==2:
					R[i]+=(int(i>j)*2-1)*A[h][i][j]
				if K[h][i][j]==3:
					L[i][i]+=1.0/A[h][i][j]
					L[i][j]-=1.0/A[h][i][j]
	for i in range(1,n):
		L[n][i]=0
	L[n][n]=1

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
def solve():#高斯消元
	for i in range(1,n+1):
		r[i]=R[i]
	for i in range(1,n+1):
		for j in range(1,n+1):
			l[i][j]=L[i][j]
	for i in range(1,n+1):
		if abs(l[i][i])<=eps:
			for j in range(i+1,n+1):
				if abs(l[j][i])>eps:
					for k in range(1,n+1):
						t=l[i][k]
						l[i][k]=l[j][k]
						l[j][k]=t
					break
		if(abs(l[i][i])<=eps):
			continue
		t=l[i][i]
		for j in range(i,n+1):
			l[i][j]/=t
		r[i]/=t
		for j in range(i+1,n+1):
			t=l[j][i]
			for k in range(i,n+1):
				l[j][k]-=l[i][k]*t
			r[j]-=r[i]*t

	for i in range(n,0,-1):
		for j in range(i-1,0,-1):
			r[j]-=l[j][i]*r[i]
			l[j][i]=0.0

求解微分方程组

分一阶二阶及三阶的情况讨论，注意求解时要化为微分方程组的标准形式，其中三阶电路我默认为电感与电容串联之后再与电阻并联的形式，如下

avartar

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
def paint():
	if tot==1:
		def diff_equation(y,x):
			return np.array((uoc-y)/(req*cap[captot]))
		x=np.linspace(0,1,num=1000)
		result=odeint(diff_equation,0,x)
		plt.plot(x,result[:,0])
		plt.grid()
		plt.show()
	if tot==2:
		def diff_equation21(y_list,x):
			y,z=y_list
			return np.array([z,(uoc-y-req*cap[captot]*z)/(res[restot]*cap[captot])])
		def diff_equation22(y_list,x):
			y,z=y_list
			return np.array([((uoc-y)/req-z)/cap[captot],y/res[restot]])

		x=np.linspace(0,0.5,num=1000)
		y0=[0,0]
		if binglian==1:
			result=odeint(diff_equation22,y0,x)
		else :
			result=odeint(diff_equation21,y0,x)
		plt.plot(x,result[:,0],label='U_c')
		plt.plot(x,result[:,1],label='I_l')
		plt.legend()
		plt.grid()
		plt.show()
	if tot==3:
		def diff_equation3(y_list,x):
			y,z,w=y_list
			return np.array([w,((uoc-z)/req-cap[1]*w)/cap[2],(z-y)/(res[1]*cap[1])])
		x=np.linspace(0,100,num=100000)
		y0=[0,0,0]
		result=odeint(diff_equation3,y0,x)
		plt.plot(x,result[:,0],label='U_c2')
		plt.plot(x,result[:,1],label='U_c1')
		plt.plot(x,result[:,2],label='I_l')
		plt.legend()
		plt.grid()
		plt.show()

样例及结果

一阶电路

avartar
avartar

二阶电路

avartar
avartar

三阶电路

avartar
avartar

源代码

因为我对python掌握的还较为浅陋，很多东西都是现学现卖，所以代码风格可能有些辣眼睛，望轻喷~

这里的输入时的节点1和节点2代表含有电感电容原件所在的支路，对于电容和电感，只需输入各个参数大小即可，其在电路中的位置不用输入
输入总节点和总支路时，不考虑电容电感内的支路，如果一个节点所连的支路都是电容或者电感，则不考虑该节点
总节点数若有n个，则始终把接地当作n节点

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
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
import tkinter as tk
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import odeint
K=[[[0 for i in range(100)] for j in range(100) ] for h in range(100)]
A=[[[0.0 for i in range(100)] for j in range(100) ] for h in range(100)]
link=[([0]*100)for i in range(100)]
conx=[[[0 for i in range(100)] for j in range(100) ] for h in range(100)]
cony=[[[0 for i in range(100)] for j in range(100) ] for h in range(100)]
L=[([0.0]*100)for i in range(100)]
l=[([0.0]*100)for i in range(100)]
R=[0.0 for i in range(100)]
r=[0.0 for i in range(100)]
num=[([0]*100)for i in range(100)]
tmpa=[0.0 for i in range(100)]
tmpk=[0 for i in range(100)]
n=0
flag=0
cnt=0
eps=0.00000001
foot=0.001
captot=0
restot=0
res=[0.0 for i in range(100)]
cap=[0.0 for i in range(100)]
tot=0
p1=0
p2=0
uoc=0
req=0
tmp=0
binglian=0
def rec():
	for i in range(1,n+1):
		R[i]=0.0
		for j in range(1,n+1):
			L[i][j]=0.0
	for i in range(1,n):
		global flag
		flag=0
		for j in range(1,n+1):
			x=min(i,j)
			y=max(j,i)			
			if flag:
				break
			if link[i][j]==0:
				continue
			for h in range(1,num[i][j]+1):
				if K[h][i][j]==-1:
					flag=1
					cx=conx[h][i][j]
					cy=cony[h][i][j]
					if K[1][cx][cy]==2:
						for q in range(1,n+1):
							L[i][q]=0
						L[i][x]+=1
						L[i][y]-=1
						R[i]=A[h][i][j]*A[1][cx][cy]
						break
					if K[1][cx][cy]==3:
						for q in range(1,n):
							L[i][q]=0
						L[i][x]+=1
						L[i][y]-=1
						L[i][cx]-=1/A[1][cx][cy]*A[h][i][j]
						L[i][cy]+=1/A[1][cx][cy]*A[h][i][j]
						R[i]=0
						break

				if K[h][i][j]==-2:
					cx=conx[h][i][j]
					cy=cony[h][i][j]
					L[i][cx]+=A[h][i][j]
					L[i][cy]-=A[h][i][j]

				if K[h][i][j]==1:
					flag=1
					for p in range(1,n+1):
						L[i][p]=0
					L[i][x]=1
					L[i][y]=-1
					R[i]=A[h][i][j]
					break
				if K[h][i][j]==2:
					R[i]+=(int(i>j)*2-1)*A[h][i][j]
				if K[h][i][j]==3:
					L[i][i]+=1.0/A[h][i][j]
					L[i][j]-=1.0/A[h][i][j]
	for i in range(1,n):
		L[n][i]=0
	L[n][n]=1

def solve():
	for i in range(1,n+1):
		r[i]=R[i]
	for i in range(1,n+1):
		for j in range(1,n+1):
			l[i][j]=L[i][j]
	for i in range(1,n+1):
		if abs(l[i][i])<=eps:
			for j in range(i+1,n+1):
				if abs(l[j][i])>eps:
					for k in range(1,n+1):
						t=l[i][k]
						l[i][k]=l[j][k]
						l[j][k]=t
					break
		if(abs(l[i][i])<=eps):
			continue
		t=l[i][i]
		for j in range(i,n+1):
			l[i][j]/=t
		r[i]/=t
		for j in range(i+1,n+1):
			t=l[j][i]
			for k in range(i,n+1):
				l[j][k]-=l[i][k]*t
			r[j]-=r[i]*t

	for i in range(n,0,-1):
		for j in range(i-1,0,-1):
			r[j]-=l[j][i]*r[i]
			l[j][i]=0.0



def paint():
	if tot==1:
		def diff_equation(y,x):
			return np.array((uoc-y)/(req*cap[captot]))
		x=np.linspace(0,1,num=1000)
		result=odeint(diff_equation,0,x)
		plt.plot(x,result[:,0])
		plt.grid()
		plt.show()
	if tot==2:
		def diff_equation21(y_list,x):
			y,z=y_list
			return np.array([z,(uoc-y-req*cap[captot]*z)/(res[restot]*cap[captot])])
		def diff_equation22(y_list,x):
			y,z=y_list
			return np.array([((uoc-y)/req-z)/cap[captot],y/res[restot]])

		x=np.linspace(0,0.5,num=1000)
		y0=[0,0]
		if binglian==1:
			result=odeint(diff_equation22,y0,x)
		else :
			result=odeint(diff_equation21,y0,x)
		plt.plot(x,result[:,0],label='U_c')
		plt.plot(x,result[:,1],label='I_l')
		plt.legend()
		plt.grid()
		plt.show()
	if tot==3:
		def diff_equation3(y_list,x):
			y,z,w=y_list
			return np.array([w,((uoc-z)/req-cap[1]*w)/cap[2],(z-y)/(res[1]*cap[1])])
		x=np.linspace(0,100,num=100000)
		y0=[0,0,0]
		result=odeint(diff_equation3,y0,x)
		plt.plot(x,result[:,0],label='U_c2')
		plt.plot(x,result[:,1],label='U_c1')
		plt.plot(x,result[:,2],label='I_l')
		plt.legend()
		plt.grid()
		plt.show()		

def window2():
	window = tk.Tk()
	window.geometry("500x300")
	window.title("Second Window")
	lb1 = tk.Label(window, text="节点一")
	lb1.grid(column=0, row=0)
	s1 = tk.StringVar()
	txt1 = tk.Entry(window, width=12,textvariable=s1)
	txt1.grid(column=0, row=1)
	lb2 = tk.Label(window, text="节点二")
	lb2.grid(column=1, row=0)
	s2 = tk.StringVar()
	txt2 = tk.Entry(window, width=12,textvariable=s2)
	txt2.grid(column=1, row=1)
	lb3 = tk.Label(window, text="元件类型")
	lb3.grid(column=2, row=0)
	s3 = tk.StringVar()
	txt3 = tk.Entry(window, width=12,textvariable=s3)
	txt3.grid(column=2, row=1)
	lb4 = tk.Label(window, text="参数大小")
	lb4.grid(column=3, row=0)
	s4 = tk.StringVar()
	txt4 = tk.Entry(window, width=12,textvariable=s4)
	txt4.grid(column=3, row=1)
	lb5 = tk.Label(window, text="电压源即1，电流源即2，受控电压源为-1，受控电流源即-2，电阻即3")
	lb5.grid(columnspan=8,row=5)
	lb6 = tk.Label(window, text="注意输入时参数时须保证电流正方向为从小节点编号到大节点编号，电压正极为小编号")
	lb6.grid(columnspan=10, row=6)
	def check():
		global cnt		
		cnt+=1
		u=int(txt1.get())
		v=int(txt2.get())
		p=int(txt3.get())
		w=float(txt4.get())	
		num[u][v]+=1
		num[v][u]+=1
		if p<0:
			window = tk.Tk()
			window.geometry("350x200")
			window.title("输入控制源")
			lb31 = tk.Label(window, text="节点一")
			lb31.grid(column=0, row=0)
			s31 = tk.StringVar()
			txt31 = tk.Entry(window, width=10,textvariable=s31)
			txt31.grid(column=0, row=1)
			lb32 = tk.Label(window, text="节点二")
			lb32.grid(column=1, row=0)
			s32 = tk.StringVar()
			txt32 = tk.Entry(window, width=10,textvariable=s32)
			txt32.grid(column=1, row=1)			
			def scan():
				x=int(txt31.get())
				y=int(txt32.get())
				conx[num[u][v]][u][v]=min(x,y)
				conx[num[u][v]][v][u]=min(x,y)
				cony[num[u][v]][u][v]=max(x,y)
				cony[num[u][v]][v][u]=max(x,y)
				window.quit()
			button31 = tk.Button(window, text='记录',width=10,command=scan)
			button31.grid(row=3, column=3)		
			window.mainloop()	
		link[u][v]=1
		link[v][u]=1
		K[num[u][v]][u][v]=p
		K[num[u][v]][v][u]=p
		A[num[u][v]][u][v]=w
		A[num[u][v]][v][u]=w
		txt1.delete(0, 'end')
		txt2.delete(0, 'end')
		txt3.delete(0, 'end')
		txt4.delete(0, 'end')
		if cnt==m:
			global uoc
			global req
			link[p1][p2]=1
			link[p2][p1]=1
			num[p1][p2]+=1
			num[p2][p1]+=1
			K[num[p1][p2]][p1][p2]=2
			K[num[p1][p2]][p2][p1]=2
			A[num[p1][p2]][p1][p2]=0
			A[num[p1][p2]][p2][p1]=0
			rec()
			solve()
			uoc=r[p1]-r[p2]
			K[num[p1][p2]][p1][p2]=3
			K[num[p1][p2]][p2][p1]=3
			A[num[p1][p2]][p1][p2]=1
			A[num[p1][p2]][p2][p1]=1
			rec()
			solve()
			u1=r[p1]-r[p2]
			req=uoc/u1-1
			paint()

	button1 = tk.Button(window, text='记录',width=10,command=check)
	button1.grid(row=7, column=7)
	window.mainloop()	

tmp=0

def window3():
	window = tk.Tk()
	window.geometry("500x300")
	window.title("Second Window")
	lb3 = tk.Label(window, text="元件类型")
	lb3.grid(column=0, row=0)
	s3 = tk.StringVar()
	txt3 = tk.Entry(window, width=12,textvariable=s3)
	txt3.grid(column=0, row=1)
	lb4 = tk.Label(window, text="参数大小")
	lb4.grid(column=1, row=0)
	s4 = tk.StringVar()
	txt4 = tk.Entry(window, width=12,textvariable=s4)
	txt4.grid(column=1, row=1)
	lb5 = tk.Label(window, text="电容即1，电感即2")
	lb5.grid(column=0, row=2)

	def check():
		global tot
		global captot
		global restot	
		global tmp
		global binglian
		tmp+=1
		p=int(txt3.get())
		w=float(txt4.get())	
		if p==1:
			captot+=1
			cap[captot]=w
		if p==2:
			restot+=1
			res[restot]=w
		txt3.delete(0, 'end')
		txt4.delete(0, 'end')
		if tmp==tot:
			window.destroy()
	button1 = tk.Button(window, text='记录',width=10,command=check)
	button1.grid(row=5, column=0)
	window.mainloop()
def init():
	window = tk.Tk()
	window.title("简陋的交互界面T_T")
	window.geometry("500x500")
	lbl = tk.Label(window, text="请输入总节点数（除去电容电感支路）")
	lbl.grid(columnspan=6, row=0)
	s = tk.StringVar()
	txt = tk.Entry(window, width=10,textvariable=s)
	txt.grid(column=7, row=0)
	lbl2 = tk.Label(window, text="请输入总支路数（除去电容电感支路）")
	lbl2.grid(columnspan=6, row=1)
	s2 = tk.StringVar()
	txt2 = tk.Entry(window, width=10,textvariable=s2)
	txt2.grid(column=7, row=1)
	s5 = tk.StringVar()
	txt5 = tk.Entry(window, width=12,textvariable=s5)
	txt5.grid(column=7, row=2)
	lb6 = tk.Label(window, text="求解支路中有并联即输入1，无并联即输入0")
	lb6.grid(columnspan=6, row=2)
	lb21 = tk.Label(window, text="节点一")
	lb21.grid(column=0, row=3)
	s21 = tk.StringVar()
	txt21 = tk.Entry(window, width=10,textvariable=s21)
	txt21.grid(column=0, row=4)
	lb22 = tk.Label(window, text="节点二")
	lb22.grid(column=1, row=3)
	s22 = tk.StringVar()
	txt22 = tk.Entry(window, width=10,textvariable=s22)
	txt22.grid(column=1, row=4)	
	lb23 = tk.Label(window, text="元件个数")
	lb23.grid(column=2, row=3)
	s23 = tk.StringVar()
	txt23 = tk.Entry(window, width=10,textvariable=s23)
	txt23.grid(column=2, row=4)	

	def clicked():
		global n
		global m
		global tot
		global p1
		global p2
		global binglian
		n=(int(txt.get()))
		m=(int(txt2.get()))
		p1=int(txt21.get())
		p2=int(txt22.get())
		tot=int(txt23.get())
		binglian=int(txt5.get())
		window.destroy()
		window3()
		window2()
		
	btn = tk.Button(window, text="Continue", command=clicked)
	btn.grid(column=0, row=5)
	window.mainloop()


init()