forn=1:
m;
fori=2:
14;
forj=2:
19;
b=u(i,j);
c=u(i,j+1);
d=u(i+1,j);
e=u(i-1,j);
f=u(i,j-1);
g=0.25.*(c+d+e+f);
u(i,j)=b+a.*(g-b);
end
end
end
mesh(u);
通过对相关资料的查询及学习,我认为运用MATLAB编程就是要将有限差分法的基本公式实现。
以此为出发点可以让解析和编写过程相对简单化。
首先,先要对电场的边界条件进行设定:
u=zeros(15,20);
i=2:
14;
u(i,20)=100;
使其满足题目的要求右侧边界电位为100,其余三边为0,由此也可以判断出此题目为第一类边界条件的问题。
同时,此处,也要特别注意,因为在MATLAB程序中,编写一个矩阵的默认顺序是从左到右,从上到下。
另外,命名矩阵时是先编写行,后编写列,因此,对于u(i,j)中,i代表的是第i行即对应纵坐标,j代表的是第j列即对应横坐标,与坐标表示有一定的差别,需要区分和注意。
之后是相对关键的编写部分,就是对电场内部的电位设定,此处,我通过分离变量法将u分解为u(x)和u(y),分别求解。
为让计算简便,我设定其内电场为线性变化的,即u(x)=kx+b,u(y)=0作为初值进行计算(应用分离变量法得到的结果是u(x)=Bsinhm∏/15
u(y)=Asinm∏/15。
可以通过此处在得出结果后进行对比和误差分析。
)又因为要在MATLAB中进行运算,所以首先要转化为i和j的形式,因此可得:
u(j)=kj+b和u(i)=0。
从而得到方程组:
u
(1)=k+b=0;
u(20)=20k+b=100。
解得:
u(j)=100/19*(j-1),即:
u(i,j)=100/19*(j-1)。
由此得到了内电场除去上下边界(衡为零)的电位方程,在MATLAB中我运用了for循环语句来实现其内电场的赋值,源代码如下:
forj=1:
20
fori=2:
14
u(i,j)=100/19*(j-1);
end
end
将电场各点的电位赋初值完成后,就需要带入公式进行计算了。
在之前先要对“加速收敛因子”α的最佳值和预定的最大允许误差进行确定。
关于加速因子的最佳值确定问题,我通过学习和查阅相关资料得到了两种方法,如下:
一是逐步搜索法;
将α的取值区间(1,2)进行M等分,α分别取1+1/M,1+2/M……1+(M-1)/M,通过上面提到的迭代公式2依次对同一个精度要求求出迭代次数k的值,并从中选出“加速收敛因子”α的最佳值,具体步骤如下:
步骤1给定等分数M和精度要求ε的值,令α的初始值为1;
步骤2令p=1,2,3,……M-1,重复步骤3-5;
步骤3αp=1+p/M;
步骤4按照如下公式迭代
找出符合精度要求ε的迭代次数kp;
步骤5比较找出kp值最小的αp为最佳的加速收敛因子α的值。
二是黄金分割法:
依据黄金分割法的思想,通过计算机自动选取最佳加速收敛因子α的近似值,具体步骤如下:
步骤1对(1,2)区间第一次0.618的分割,区间分界a1=1,b1=2,在(a1,b1)区间分割出黄金点point1=a1+0.618(b1-a1),进行迭代公式2的迭代,求出迭代次数K值,如果迭代次数没有超出规定的发散常数,迭代结束,否则转为步骤2。
步骤2在(1,1.618)和(1.618,2)之间进行第二次的黄金分割,找出分割点point2=a2+0.618(b2-a2),其中a2和b2是新分割区间的左右边界。
找出迭代次数最少的α。
以上两种方法是比较具有实践和研究价值进行的,在这里介绍和分享给大家探讨。
在我分析过后认为第二种较好,因为其可以通过计算机自动选取到最佳的加速收敛因子,但是实施起来比较复杂,要通过MATLAB或其它语言进行编码实验得到,我初步试验过,没能成功,就运用了第一种方法进行了计算,但是,我发现其运算过于复杂,笔算正确的成功率比较低,而运用电脑编程要重复多次运算,虽然,计算机运算快速且方便,但是,这样的源代码同样过于复杂,可以进一步进行深入的研究。
本次的题目解答中,由于是应用计算机计算,所以计算速度和计算效率在不同的α取值中没有明显区别。
因此,我就只是设定了一个1~2之间的一个数字:
1.8为“加速收敛因子”α的值了,进行了近似计算,实验误差在可接受的范围。
然后,是对于迭代次数的确定。
当迭代满足停止条件:
时,对应的迭代次数即为最终的迭代次数,这里的W为设定好的允许误差,为从简方式,我选择运用输入迭代次数,并通过多次输入对比输出结果选取最佳迭代次数的方法来选取,存在了一定的误差。
但是,相比于我运用上面的公式求解计算迭代次数,进而得出的图形的误差比要小许多,我一直尝试进行改进,但一直没能很好的解决误差过大的问题:
最主要的问题是边界值一直不满足我的设定值。
希望通过进一步的研讨和大家的交流能过得到最优方式。
通过以上的分析,我得出了进一步的编程代码:
a=input('pleaseinputa(1m=input('pleaseinputm(1其中a和m分别为最佳“加速收敛因子”和迭代次数,均为我们输入数值计算的方式。
先要预算出一些取值范围,然后通过分别输入,对比输出结果进而来选取最优解。
之后的源代码的编写思想就是通过迭代计算公式来实现对电场内部各点电位值的计算,因此,在输入“加速收敛因子”和“迭代次数”的数值后,我通过for的循环语句实现了对电场内部的电位迭代求解,并用图形形式表示出来。
其源代码如下:
forn=1:
m;
fori=2:
14;
forj=2:
19;
b=u(i,j);
c=u(i,j+1);
d=u(i+1,j);
e=u(i-1,j);
f=u(i,j-1);
g=0.25.*(c+d+e+f);
u(i,j)=b+a.*(g-b);
end
end
end
mesh(u);
上述代码中,b,c,d,e,f分别场内一点及其左右上下各点的电位值,之后的g和u(i,j)则实现了逐次超松弛法的迭代,最后用mesh实现图形的输出。
六、结果和误差的对比分析
由分离变量法得到的初步结果是u(x)=Bsinhm∏/15和u(y)=Asinm∏/15,从而应有u=Csinh(m∏/15)sinm(∏/15)的结果,从其表达式的形式,我们可以简单推断出电场的电位在x轴上的分布是成虚指数函数的形式呈现,而在y轴上是成正弦函数呈现的,这和通过有限差分法的计算机仿真的图形基本符合,如下图:
因此,可以基本确定通过计算机仿真得到的结果的正确性。
另外,可能存在较大误差的方面是跌代次数,因为迭代次数是我们人工输入的,所以判断其正确与否要通过输出的结果来判定、来选择,下面几幅图是当“加速收敛因子”α=1.8时,不同的迭代次数m的对应图形:
m=1时:
m=10时:
m=20
m=30
m=32
m=35
m=40
通过上面迭代次数m取不同的值对应的不同图形分析可以看出,在m=1,10,20时期计算结果都存在较大误差,图形与分离变量法得到的结果相差很大,有明显的不规则部分,不成现光滑的曲线。
而当m=30后图形与分离变量法计算到的图形近似,而且图形规则,曲线光滑,而且随着m值的增大图形的变化非常微小了,因此,我选择m=32为迭代次数,所得到的结果误差比较小,结果相对准确。
七、研究结论
学习了运用有限差分法计算第一类边界条件的电位值,通过计算机仿真进行了实际例题的求解,通过与分离变量法求解的对比,证明可行性和有效性,适用于解决第一类边界问题,而且更加方便和快捷。
同时,“加速收敛因子”最佳解问题比较复杂,不过通过计算机仿真的过程中,不同数值“加速收敛因子”的加速效果不易体现出来。
跌代次数对计算结果又较大影响,迭代次数越多,结果越准确,误差越小,但相应的计算也会更加复杂,因此只要迭代次数得到得结果满足一定的允许误差范围即可。
附录:
1三类边界条件:
第一类:
已知电位在边界上的数值。
第二类:
已知电位的导数在边界上的数值。
第三类:
已知电位在一部分边界上的数值和在其余边界上的导数的数值。
2解析法包括:
分离变量法,镜像法,复变函数法,保角变换法,格林函数法。
3数值法包括:
数值积分法,有限差分法,有限元法,矩量法。
4解析法的优点:
由有限个项构成闭合解或无穷级数,通常具有鲜明的物理意义。
缺点:
解题范围窄小,对边界形状十分挑剔。
数值法的优点:
解题范围宽广,对边界的形状没有限制。
缺点:
数据离散,物理意义深藏其中。
心得体会:
有限差分法的公式套用很简单,记住结果很简单,但是要想深刻理解其内容和灵活应用则比较困难,需要一定时间的深刻学习和做题应用,同时,这个方法确实是一个非常方便和有效地解决三类边界条件问题的工具。
MATLAB在电磁场方面的应用相当广泛,其内拥有的很多函数都可以模拟静电场中的实例,就比如这次的研讨就应用了矩阵,通过有限差分法来模拟静电场的电位分布,相信将来还会有更多的有意思,有研究价值计算和模拟电磁场的相关知识。
是学好电磁场必须要牢固掌握的工具。