电磁场第三次数值仿真报告

1、2012年秋季学期电磁场第三次数值仿真电磁场第三次数值仿真2012/12/9wy电磁场第三次数值仿真1. 题目叙述为了获得一定区域上的匀强磁场，可采用多组Helmholtz线圈结构。一种两对线圈的结构如图1所示。线圈半径a1，a2，线圈间距h1，h2，以及线圈中通过电流i1，i2可变化量，如图1（a）所示。为了定量衡量关注区域的磁场均压程度，过轴线做截面ABo1o2，取CD=0.8×AB和Eo1 = 0.8×o1o2，在CD和Eo1线段上每边均匀取20采样点，从而形成如图1（b）所示的采样节点，定义z方向B的不均压系数为：2.其中，为所有采样点的z方向磁感应强度平均值；为第

2、n个采样点的z方向磁感应强度值。N为采样点总数。如给定a1 = 1 m， i1=1 A，且h1 > 0.2 m， 情况下，如果要使得d达到最小，线圈半径a2，线圈间距h1，h2，以及线圈中通过电流i2应如何取值？(a)线圈结构示意图（b）磁场采样节点示意图2. 初步工作2.1. 推导和计算直线电流元在任意观察点的磁感应强度电磁场课上曾讲过相关推导过程与运算结果。柱坐标系下，对于下图：在观察点A出的磁感应强度为：在直角坐标系下，相应的公式是：对于以上两个公式，在直角坐标系下，有：cos4=AC·APACAP cos3=CA·CPCACP 相应的MATLAB代码如下：（记

3、得说只保留了Z方向分量）2.2. 积分法计算线圈组在其内任一点的磁感应强度要算线圈组在其内任意一点的磁感应强度B，必须先知道一个载流圆线圈在空间任一点产生的磁场情况。建立空间直角坐标系，作图如下：对于空间任意观察点P，有以下推导：rs=acos' ex+ asin' eyr=xex+zey x=rsinz=rcosR=r-rs=x-acos'ex-asin' ey+z ez R2=z2+a2+x2-2axcos'dl=ad'-sin'ex+cos'eyB=04Idl×RR3=0Ia402(-sin'ex+cos&

4、#39;ey)×(x-acos'ex-asin' ey+z ez)(z2+a2+x2-2axcos')32d'=0Ia402zcos'ex+zsin'ey+x-acos'ex(z2+a2+x2-2axcos')32d'进而得出：Bz=0Ia402a-xcos'(z2+a2+x2-2axcos')32d'相应matlab代码如下：*请简述你编的几个函数的功能每个函数都要说，每个函数用“*”隔开2.3. 数值法计算线圈组在其内任一点的磁感应强度一方面，考虑到用积分法计算磁感应强度，运算步数多，

5、运算时间长；另一方面，我们已经得出了直流电流元在任意观察点磁感应强度的表达式，基于相关表达式与函数，我们能够较为容易的得出正多边形线圈组在其内任一点的磁感应强度，进而用正多边形代替圆，在误差允许的范围内，进一步减少运算步数，缩短运算时间。在直角坐标系下，对于下图的正多边形（以正六边形代替）：它第i条边的起点坐标是（acosi-1,asini-1），终点坐标是（acosi,asini）。利用2.1中已经得出的calculate_b_line函数，我们循环算出每条边在观察点的B，并将其叠加，就得到边数为M的正多边形线圈组，在它们内部任意一点的磁感应强度。相应MATLAB代码如下：*function

6、 B=multi_solve_point(a,M,destinationpoint,h1,C) b=0;for i=1:M sita=2*pi/M; x1=a*cos(i-1)*sita); y1=a*sin(i-1)*sita); x2=a*cos(i)*sita); y2=a*sin(i)*sita); bn1=calculate_b_line(x1 y1 h1,x2 y2 h1,destinationpoint,C); bn2=calculate_b_line(x1 y1 -h1,x2 y2 -h1,destinationpoint,C); b=b+bn1+bn2;endB=b;*这种方

