【电磁场与电磁波】有限差分法求解边值问题--PYTHON代码实现

sobreee

于 2024-05-07 15:09:07 发布

阅读量159

点赞数 1

文章标签： python 矩阵线性代数

本文链接：https://blog.csdn.net/weixin_53338444/article/details/125244963

版权

import numpy as np
import matplotlib.pyplot as plt

# 定义待求解9个节点电位矩阵 是一个5*5矩阵
phi_array = np.zeros((5, 5))
# print(phi_array)

# 定义缓存矩阵5*5
cache_array = np.zeros((5, 5))

# 定义误差限δ  可修改此处参数实现改变停止迭代的误差门限
Delta = 0.001
print("误差门限为：", Delta)


# 修改边界已知的16个电位值 基于本程序只需修改上边界电位为40V即可 可修改此处实现改变上盖板的电位
def phi_initial():
    for i in range(5):  # 设置上电位板的电位
        phi_array[i, 4] = 40
    for_list = [1, 2, 3]
    for i in for_list:
        phi_array[i, 1] = 10  # 设置第二列初值
        phi_array[i, 2] = 20  # 设置第三列初值
        phi_array[i, 3] = 30  # 设置第四列初值
    print("迭代初值为：", phi_array)
    return


# 定义迭代函数Iteration_fun:
def Iteration_fun(phi_im_j, phi_i_jm, phi_ip_j, phi_i_jp):
    phi_i_j_new = (1/4)*(phi_im_j+phi_i_jm+phi_ip_j+phi_i_jp)
    return phi_i_j_new


# 迭代主循环函数
def cycle_fun():
    for_list = [1, 2, 3]  # 定义迭代范围
    for i in for_list:
        for j in for_list:
            phi_array[i, j] = Iteration_fun(phi_array[i-1, j], phi_array[i, j-1], phi_array[i+1, j], phi_array[i, j+1])
        # 迭代函数执行↑
    return


# 获得迭代值缓冲列表， 目的是做出优化曲线
def suffer_list(x, y):
    new_value = Iteration_fun(phi_array[x-1, y], phi_array[x, y-1], phi_array[x+1, y], phi_array[x, y+1])
    return new_value


# 该函数功能是将更新值载入缓冲矩阵cache_array，用于判断
def suffer_save():
    for i in range(5):
        for j in range(5):
            cache_array[i, j] = phi_array[i, j]
    return


# 求解差值最大值函数，用于判断循环条件
def find_dif_max():
    dif_max = []   # 定义搜寻一维矩阵缓冲区
    dif_array = phi_array-cache_array  # 将迭代前后矩阵相减
    for_list = [1, 2, 3]  # 定义搜寻范围
    # 取绝对值
    for i in for_list:
        for j in for_list:
            dif_array[i, j] = abs(dif_array[i, j])
    # 取二维向量每一行最大值  变为一维  然后利用max函数取得每一行最大值其中的 最大值并作为函数的返回值
    for i in for_list:
        dif_max.append(max(dif_array[i]))
        # print('差值的最大值为：', max(dif_max))
    return max(dif_max)

# 定义迭代条件缓存值 即max小于误差限
# dif_max = max(1, 2)
# print(dif_max)


# 电位矩阵初始化
phi_initial()

# 定义迭代次数t
t = 0

# 定义迭代值缓冲区
y11 = [phi_array[1, 1]]
y12 = [phi_array[1, 2]]
y13 = [phi_array[1, 3]]
y21 = [phi_array[2, 1]]
y22 = [phi_array[2, 2]]
y23 = [phi_array[2, 3]]
y31 = [phi_array[3, 1]]
y32 = [phi_array[3, 2]]
y33 = [phi_array[3, 3]]
# cycle_fun()
# print(len(phi_array))
# 优化主循环
while 1:
    suffer_save()
    cycle_fun()
    t = t+1
    y11.append(suffer_list(1, 1))
    y12.append(suffer_list(1, 2))
    y13.append(suffer_list(1, 3))
    y21.append(suffer_list(2, 1))
    y22.append(suffer_list(2, 2))
    y23.append(suffer_list(2, 3))
    y31.append(suffer_list(3, 1))
    y32.append(suffer_list(3, 2))
    y33.append(suffer_list(3, 3))
    dif_max = find_dif_max()
    if dif_max < Delta:
        break
print(phi_array)
print("迭代次数为:", t)

'''
# plt测试
x = range(9)
y = range(9)
plt.plot(x, y)
plt.show()
'''
# 显示其中某点的优化曲线  可通过修改主迭代循环中的suffer_list函数的参数，实现展示不同点位的迭代曲线
x = range(t+1)
plt.plot(x, y11)
plt.plot(x, y12)
plt.plot(x, y13)
plt.plot(x, y21)
plt.plot(x, y22)
plt.plot(x, y23)
plt.plot(x, y31)
plt.plot(x, y32)
plt.plot(x, y33)
plt.legend(["phi(1,1)", "phi(1,2)", "phi(1,3)", "phi(2,1)",
            "phi(2,2)", "phi(2,3)", "phi(3,1)", "phi(3,2)", "phi(3,3)"], loc=0)
plt.xlabel("time")
plt.ylabel("phi")
plt.title("Potential changes with time")
plt.show()