《仿真技术在电气工程中领域的应用》课程作业.docx

《《仿真技术在电气工程中领域的应用》课程作业.docx》由会员分享，可在线阅读，更多相关《《仿真技术在电气工程中领域的应用》课程作业.docx（38页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、仿真技术在电气工程中领域的应用课程作业（2019-2020 学年第一学期）作业0202稀疏线性方程组的求解作业0203微分方程的求解作业0302多机系统短路故障后时域仿真学院：电力学院任课老师：武志刚作业0203微分方程的求解一、问题重述求解具有如下初值的微分方程粤=侬)x ea,bil):Bi2U=Bi2U-Bi2il*UilUif(i2=il+l):Di2i2 = Bi2i2if (i2il):Li2il=Bi2il/Dililreturn LDUdef LDU_store(A): #三角存储检索，不是实际的坐标L,D,U=LDU(A)size=A.shapeOL_value=L_row=

2、L_column=D_value=U_value=U_row=U_column=for i in range(l,size):for j in range(i):L_d(LiUDL_d(i)L_d(j)for i in range(size-l):for j in range(i+l,size):if(UiU!=0)：U_d(UiUDU_d(i)U_d(j)for i in range(size):D_d(Dii)return (L_value,L_row/L_column),(D_value),(U_value/U_row/U_column)def front(L,b): #前代，用LDU系

3、数存储size=b.shapeOz=np.zeros(size,l)for i in range(len(b):sum = 0if Ll.count(i)!=0:for j in range(Ll.index(i),Ll.index(i)+Ll.count(i):sum+=L0j*zL2j,0zi,0=bi-sumreturn zdeffrontl(L,b): #前代，没用LDU系数存储size = b.shapeOz=np.zeros(size,l)sum_list=for i in range(len(b):sum=0for j in range(i):if(i0):sum+=Lij*zj

4、,Osum_d(sum)zi,0 = bi,O-sum_listireturn z def mid(D,z): #规格化，用LDU系数存储size = z.shape0y = np.zeros(size,l)for i in range(len(D):yi,0=zi,0/Direturn ydef midl(D,z): #规格化，没用LDU系数存储size = z.shape0y = np.zeros(size,l)for i in range(D.shapeO):yi,0=zi,0/Diireturn ydef back(U,y): #回代，用LDU系数存储size=y.shapeOx=np

5、.zeros(size,l)for i in range(size):sum = 0if Ul.count(size-i-l)!=0:forjinrange(Ul.index(size-i-l),Ul.index(size-i-l)+Ul.count(size-i-l):sum+=U0j*xU2jL0xsize-i-l,O=ysize-i-l,O-sum #现在 x 是从 xsize-l到 x0x_reverse=np.zeros(size,l)for i in range(size):x_reversei,O=xsize-l-i,Oreturn x_reversedef backl(U,y)

6、: #回代，没用LDU系数存储size=y.shapeOx = np.zeros(size, 1)sum_list=for i in range(U.shapeO):sum=0for j in range(i):if(i0):sum+=Usize-i-lsize-j-l*xsize-j-l,Osum_d(sum)xsize-i-l,0 = ysize-i-l,O-sumJistix_re verse = np.zeros(size, 1)for i in range(size):x_reversei,O = xsize -1 - i,0return x_reversedef solve(A,b

7、): #求解 AX=bL/D,U=LDU_store(A)X=back(U,mid(D,front(L,b)return X def solve_Sparse_vector(A,b,xJndex): #稀疏向量法，xjndex 为所需求解的特定 变量L,D,U=LDU(A)y=midl(D, frontl(L, b)list_delete_x=for i in range(x_indexjen(b):if not any(Ux_index-l:i,i): #若第 i 个 元素和 xjndex 无关，则 Ux_index-l:ii应全为 0;list_delete_x.append(i)for

8、i in list delete x:U = np.delete(U, i, axis=O)U = np.delete(U, i, axis=l)y = np.delete(y, i, axis=O)for i in range(xjndex-l):U = np.delete(U, 0, axis=0)#每次循环都删第一行第一列，循环后则把原矩阵的前k行k列都删去了U = np.delete(U, 0, axis=l)y = np.delete(yz 0, axis=0)x=backl(U,y)return x0A=np.array(2/0,0/0/0/0,l,0,0,0,0,l,10,2,0

9、,0,1,0,0,0,1,0,0,0,0,0,2.,0,0,1,0,0,1,0,0,0,10,0,0,2,0,0,1,0,0,1,0,0,04,0,0,2,0,0,0,0,0,0,1,0,04,0,0,2,0,0,0,1,0,0,11,0,0,1,0,0,2,0,0,0,0,0,o,o,o,o,o,ozoz2,i,o,i,o,0,1,1,000,0,1,2,0,0,0,0,0,0,1,0,100,0,2,1,0,0,0,0,0,000,1,0,1,2,1,1,0,0,0,1,0,0,0,0,04,2)b=np.array(4U4U4U4L4U4U4U4U5U5U545)print(解得的X为：

10、。print(solve(A,b)print(需要求解的特定变量的值为：)print(solve_Sparse_vector(A,b,3)#solve_Sparse_vector函数的第三个参数则为需求解的某个特定变量的下标；如需 求x3,则设置为3二、作业203python代码import numpy as npimport t as pitimport sympy as syimport pandas as pd#定义欧拉法、隐式梯形法、改进欧拉法的函数def calculate_f_xy(derivative, x=None, y=None):”根据f(x,y)参数个数采用不同赋值”，i

11、f x and y:return Noneelif x in derivative.code.co_varnames:dy_dx = derivative(x=x)elif y in derivative.code.co_varnames:dy_dx = derivative(y=y)else:dy_dx = derivative(x=x, y=y)def base_euler(derivative, n=50, limit=(0,1), begin=(0,1):III欧拉法求微分方程Parameters(derivative 导数 dy/dxn步长limit 函数区间(start, end)

12、begin 初值(xO, yO) instep = (limitl - limit0) / nx = np.linspace(limitO, limitl, n)y = np.full(n, beginl, dtype=np.float32)for i in range(l, n):yi = yi -1 + step * derivative(yi -1) return x, ydef improve_euler(derivative, n=50, limit=(0, 1), begin=(0,1):III改进欧拉法求微分方程Parametersderivative 导数 dy/dxn步长li

13、mit 函数区间(start, end)begin 初值(xO, yO)instep = (limitl - limit0) / nx = np.linspace(limitO, limitl, n)y = np.full(n, beginl, dtype=np.float32)for i in range(l, n):yi = yi -1 + step / 2 * (derivative(yi -1) +derivative(yi -1+ step *derivative(yi -1)return x, y def implicit_trapezoidal(derivative, deriv

14、ative_again, n=50, limit=(0/ 1), begin=(0, 1), e=0.00001):in隐式梯形法求微分方程Parametersderivative 导数 dy/dxderivative_again d(dy/dx)/dyn步长limit 函数区间(start, end)begin 初值(xO, yO)e收敛判据|y-y| eIllstep = (limitl - limitO) / nx = np.linspace(limitO, limitl, n)y = np.full(n, beginfl, dtype=np.float32)for i in range

15、(l, n):nocoverage = Truewhile(nocoverage):F_y = yi -yi-l - step / 2 * (derivative(yi -1) + derivative(yi)dF_dy = 1 - step / 2 * derivative_again(yi) yi -= F_y/dF_dyF_y_iter = yi -yi-1 - step/2 *(derivative(yi -1) + derivative(yi)nocoverage = False if abs(F_yJter - F_y) e else Truereturn x, y#定义输入的导数

16、、步长、函数定义域、原函数#sympy定义x、y数学符号 x, y = sy.symbols(x y)# 用 x、y 定义 f(x,y) deriv = -20 * yderiv_2 = sy.lambdify(y, sy.diff(deriv, y), numpy)deriv = sy.lambdify(y, deriv, numpy)ESS = pd.DataFrame(columns=Euler, Improve Euler, Implicit Trapezoidal).name = N#图形绘制fig, axes = ots(2, 2, figsize=(20,10)N = 10, 1

17、5, 20, 100# titles = $Steps is 100$, $Steps is 50$ $Steps is 20$ $Steps is 10$for ax, n in zip(en()z N):n, limit = n, (0,1)x_r = np.linspace(limitO, limitl, n)y_r = sy.lambdify(x, np.e*(-20 * x), numpy)(x_r)x_b, y_b = base_euler(deriv, limit=limit, n=n)x_i, y_i = improve_euler(deriv, limit=limit, n=

18、n)x_it, y_it = implicit_trapezoidal(deriv, deriv_2, limit=limit, n=n)ax.plot(x_r, y_r/-g。label=original)ax.plot(x_b, y_b, -b, label=euler)ax.plot(x_i, y_i, -.r, label=improve euler)ax.plot(x_itz yjt, L.y) label=implicit trapezoidal)ax.set_title(fStep is l/n:.3f, fontsize=20)ax.legend(loc=best, prop=

19、size: 16)#三种方法的误差平方和str(n), Euler = sum(y_r - y_b)*2)str(n), Improve Euler = sum(y_r - y_i)*2)str(n), Implicit Trapezoidal = sum(y_r - y_it)*2)layout()0# ig(differential_solve_, bbox_inches=tight,)三、作业302python主函数代码import numpy as npimport pandas as pdimport sympy as syimport timefrom g import solve

20、import t as pitfromimpedanceimportZ_value,alpha,branch_data_for_imp_SC,branch_data_forJmp_afterCutimport Ydef subs(func, value):”将表达式用value字典替换“if type(func) in (int,float):return funcelif callable(func):return func(*s()else:return (value)def Pe(E=None, Delta=None, Z_value=None, Alpha=None):三机系统Pe方程

21、Pe_list =for i in range(3):j,k = 0,1,2-ipe=Ei*2/Z_valuei,i*sy.sin(Alphai,i)+Ei*Ej/Z_valueij*sy.sin(Deltai-Deltaj)*sy. pi/180-Alphaij)+Ei*Ek/Z_valuei,k*sy.sin(Deltai-Deltak)*sy.pi/180-Alphai,k )Pe_d(pe)return Pejistdef newton_laphson(Functions, init_valuesz e=0.00001):III牛顿拉夫逊法parametersfunctions多个形如

22、f(x,y,z,.)=0的待迭代函数的迭代类型(向量)init_values函数待求解变量的初始值，其中元素个数与所有函数的总变量个 数相同h步长e 收敛满足 max(|F(n+D-Fn|) eIII#初始化雅可比矩阵iter_keys, iter_values = np.array(list(init_(),np.array(list(init_s()Jacobian = sy.diff(Functions, iter_keys).transpose()Functions = sy.lambdify(all_var, Functions/numpy)nocoverage = True作业02

23、02稀疏线性方程组的求解一、问题重现利用课堂上介绍的稀疏线性方程组求解方法编写通用程序，求解下述方程组。Ax = b其中:1 0 00 0 10 0 11 0 00 0 00 0 02 0 00 2 10 1 20 0 00 1 00 0 0OOP0 0 00 0 01 0 00 0 11 0 00 0 00 1 00 0 02 1 01 2 10 1 2) =(4 4444444555 5)丁二、解题思路1.基本作业要求由于本文需解决的问题是：编写求解稀疏线性方程组的通用程序，因此可采 用基于因子表分解的求解线性方程组的方法。假设需要求解的方程组为：Ax=b又因为每个方阵都可以进行LDU分解，因此矩阵A可分解为下三角矩阵L、while(nocoverage):#前一次所有faunction的值构成的arrayFunctions_values = np.array(Functions(*iter_values)Jacobian_values = (dict(zip(iter_keys,iter_values)#将雅可比矩阵转成np.float32类型，否则无法求解Jacobian_values = np.array(Jacobian_rix(),dtype=np.float32)it