2D-Possion方程的FDM-齐次边界-2

该代码示例使用Python实现了基于有限元方法(FEM)求解二维偏微分方程(PDE)的狄利克雷条件。通过定义不同边界条件和函数,构建和求解线性系统的系数矩阵,得出数值解,并进行误差分析。代码中包含了矩阵组装、求解和误差可视化的过程。
摘要由CSDN通过智能技术生成
# -*- coding: utf-8 -*-
"""
Created on Wed Apr  1 15:41:05 2020

@author: Mr.Yang
"""


""" 

函数形如
                Uxx+Uyy=f(x)   狄利克雷条件      步长不一样的时候 h1 != h2   任意的时候
                ▽•△U=f(x)     注意符号  左边是没有 - 的
                U(a)=A
                U(b)=B
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy.linalg as la
import time
#---------------------------------------------------------------------
#定义函数   先定义后调用
def F_f(x,y):
    #f=0
    #f=-(np.pi*np.pi-1)*np.exp(x)*np.sin(np.pi*y)     #方程 1
    #f=-2*np.pi*np.pi*np.sin(np.pi*x)*np.sin(np.pi*y)  #方程2
    #f=1.25*np.exp(x+y/2)   #作业
    f=-1  #测试有限元
    return(f)
    
def e_fn(x,y):
    #f1=np.exp(x)*np.sin(np.pi*y)                     #方程1
    #f1=np.sin(np.pi*x)*np.sin(np.pi*y)                  #方程2          
    #f1=np.exp(x+y/2)             #作业
    f=0  #测试有限元
    return(f)
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
def U_f(x,y):    #上边界
    #f=0                                          #方程 1
    #f=0                                          #方程2
    #f=np.exp(x+1/2) #作业
    f=0  #测试有限元
    return(f)
def D_f(x,y):    #下边界
    #f=0                                          #方程 1
    #f=0                                          # 方程2
    #f=np.exp(x+0)  #作业
    f=0  #测试有限元
    return(f)
def L_f(x,y):    #左边界
    #f=np.sin(np.pi*y)                            #方程 1
    #f=0                                         #方程2
    #f=np.exp(0+y/2)              #作业
    f=0  #测试有限元
    return(f)
def R_f(x,y):    #右边界
    #f=np.exp(1)*np.sin(np.pi*y)                  #方程 1
    #f=0                                         #方程2
    #f=np.exp(1+y/2)    #作业
    f=0  #测试有限元
    return(f)
#--------------------------------------------------------------------------------
def Fdm_cd_1d(a,b,n):
    m_a=a
    m_b=b
    m=n #区间【a,b】之间有10个点即有U1 U2 U3 U4 U5 ...U8 U9 U10
    h=1/(m+1)
    x = np.linspace ( m_a, m_b, m +2)  #生成等差数列      x=jh
                                   #第一项为 m_a  
                                   #最后一项为 m_b
                                   #共有 m+2 项
    #m1=100
    m1=m
    h1=1/(m1+1)                               
    y = np.linspace ( m_a, m_b, m1 +2)
    #print(x)
    #系数矩阵A的赋值处理方法如下:
    
    T=np.zeros((m,m))
    
    b=np.zeros((m*m1))
    #print(b)
    #print('原始矩阵T如下\n')
    #print(T)
    for l in range(0,m-1):
    #for i in range(0,2):
            T[l][l] = T[l][l] -2
            T[l][l+1]=T[l][l+1]+1
            T[l+1][l]=T[l+1][l]+1
            T[l+1][l+1]=T[l+1][l+1]+0
    print('赋值后的矩阵T如下\n')        
    T[m-1][m-1]=-2
    print(T)
    
    F1=np.zeros((m*m1, m*m1))
    for j in range(0,m1-1):
            F1[j*m:(j+1)*m,j*m:(j+1)*m]=F1[j*m:(j+1)*m,j*m:(j+1)*m]+T
            
            F1[j*m:(j+1)*m,(j+1)*m:(j+2)*m]=F1[j*m:(j+1)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))
            
            F1[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=F1[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+np.zeros((m,m))
            
            F1[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]=F1[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))     
    F1[m*m1-m:m*m1,m*m1-m:m*m1]=T       
    #print(A2)
    print(F1)   # 得到 左边 1/(h**2)的矩阵  还有右边
   
    
    F2=np.zeros((m*m1, m*m1))
    for j in range(0,m1-1):
            F2[j*m:(j+1)*m,j*m:(j+1)*m]=F2[j*m:(j+1)*m,j*m:(j+1)*m]+np.eye(m,m)*(-2)
            
            F2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]=F2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]+np.eye(m,m)
            
            F2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=F2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+np.eye(m,m)
            
            F2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]=F2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))     
    F2[m*m1-m:m*m1,m*m1-m:m*m1]=np.eye(m,m)*(-2)       
    #print(A2)
    print(F2) 
    print(F1+F2)
    
    A2=F1/(h**2)+F2/(h**2)
    
    print('最终所得的矩阵A如下:\n')        
    print(A2)

    #方程右端系数矩阵b的赋值处理方法如下
    #  make a mistake:can't multiply sequence by non int of type float 
    #in R_f(x[m+1],[m])  no y
    for k in range(0,m1):
        for j in range(k*m,(k+1)*m):
            if(k==0):   #第一层的 b
                if(j==k*m):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)-D_f(x[1],y[0])/(h1**2)
                elif(j==m*k+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)-D_f(x[m],y[0])/(h1**2)
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-D_f(x[j-m*k],y[0])/(h**2) #这个地方有问题 上下边值是否一样 不一样需要修改
            elif(k==m-1): #最后一层的b
                if(j==k*m):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)-U_f(x[1],y[m+1])/(h1**2)
                elif(j==m*k+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)-U_f(x[m],y[m+1])/(h1**2)
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-U_f(x[j-m*k],y[m+1])/(h1**2) #4月5号 由 D_f -.>U_F()
            
            else:     #中间几行的b
                if(j==k*m):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)
                elif(j==k*m+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])
    
    #print('最初的矩阵b如下')
    #print(b)
    #b=b*h*h

    print('最终所得的矩阵b如下:\n')
    print(b)
    
    b=b.T
    u = la.solve ( A2, b)
    
    #print('最初的矩阵b如下')
    #print(b)
    #b=b*h*h

    print('最终所得的矩阵b如下:\n')
    print(b)
    b=b.T
    u = la.solve ( A2, b)
    #b1=b.T   转置与否,结果 都一样 
    print('------------------------------------------------------------')
    print('数值解计算结果u顺序依次如下:')
    print(u)
    print('------------------------------------------------------------')
    print('数值结果排列如下:')
    print ( '' )
    print('U11 U12 U13 U14...')
    print('U21 U22 U23 U24...')
    print('U31 U32 U33 U34...')
    U=u.reshape([m1,m,])
    U=U.T        # 需要转置 求解结果顺序为 U11 U21 U31...,U12 U22 U32...,U13 U23 U33..
    print('------------------------------------------------------------')
    print('即所计算的(%dx%d)个数值解:'%(m,m1))
    print(U)
    '''
    #------------------------------------------------------------------------------------------  
    #精确解的3D图形程序如下
    fig=plt.figure()
    ax=Axes3D(fig)
    x1,y1=np.meshgrid(x,y)
    ax.plot_surface(x1,y1,e_fn(x1,y1),rstride=1,cstride=1,cmap='rainbow',alpha=0.9)   #绘面
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('Real')
    ax.set_title('Figure Real')
    ax.grid(False)
    
    #--------------------------------------------------------------------------
    #近似解的3D图形程序如下
    print ( '' )
    print('------------------------------------------------------------')
    print('精确解结果如下:')
    print(e_fn(x1,y1).T)
    # 为了画出 数值解的 三维形状图  将程序中所求出的 数值mXm的结果替换  精确解中 中心位置mxm的结果 因为边值条件是一样的
    P_U= np.zeros((m+2,m1+2))  # 定义完整的数值解
    #P_U1= np.zeros((m+2,m+2)) #完整的数值解
    P_U=e_fn(x1,y1).T          #先把精确解 赋值 给他
    '''
    #could not broadcast input array from shape (3,5) into shape (3,4)  
    '''
    print(P_U[1:m+1,1:m1+1])
    P_U[1:m+1,1:m1+1]=U   #采取 内部  用数值解 替换 精确解  得到(m+2)X(m+2)的所有点的数值结果
    print ( '' )
    print('------------------------------------------------------------')
    print('完整的数值解如下:')
    print(P_U)
    #ax.view_init(elev=300,azim=300) #旋转曲面  参数1: 从哪个高度看曲线 参数2:旋转角度
    
    
    fig1=plt.figure()
    ax1=Axes3D(fig1)
    #x2,y2=np.meshgrid(x,y)
    Z1=np.array(list(P_U),dtype='float64')
    ax1.plot_surface(x1,y1,Z1.T,rstride=1,cstride=1,cmap='rainbow',alpha=0.9)   #绘面
    ax1.set_xlabel('x')
    ax1.set_ylabel('y')
    ax1.set_zlabel('Numerical')
    ax1.set_title('Figure Numerical')
    ax1.grid(False)
    #------------------------------------------------------------------------------
    #误差分析 结果以及  三维图形 程序如下
    err=np.zeros((m+2,m1+2))
    err=abs(e_fn(x1,y1).T-P_U)
    print ( '' )
    print('------------------------------------------------------------')
    print('误差结果如下:')
    print(err)
    fig1=plt.figure()
    ax2=Axes3D(fig1)
    #x2,y2=np.meshgrid(x,y)
    Z2=np.array(list(err),dtype='float64')
    surf=ax2.plot_surface(x1,y1,Z2.T,rstride=1,cstride=1,cmap='rainbow',alpha=0.9)   #绘面
    ax2.set_xlabel('x')
    ax2.set_ylabel('y')
    ax2.set_zlabel('Error')
    ax2.set_title('Figure error')
    #ax2.set_zlim(0,0.1)
    ax2.grid(False)
    cbar=plt.colorbar(surf,shrink=0.5,aspect=5)
    #cbar.set_ticks([.01,.02,.03,.04,.05,.06,.07,.08,.09,.1])
    #------------------------------------------------------------------------------------------
    #  Compare the solution and the error at the some nodes.
    #
    print('------------------------------------------------------------')
    print ( '' )
    print('只打印第三行的结果,步长为:h=%g' %(1/(m+1)))
    print ( '  Node          Num           True          Error' )
    print ( '' )
    t=e_fn(x1,y1).T
    for i in range ( 0, m1 + 2 ):
        print ( '  (%d,%d)  %14.6g  %14.6g  %14.6g' % ( 2,i, P_U[2][i], t[2][i], err[2][i] ) )

    '''
a,b,n= map(int,input('请输入左端点a,右端点b,单元格个数m:等3个系数:各数字之间空格键隔开\n').split())
start = time.clock()  #输入系数 Enter 后开始计时

Fdm_cd_1d(a,b,n)

end = time.clock()
print('计算程序运行的时间\n')
print(end-start)
print('')



  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值