7、法的优点是在精度要求不太高的情况下，能用M较小的正多边形来代替圆形，计算时间较短。缺点是本方法精度较小，如果要提高精度，就必须要增加边数M，进而计算时间会提高。现用单个线圈轴线上的B来讨论M取不同的值时，运算精度和运算时间时间的变化。理论计算公式是B=0IR22R2+x232 ，来源于课件有关题目。运算结果如下：边数M10501001000数值结果B1/T5.6610e-0075.5268e-0075.5227e-0075.5213e-007理论结果B2/T5.5213e-0075.5213e-0075.5213e-0075.5213e-007相对误差2.5314%0.0991%0.0248%

8、0%运算时间 t/s0.0032790.0038570.0052060.031530上表的结果验证了数值法的优缺点。由表中数据可得，在边数M=50时，运算精度和运算时间都较好。因此，在之后的运算中，均用M=50的正多边形来代替圆。2.4. 不均压系数的计算。我们现在已经得到两个线圈中任一点磁感应强度B的计算公式，下面只需按照题目要求进行取点，再代入公式计算即可。l 利用对称法进行优化宏观的看，求解区域中磁场分布是对称的。由对称法的思想知如果我们只取上半部分的点，运算量会相应减少一半，而最终结果不变。证明如下：首先，在CD和Eo1线段上取点的个数均是偶数。令N是总点数400，则有：Bz=n=1n

9、=400BznN=2×n=1n=200BznN=n=1n=200BznN2=Bz1=1Nn=1n=400Bzn-Bz2Bz=2Nn=1n=200Bzn-Bz12Bz1=1N2n=1n=200Bzn-Bz12Bz1=1编写新函数fun_half用来计算不均压系数。程序思路比较简单，主体就是一个用来计算所有点B的for循环。但是在初次编程中出现了错误，原因是求解区域是随着h1的不同而变化的，最开始的时候我们把h1固定化了。MATLAB代码如下（只截取for循环和我们出错的地方）：* function d=fun_half(x)N=20;r=linspace(0,0.8,N);ini=0.

10、8*x(1)-1.6*x(1)/19*9;z=linspace(ini,0.8*x(1),N/2);B=zeros(N,N/2);for i=0:N-1 for j=0:N/2-1 b=multi_solve_point(1,50,-r(i+1) 0 z(j+1),x(1),1)+multi_solve_point(x(4),50,-r(i+1) 0 z(j+1),x(2),x(3); B(i+1,j+1)=b; endend*自此，我们已经得到用来计算的函数。下面的工作就是讨论，在h1,h2,a2,I2为何值时，取得最小的不均压系数。3. 采用数值法与fmincon进行优化考虑到实际情况下，

11、我们更倾向于在短时间优化出一个较为满意的结果。考虑计算时间因素，我们选用正多边形模型和fmincon算法对原问题进行优化。目标是在误差允许的范围内，用较短的时间优化出较好的结果。使用x,fv,ef,out,lag,grad,hess=fmincon(fun_half,x0,A,b,v1,v2)对原问题进行优化。3.3.1. 初步尝试将优化程序编好之后，我们随机选定了一些初值和约束条件进行尝试，但是结果并不理想。例如，当初值和约束条件为：x0=0.2 0.8 -2.5 1.1;v1=0.2 0.4 -3 0;v2=1 2 3 2;A=1 -1 0 0;b=0;迭代一定步数后，会出现中间变量mea

12、nB=NaN的情况，超出范围，无法计算。除此之外，某些初值和运输条件虽然可以得出结果，但是一方面迭代次数多，运算时间长；另一方面，得到的优化结果又偏大，猜想很可能是局部最优解，并不是我们要的全局最优解。综上所述，由于函数fun_half本身的复杂性和fmincon本身功能的局限性，我们必须先人为地缩小初值和约束条件的范围，使二者的选取更有方向性，更趋近于全局最优解的位置，从而缩短运算时间，提高运算精度，得到更好地结果。3.2. 对fmincon初值和约束条件的优化我们要优化的一共有4个量：h1,h2,I2和a2。如果一起讨论这4个参数，难度过大，而且，这4个参数对结果的“影响力”很有可能也是不

13、同的。很有可能某些参数对结果其主要约束作用，另一些参数对结果起次要约束作用。我们可以简化次要参数，并减小主要参数的取值范围，从而实现对fmincon的优化。l 对半径a1的讨论首先在h1,h2,I固定的情况下绘制a2曲线，结果如下：结果显示，其他三个参数取不同值时，的最小值都在a2=1.2左右取得。因此，我们能够认为a2对于的影响较为独立，在使用fmincon时，可以将初值和约束都选择在1.2附近。另外，在以下的讨论中，如不特殊指出，a2均取1.2。l 对电流I2的讨论用matlab做出电流I2变化时，I2的图像如下：根据结果分析，在其他参数固定的情况下，I2和的相互关系呈现如下特点：a) 在

