作业十Potential(matplotlib, scipy)

2013301020084 许晗
作业 Gauss-Seildel方法 SOR方法

---

摘要

这次作业计算部分用Fortran完成, 估计python店scipy或者numpy包里有处理数组的方法; 画图部分用python完成. 解决了二维方形箱(V=0)里, 两根电势恒定的杆子(势分别为-1和+1)情况下电势的分布, 并且探究了Gauss-Seidel和SOR计算方法的与系统大小(L+1)的关系.

---

正文

计算原理

边界条件

$V(x=\pm L,y=\pm L)=0$
电势恒定的杆子造成: $V(x=\pm 3, -3<y<+3)=0$

Gauss-Seidel方法

计算格点$(i,j)$的电势:

$$V_{new}(i,j)=\frac{1}{4}[V(i-1,j)+V(i+1,j)+V(i,j-1)+V(i,j+1)]$$
每次计算的V(i,j)要用到计算V(i+1,j), V(i,j+1)中.
当每个格点计算一遍时, 求总势能的变化$\Delta V_k$,如果 $$\Delta V_k>\varepsilon$$ , 用得到的$V_k$继续计算, 否则跳出循环, 且最终电势的分布 $$V=V_k$$ , 其中$\varepsilon$为给定的精度, 这次选择$\varepsilon=0.0001$.

SOR方法

因为Gauss-Seidel方法收敛到$\varepsilon$比较慢, 为了解决这个问题引入这种方法.
先按照Gauss-Seidel方法计算得到$V_{temp}(i,j)$, 根据该点旧的电势算出

$$\Delta V(i,j)=V_{temp}(i,j)-V_{old}(i,j)$$ 新的电势根据下式求出: $$V_{new}(i,j)=\alpha \Delta V(i,j) + V_{old}(i,j)$$
判断是否符合期望的精度的方法与Gauss-Seidel方法一致.

结果展示

算出课本上的电势图

• L=31, Gauss-Seidel方法
  

• L=31, SOR方法
  

改变计算次数

• L=31, SOR方法(不限制$\varepsilon$店的值)
  N = 30
  N = 18
  N = 9
  可以看出计算次数对计算结果的影响很大, 但是在计算次数足够多时, 即便没有满足精度值, 再增加计算次数对势分布的形状影响不大

改变计算格点数

• $\pm 1$的势能杆子仍然放在$x=\pm 3$位置, 且长也为$\pm 3$
  改变L的大小, 算出两种计算方法满足精度的计算次数.
  
  用二次和一次函数拟合. 其中拟合时把左上角的图中前4个点去掉, 右上角去掉倒数第3个点.
  $N = 0.2385 L^2 - 3.846 L + 15.5 \tag 1$
  $N = 2.081L - 5.274 \tag 2$
• 左边的图里在$L=22$和$L=24$明显得不符合二次曲线, 右图在L=46处计算次数也有骤降. 接下来探究原因.
  在$L=22$和$L=24$两边看起来满足二次方. 猜想是否与杆子的长短和位置没有与格点数等比变化有关, 这次设杆长和到原点的距离为$L/10$, 而$(L+1)^2$为格点数目. 用SOR方法计算得到

L	10	20	30	40	50	60	70	80	90	100
N	15	29	56	72	110	113	91	118	168	185

在60到70仍有一个下降, 其余位置基本以1.6为斜率增加. 不懂...
当然显然的结论是: SOR计算方法的计算次数与具体系统的性质有关.