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

仿真技术在电气工程中领域的应用课程作业（2019-2020 学年第一学期）作业0202稀疏线性方程组的求解作业0203微分方程的求解作业0302多机系统短路故障后时域仿真学院：电力学院任课老师：武志刚作业0203微分方程的求解一、问题重述求解具有如下初值的微分方程粤=侬)x ea,b< dxy(%0)=刈方程中f(y)不显含x,求解的结果为x与y在定义域a,b上的数值序列二、程序涉及的微分方程求解的流程图StartEnd三、程序运行后的结果对比求解的具体微分方程为:= -20x x g0,1 dxy（o）= i当步长Ax ='中口分别取10, 15, 20, 100对应的欧拉法、改进欧拉法、隐 n式梯形法的图像分别如下：.当步长为0.1时(n=10)欧拉法求得的图像上下波动，原因是f(y)=20y与y负相关，过大的步长产 生较大的误差，甚至使与其叠加后的yi为负值，于是f(yi)变为正值，叠加后导 致y2又变为正值，如此往复，甚至形成一个无法收敛的波动。改进欧拉法为一条与x轴平行的直线，其原因在于yo=l , f(yo)+f(yo+O.l-f(yo)=O,则叠加后 yi =yo= 1,于是根据数学归纳法，yn=yn-i=.=yo= 1, 即改进欧拉法的图像为一条直线。隐式梯形法由于每一轮循环都要迭代至该轮收敛，因而在step=O. 1时就能很 快收敛至原函数。1 .当步长为0.067时(n=15)欧拉法的图像为衰减地收敛到与真实曲线重合，而改进欧拉法较真实解偏大, 隐式梯形法的曲线依然与真实曲线高度吻合。3当步长为0.05时(n=20)欧拉法图像在yi后就保持为0值，原因是yi=yo+O.O5*f(yo)=O,则 y2=y 1 +0.05*f(y 1 )=0,由数学归纳法 yn=yn-1=.=yi =0改进欧拉法依然较真实解偏大。4当步长为0.01时(n=50)三种方法的解皆与真实解贴近，且三者的精确度为8排序为£隐£改£欧。四、微分方程求解方法的数值评价由于方程为非线性方程，不能选取拟合优度R2作为评价指标，这里直接选取误 差平方和ESS作为评价指标，ESS指标如下：几ess = Z(B %)21=1200.0.0.1000.0.5.57732E-05注：蓝色为低精度，粉红色为高精度由色阶图，欧拉法和改进欧拉法在步长很大时的精度远低于隐式梯形法，随 着步长减小，欧拉法和隐式梯形法的精度均有明显提升，而在步长变化过程中隐 式梯形法始终能保持很高的精度。作业0302多机系统短路故障后时域仿真一、问题重述给定初始状态稳定和元件参数已知的多机电力系统，对其短路故障和故障切 除后的状态进行时域仿真。二、问题建模.主要模型研究电力系统运行的稳定性，主要是确定各同步电机是否处于同步运行状态, 即各转子间有无相对摇摆。功角S随时间变化描述了各发电机转子间的相对运动, 恰好反映了各发电机之间是否同步。根据算例所给的条件特点，考虑使用同步电机转子二阶模型，即警=（4一久）/巧由隐式梯形法得，tn至31时步的差分方程为“=0i+l -一“=0i+l -一(cuyi+i + (Un)0盗 (rn,n+l K,n+1) + (Bi,n - e,n)=。对于三机系统，任意一台发电机的电磁功率如下Pei = sin an +叔IPei = sin an +叔I而 j sin (如otij) + / j sin (左a语)其中= & - %觎=&一淳&jk = 5j 6k1 .方程组中未知参数的求解本小节主要针对电磁功率Pe表达式中的未知参数进行计算。其中，未知参 数主要为各台发电机的输入阻抗和之间的转移阻抗，Pe表达式中需要用到上述阻 抗的模值及阻抗角。计算输入阻抗和转移阻抗的关键是生成节点阻抗矩阵，根据计算公式即能算 得所需阻抗值。为了正确生成故障发生时和故障切除后的节点阻抗矩阵，需要按 以下步骤进行操作：（1）定义通用的数据格式（2）生成正负零序的节点导纳矩阵（3）计算附加阻抗（4）生成节点阻抗矩阵，计算发电机节点的输入阻抗和之间的转移阻抗具体程序及实现方法将在第三节进行详细说明。三、程序流程本程序的实现流程及具体实现方法将在本节详细介绍。1.主程序主程序，即利用牛顿-拉夫逊法及隐式梯形法求解模型的流程图，如下图所 示：t = t + h（h为步长）更换Pe表达式中的参数 (阻抗模值及阻抗角)结束仿真Pe表达式中的参数不变(阻抗模值及阻抗角)更换雅可比矩阵雅可比矩阵形式不变利用牛顿拉夫逊法 进行步长内的循环 求解下一时刻的值利用牛顿拉夫逊法 进行步长内的循环 求解下一时刻的值该部分由power_system_文件内程序实现。(1)接口函数 newton_laphson(Functions, init_values, e=0.00001)牛顿拉夫逊法的接口函数，其参数如下：Functions需要通过牛拉法迭代的形如f(x,y,.,z)=0的一组代数方程-init_values初始值的字典-e前后两次迭代的差值，所有方程组一次迭代前后差值最大者e该函数将传入的方程组迭代收敛至给定误差后，返回迭代后值的字典。 hiding_trapezium(derivative, n=50, limit=(O,l), e=0.0001)隐式梯形法的接口函数，其参数如下： derivative含所有方程(包括微分方程和代数方程)及相应初值的字典n绘制点个数 limit绘图定义域e传入牛拉法中的迭代精度该函数将微分方程先变形为代数方程，并与代数方程一同传入牛拉法函数迭 代。 Pe(E=None, Delta=None, Z_value=None, Alpha=None)E算得得三台同步电机电势向量-Delta三台电机功角sympy符号-Z_value每一次系统状态改变后，相应的发电机组输入、转移阻抗模值 矩阵。-Alpha每一次系统状态改变后，相应的弧度制阻抗角矩阵返回值为三台电机的电磁功率向量。(2)程序编写思路A.隐式梯形法a.生成定义域内对应步长的numpy序列；b.将所有方程中属于微分方程的定义转换为隐式梯形法对应的代数方程；c.利用给定初始值生成每一时间因变量的numpy序列；d.将每一次牛拉法迭代后的初值和代数方程一同传入牛拉法迭代，产生下一时 间点数据，修改因变量序列对应时间索引的值。B.牛顿拉夫逊法为了发挥numpy、sympy数值计算的最大优势，使用sympy自动求导功能, 并且成功尝试使用sy.diff函数，让代数n个方程组成的向量对由n个自变量组成 的向量求导，从而能一步产生n*n的求导后矩阵，大大简化了程序编写，具体过 程如下：a.初始化雅可比矩阵及各参数；b.进入迭代过程；c.将每一次迭代初值代入代数方程；d.将每一次迭代初值代入雅可比矩阵，求解代数方程组计算因变量变化量；e.利用修改每一次迭代值，重复步骤3；f.判断两次的代数方程差是否小于给定误差，若小于则返回最终迭代值，否则 回到过程3。2.计算Pe表达式参数流程:该部分由Ypy和文件中的程序实现。具体程序实现如下所示。（1）定义通用的数据格式本程序定义了通用的数据格式，使得对于任何满足格式的案例，不需要修改 程序代码就能得到计算结果。本程序定义了三个数据文件：描述节点参数的Node_、描述支路参数的branch.、描述变压器参数的trans_o它们定义的皆为numpy的array格式，例题 的数据保存在附件中的data文件夹中。具体定义格式分别如下所示：Node_的数据格式:组员信息及分工组员：阮艺源 胡嘉铭周欣缘 汪笑雨阮艺源17级电气工程卓越班.问题1部分Python程序编写1 .问题2基于欧拉法、改进欧拉法、隐式梯形法的微分方程求解的Python程 序编写、绘图、中文文档写作.问题3通用牛顿拉夫逊法、隐式梯形法接口的Python程序编写、绘图、相 应部分文档写作胡嘉铭17级电气工程卓越班.问题1部分Python程序编写、中文文档写作1 .问题2讨论并共同设计编程思路.问题3输入阻抗、转移阻抗的通用计算（包括Y.py和文件）、通用数据格式 接口的设计及程序编写、相应部分中文文档写作周欣缘17级电气工程中澳班.问题1 MATLAB程序版本1 .问题2 MATLAB程序版本.问题3 MATLAB程序版本2 .英文版本的论文写作汪笑雨16级电气工程五班1 .问题1讨论程序编写思路、修改论文.问题2讨论程序编写思路、修改论文2 .问题3节点导纳矩阵的程序编写、修改论文列数含义1节点序号2节点类型（短路节点记为0,负荷节点或联络节 点记为1,电源节点记为2）3电源节点电动势E，4电源节点的XJ （正序）5电源节点X2 （负序）6电源节点Xo （零序）（1）列3-6：不是电源节点的节点标为。（2）本题中，节点0：地；节点1： A；节点2： B；节点3： C；节点4： D；节点5： E1；节点6： E2；节点7： E3。（3）电源节点定义是在XJ外。branch_的数据格式:列数含义1起始节点编号2终止节点编号3切除故障前支路上等价的电阻4切除故障前支路上等价的电抗5切除故障后支路上等价的电阻6切除故障后支路上等价的电抗7零序电路中的电抗(1)若终止节点为地(负荷)，则将其置零。(2)将短路节点定为节点lo(3)原始数据中不需要加入Xd'。trans_的数据格式:列数含义1变压器两端节点中靠近短路点的节点序号2变压器两端节点中远离短路点的节点序号3变压器在列1对应一端的绕组接法4变压器在列2对应一端的绕组接法(YN： 0, Y:l, D:2)(2)生成节点导纳矩阵利用支路追加法，根据支路参数的数据可生成节点导纳矩阵。该部分程序在Y.py文件中。文件中定义了两个函数，Y_generation()和 Y_generation_ZS(),前者可以生成正负序的节点导纳矩阵，后者可生成零序的节 点导纳矩阵。由于正负序电路中除XJ外的参数是相同的，而零序电路中线路电 抗与正负序电路不同，故设计两个函数分别计算。(3)计算附加阻抗根据正序等效定则，不对称故障分析时，可将负序、零序电路的效果等效为 附加阻抗叠加在正序电路上，使得分析不对称故障时只需分析扩展正序电路。I .负序、零序短路点处的输入阻抗向branch_data中增加发电机节点的Xd2,及改变负荷等效阻抗值为负序电抗 值，即能生成新的temp_branch o将temp_branch传入节点导纳矩阵生成函数 Y_generation()中，能够生成负序电路的节点导纳矩阵丫2。将丫2求逆即可得到 节点阻抗矩阵Z2,从中便能读出负序电路中短路点处的输入阻抗Z2ff o零序电路的求解方法类似，不同之处为：零序电路需要考虑零序电流的流通 情况，此时则需要变压器参数的数据trans_data。对于两侧不都是星型接地的变 压器，零序电流不能流经其下游电路。并且，零序电路中不考虑负荷。程序中，通过NSI()和ZSI()分别计算得到负序、零序短路点处的输入阻抗。II .附加阻抗的生成根据正序等效定则，对于不同的不对称故障类型，附加阻抗的计算方式不同0 本程序由函数X_star()生成附加阻抗X型其第四个输入参数为短路类型，输入 不同的数值则能计算得到不同的。其定义如下所示：短路类型对应的第四个输入参数附加阻抗的计算三相短路10两相接地短路2X2ff / X()ff两相短路3X2ff单相短路4X2ff + Xoff(4)计算发电机节点的输入阻抗和之间的转移阻抗此步需要计算出节点阻抗矩阵，并由此得到输入阻抗和转移阻抗。生成节点 阻抗矩阵前，首先必须得到对应的正确的支路数据。本程序通过函数branch_data_for_imp_SC(),向branch_data中增加发电机节 点的Xd2,及短路节点处的X2生成故障发生时的新支路数据。利用函数 Impedance()调用新支路数据，生成发电机节点的输入阻抗和之间的转移阻抗。 输出结果为mXm矩阵Z,对角元为输入阻抗，非对角元为转移阻抗(m为发电 机节点的个数)。通过函数 branch_data_for_imp_afterCut(),向 branch_data 中增加发电机节点 的Xd2,并更新切除故障处附近的支路数据，生成故障切除后的新支路数据。类似 上段做法，生成故障切除后的发电机节点输入阻抗和转移阻抗组成的mXm矩阵。通过函数Z_value()和alpha()，可以生成上述矩阵Z各元素的模值和阻抗角， 将其代入模型的Pe表达式中。四、程序求解的算例G 1 0.1, X2 0.1,7j = 10sG-2 :xfd = 0.15/2 = 0.15,2j = 7sG 3 :% = 0.06m2 = 0.06,= 15s变压器电抗：xT1 = 0.08, xT2 = 0.1, xT3 = 0.04, xT4 = 0.05；线路电抗:AB 段双回Xli = 0.2, xL0 = 3.5xL1；BC 段双回Xli = 0.1, xL0 = 3.5xL1o 系统的初始状态：VD0 = l.o, SLD0 = 5.5 + jl.25, S20 = 1.0 + j0.5, S30 = 3.0 + j0.8。扰动事件描述：在线路AB段首端f点发生两相短路接地，经0.1s切除故障线路。五、程序运行结果及分析.故障各阶段仿真时间的分配仿真过程中0-0.05S为稳定状态，0.05sf点发生两相接地短路，故障持续至 0.15s, 0.15s后故障切除。1 .仿真0.5s后各状态量，及同步电机相对功角变化曲线图由曲线图1,可得出以下结论：(1)由系统中各同步电机相对功角变化，可知该系统在电机第一摆能保持 稳定，而未发散。(2)三台机组中，仅有电机一在故障切除，即0.15s左右速度达到第一摆 动最大，电机二、电机三分别在0.44s及0.5s后才至第一摆动最大值。图1仿真0.5s同步电机各状态量及相对功角变化.仿真3s后各状态量，及同步电机相对功角变化曲线图为了探究各台机组在故障切除后最终能否达到相对稳定，讲仿真时间延长至 3s,从图2有如下新的发现：（1）三台同步电机在故障切除后都无法达到相对稳定状态，而是在上下反 复波动。原因是模型中没有考虑阻尼项导致摆动无法衰减。（2）三台机组的转速是波动上升的，说明转速变化同时具有波动分量和升 高分量。图2仿真3s同步电机各状态量及相对功角变化六、程序创新点.本程序定义了通用的数据格式，使得对于任何满足格式的案例，不需要修改程 序代码就能得到计算结果。1 .本程序利用了 sympy符号库的自动求导功能与优化的latex输出，免去了求解 雅可比矩阵时繁重的求导任务。2 .本程序编写了许多通用的函数，如：牛顿-拉夫逊法、生成节点导纳矩阵、求 解输入阻抗及转移阻抗等，日后针对其他问题均能调用这些函数进行求解。一、作业202python代码import numpy as np def LDU(A): #LDU 分解函数size=A.shapeOB=A.copy() #直接的=会将地址赋值，要重新建立一个矩阵，要用copy()L= np.zeros(size, size)D= np.zeros(sizez size)U= np.zeros(size, size)D00=B00for i in range(size):Lii = lfor il in range(size):for i2 in range(ilsize):for j in range(il+l,size):if(i2=il):Ui2U=Bi2U/Di2i2if(i2>il):Bi2U=Bi2U-Bi2il*UilUif(i2=il+l):Di2i2 = Bi2i2if (i2>il):Li2il=Bi2il/Dililreturn LDUdef LDU_store(A): #三角存储检索，不是实际的坐标L,D,U=LDU(A)size=A.shapeOL_value=L_row=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系数存储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(i>0):sum+=Lij*zj,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.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): #回代，没用LDU系数存储size=y.shapeOx = np.zeros(size, 1)sum_list=for i in range(U.shapeO):sum=0for j in range(i):if(i>0):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): #求解 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 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,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为：。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)参数个数采用不同赋值”，if 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)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步长limit 函数区间(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, derivative_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(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#定义输入的导数、步长、函数定义域、原函数#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, 15, 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=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(f"Step is l/n:.3f", fontsize=20)ax.legend(loc='best', prop='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 solveimport 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方程"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多个形如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作业0202稀疏线性方程组的求解一、问题重现利用课堂上介绍的稀疏线性方程组求解方法编写通用程序，求解下述方程组。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