14、I2为某个负值时有一个极大值出现，且取得极大值时对应的I2值随着其他参数的不同而变化。b) 其他参数固定的情况下，的最小值在I2=±3时取得，考虑fmincon迭代的需要，我们可以将I2的初值分别取在-2.8和2.8。c) 在h1=0.4，h2=0.6时，减小幅度较剧烈，猜想的最小值可能在h1=0.4，h2=0.6左右取得。l 对h1和h2的讨论因为此时只剩下h1和h2两个参数要讨论，参数数量已经减少很多；另一方面，相较于a2，I2间的区别，h1和h2的性质更为近似，他们可能会作为一个因子来影响，不方便单独讨论。所以，我们选择做h1-h2-的三维图来讨论这两个参数对结果的影响。结果显

15、示，关于h1和h2的变化关系较为复杂，但是通过使用matlab的图像旋转功能，能较为清晰的看出的最小值在h1=0.4，h2=0.5处取得。考虑肉眼观察的模糊性，利用fmincon做优化时，约束范围可稍微大一些。综上所述，我们最终选取初值和约束条件如下：x0=0.4 0.5±2.8 1.2v1=0.21 0.4-3 1v2=0.5 0.8 3 1.5A=1 1 0 0b=03.3. 用fmincon进行优化的结果考虑到用正多边形代替圆存在固有的误差，因此，我们选取边数分别为M=50和M=100的两种情况来进行优化。结果如下表所示：边数MI2的初值最优参数优化所得计算用时th1h2I2a

16、250+2.80.21000.80002.87041.06930.010431.35s50-2.80.48230.5717-3.00001.18880.005856.91s100+2.80.21000.80002.86971.06970.010357.09s100-2.80.48330.5728-3.00001.19020.005899.36s表中结果显示，当边数M=50之后，再增加边数，对最终优化所得的已经没有影响，所以可以认为在大多数情况下，均可用正50边形来代替原型。另一方面，优化所得的最优参数为h1=0.4823,h2=0.5728,I2=-3,a2=1.1888,计算用时也在可接受的

17、范围内。实际上，如果初值和约束条件进一步靠近最优参数处，计算用时能进一步缩小。经多次尝试，最小用时大概在20s左右。但是这样做的前提是我们已经知道最优解的精确位置，在实际操作中没有多大意义。4. 验证优化所得的结果对于原问题，由于数学关系非常复杂，难以用解析法精确的求出最优解，我们在以上的优化步骤中选用了作图法进行优化。但是，用作图法求解，本质是根据趋势找规律，是建立在“描述”的基础上，“逻辑推理”的成分较少，说服力不是很强。另一方面，由于函数fmincon本身的缺陷，只能求出局部最优解，在原函数存在震荡的情况下难以保证结果的正确性。因此，在优化完毕之后，我们还应用不同的方法来验证我们所得的优

18、化结果。4.4.1. 用作图法验证优化结果在用其他方法验证之前，我们最好用作图法来观察优化结果的情况，对其进行“毛估”，以防止发生方向性的错误，为其他的验证方法奠定基础。从不均压系数的数学表达式来看，这是一个用来表示磁场分布均匀性的参数。因此，我们可以做出线圈组间磁场分布的图像，通过判断磁场是否恒定来检验优化结果的正确性。相应matlab代码如下：*clc;clear all;tic;C=1.2;a=1.2; h1=linspace(0.2,2,40);h2=linspace(0.2,2,40); d1=zeros(40,40);for i=0:39 for j=0:39 d=fun_for_

19、test(h1(i+1),h2(j+1),C,a); d1(i+1,j+1)=d; endend mesh(h1,h2,d1)toc;*所得图像如下：由图观之，在XOZ平面上，磁感应强度B的分布已经相当均匀，仅在靠近线圈组附近存在一些变化。因此，我们可以认为所得的优化结果是正确的。可以进行其他方法的验证。值得注意的是，如果用遗传算法+积分法进行验证，在积分运算精度较高的情况下，运算时间过长；而在积分运算精度较低的情况下，其精确性甚至不如用正多边形模型来得高，因此，我们选择单独用这两个方法进行验证。4.2. 用遗传算法来验证优化结果考虑到函数fmincon本身只能求出局部最优解，为了消除其算法本

20、身可能对优化结果造成的误差，我们决定用遗传算法进行验证。为了消除自变量范围对最终结果的影响，我们遗传算法的自变量范围与fmincon中的相同。另，种群数量是100。运算结果如下：实际上，由于遗传算法是基于概率的一种算法，每次运算的结果都不一样，所以我们用其进行了多次运算，并取了最小结果如上图所示。用遗传算法每次用时大概在40分钟，所得最小结果是=0.00912。可见，这个结果大于用fmincon优化的出的结果，并不够优异；而且运算时间是用fmincon的几十倍。因此，从各方面来看，我们用fmincon优化出的结果都可以说是正确而“优秀”的。4.3. 用积分法进行验证由于我们实际上是用一个正50

21、边形来代替圆，与实际线圈的形状存在误差。为了减小这个误差，我们用积分法+fmincon进行验证。积分中的精度取默认精度1e-6。所得结果如下表所示：h1h2I2a2运算时间t正50边形0.48230.5717-3.00001.18880.005856.91s积分法0.44420.5420-2.80981.15862.9441e-0041968s相对误差8.57%5.48%6.77%2.61%/由表中结果显示，积分法的优化结果与正50边形的优化结果存在一定误差，下面分别来分析。首先对于不均压系数，两者差别很大，这可能是用多边形替代圆形的固有缺陷。实际上，多边形的边数M越少，不均压系数的值就越大，

22、而且，我们最终关注的其实不是，而是h1,h2,I2,a2这4个量该如何选取。从这个表中我们看出，用积分和用正多边形对这4个量的优化结果并不同，其中，h1的相对误差为8.4%，已经偏大。另外，把用积分法算的最优值带入到正多边形模型中，所得的不均压系数是0.0106，大于0.0058。综上所述，正多边形模型并不如积分模型精确，这就启示我们要对积分模型进行进一步探究。5. 对积分法的探究及优化5.5.1. 两种积分方法的比较经过一定的阅读和讨论，我们发现，所谓的正多边形法其实思路和梯形积分一样。因此，我们首先要讨论梯形积分公式和辛普森积分公示的差别。令函数为y=f(x),步长是h，积分区间是（a，b

23、），分割份数是n，梯形公式某区间上计算结果是Tn，辛普森公式同样区间上计算结果是Sn 。则有如下推导：h=(b-a)/n梯形求积公式：Tn=hk=1n-1fk+h2(f0+fn)令n=2m，则辛普森求积公式：Sn=h3(f0+f2m+4k=1m-1f2k+1+2k=1m-1f2k)由差值相关知识和积分中值定理可得：xxk+1fx-Txdx=f''(k)2xxk+1x-xkx-xk+1dx=-h312f''k, k(xk,xk+1)则梯形求积公式的误差：P1=abfxdx-Tn=k=0n-1xkxk-1fx-T(x)dxh312k=0n-1f''(

24、k)令：M2=maxf''(x)则：P1h212M2(b-a)即梯形公式是2阶收敛的。用类似的方法能够求出辛普森求积公式的误差为：P2h4180M4(b-a)即辛普森公式是4阶收敛的。而这次仿真作业里用的是自适应求积方法，步长会自动调整，精度会比一般的辛普森公式再高一些。因此，在分割次数较小的时候，用正多边形法（梯形积分法）和辛普森积分法误差会比较大。两者误差较大，则说明仿真所得磁场的分布有较大差别，用同样的算法计算存在误差在情理之中。解释了二者有误差之后，我们就应该选用精度更高的辛普森求积公式来做优化。由于运算精度和运算时间往往是一对矛盾的量，我们接下来就应讨论，在使用辛普森

25、求积公式运算时，如何做到运算精度和运算时间的最佳平衡。5.2. 讨论运用辛普森求积公式时的运算精度和运算时间l 单个线圈的情况l 用fmincon做优化时的情况由我们在4.3中的讨论可知，虽然辛普森积分法和正多边形法在h1，h2，I2，a2的优化结果上存在误差，但是二者都在3.3确定的范围内，且二者的均在初值附近，二者的差别相对于整个约束条件的范围也较小。因此，我们不需要改动fmincon的初值和约束条件。 结果如下表：绝对误差tolh1h2I2a2运算时间t/s1e-20.41770.5438-2.81071.18000.006222s0.5*1e-20.46600.5476-2.82181

26、.16660.003539s1e-40.50000.6029-2.81521.17764.7727e-004274s1e-60.44420.5420-2.80981.15862.9441e-0041968s从表中结果可见，除了tol=1e-2时的大于0.058之外，其余状态下的不均压系数均在0.058之内。而且，当tol0.5*1e-2时，实际的分割份数已经远大于50，精度应高于正50边形的情况。与此同时，运算时间为39s，小于正50边形的运算时间。因此，用辛普森公式优化的话，运算精度更高，运算时间更短，所以这是一个更为优秀的方案。经多次实验，绝对误差的范围选在0.5*1e-2到1e-3时，运

27、算精度和运算时间能取得较好的平衡。采用遗传算法对tol=0.5*1e-2的情况进行验证，所得结果如下：采用遗传算法所得比用fmincon算得的还要大，且用时很长，明显劣于用fminc优化。当然，也从一个侧面证明了用fmincon优化的正确性。l 对优化的总结至此，我们对本次任务能够做的优化已经全部完成。由于我们优化做的比较详细，最终的计算时间和计算精度均很优秀。现将我们所做的优化总结如下。a) 在算磁感应强度B的时候，我们只计算在Z方向上的分量，减少了一半以上的计算量。b) 在计算的时候，我们利用对称法，在保证结果绝对正确的情况下减少了一半计算量。c) 在对fmincon初值和约束条件进行优化

28、的时候，我们没有盲目的尝试，而是采用作图法估计出了初值和约束的大概范围，缩短了计算时间，提高了计算精度。d) 我们使用了正多边形模型和辛普森积分模型，并对二者进行了评价，最终选定了更为优秀的辛普森模型。之后，我们又对辛普森积分模型进行了讨论，找到了计算时间和计算精度都较为优秀的tol值范围。6. 对一些问题的进一步探究6.1. 尝试用解析法求最小不均压系数如果按照题中所给公式一步步推的话，公式过于复杂，难以计算。因此，我们必须要找到一个简化的模型进行求解。载流圆线圈在轴上的磁场B相对容易求，我们可以通过评估线圈组轴线上磁场的均匀程度来估计整体磁场的均匀程度。知道对于单个载流圆线圈而言，它轴线上

29、的磁感应强度是：B=0IR22R2+x232易知B是x的偶函数，那么它奇数次导数都是奇函数。因此，对于偶数个对称排列的线圈组而言，它们的奇数次导数全是0 。所以，评价其轴线上磁场的均匀程度可以通过讨论它们偶次导数来实现。如果线圈组轴线上磁场分布是理想均匀的，则必然有：Bxx=x0=2Bx2x=x0=limnnBxnx=x0=0但实际情况下往往不可能有无数个变量让你列出无穷阶导数，也就是说可控制的导数阶数是受变量个数约束的。对于题目中所给情况，电流I1是确定的，ai和hi之间有约束，只能看做两组变量，再加上电流I2，可以决定的等式个数是：1+2=3个。因此，我们需要列出本题4个线圈轴线上B的2阶、4阶、和6阶导数，选择不同的ai、hi和I2，让三个等式均为零即可。为了便于讨论，令：cos=Ra2+R2则相应的2、4、6阶导数表达式分别是：I1cos71-4sin21cos51R13+I2cos72-4sin22cos52R23=0I1cos111-12sin21cos91+8sin41cos71R15+I2cos112-12sin22cos92+8sin42cos72R25=0i=12Ii-5cos15i+120sin21cos131-240sin4icos11i+64sin6icos9iRi7=0下面就是解以上三个方程了。由于我们数学能力有限，解了