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

# -*- coding: utf-8 -*-
"""
Created on Wed Apr  1 15:41:05 2020

@author: Mr.Yang
"""


""" 
                  -△U=-(Uxx+Uyy)
函数形如
                -(aUx)x+(aUy)y+(b1U)x+(b2U)y+cU=f(x,y)
                即为：
                -▽•(a(x,y)▽U)+(b1U)x+(b2U)y+cU=f(x,y)
                
                
                U|F=q(x,y)   边界条件
        五点差分格式化为如下形式：（程序设计按此顺序 改变系数函数顺序  将需要 重调程序  否则 会错误）
        
        A(i,j)U(i-1,j)+B(i,j)U(i+1,j)+C(i,j)U(i,j-1)+D(i,j)U(i,j+1)-E(i,j)U(i,j)=f(i,j)
        
        此方程的形式 常带有 混合边界条件 的形式
        
        
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy.linalg as la
import time
#---------------------------------------------------------------------
#定义函数   先定义后调用
def a_f(x,y):             #原方程 系数函数 a(x,y)
    a=1
    return(a)

def b1_f(x,y):            #原方程 系数函数 b1(x,y)
    b1=0
    return(b1)

def b2_f(x,y):           #原方程 系数函数 b2(x,y)
    b2=0
    return(b2)

def c_f(x,y):             #原方程 系数函数 c(x,y)
    c=0
    return(c)
#-----------------------------------------------------------------
def F_f(x,y):  #偏微分方程 的右端函数
    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=0                                              #方程3
    #f=2*np.pi*np.pi*np.sin(np.pi*x)*np.sin(np.pi*y)  #方程 4
    #f=-(y*y*np.exp(y)+x*x*np.exp(x))*np.exp(x*y)+(y*y+x*x+1)*np.exp(x*y)  #变系数方程1
    return(f)
    
def e_fn(x,y):  # 偏微分方程的解析解 函数 精确解函数
    f=np.exp(x)*np.sin(np.pi*y)                     #方程1
    #f=np.sin(np.pi*x)*np.sin(np.pi*y)                  #方程2   
    #f=np.exp(x)*(np.sin(y)+np.cos(y))               #方程3
    
    #f=np.exp(x*y)                                    #变系数方程1
    return(f) 


#-------------------------------------------------------------------------
                 #差分格式系数 都乘上了 h*h 左端的h*h乘到了  右端函数中去了
                 #五点差分格式的系数函数A(x,y)  顺序不能 乱调
                 # 此部分 需要改的就是 h的大小  其他为固定的
def A_f(x,y):    
    h=1/4
    f=a_f(x-(1/2)*h,y)-b1_f(x,y)*h/2
    return(f)

def B_f(x,y):       #五点差分格式的系数函数B(x,y)  顺序不能 乱调
    h=1/4
    f=a_f(x+(1/2)*h,y)+b1_f(x,y)*h/2
    return(f)

def C_f(x,y):   #五点差分格式的系数函数C(x,y)  顺序不能 乱调
    h=1/4
    f=a_f(x,y-(1/2)*h)-b2_f(x,y)*h/2
    return(f)

def D_ff(x,y):  #这个是系数f函数  #五点差分格式的系数函数D(x,y)  顺序不能 乱调
    h=1/4
    f=a_f(x,y+(1/2)*h)+b2_f(x,y)*h/2
    return(f)
def E_f(x,y):     #五点差分格式的系数函数E(x,y)  顺序不能 乱调
    f=A_f(x,y)+B_f(x,y)+C_f(x,y)+D_ff(x,y)+c_f(x,y)
    return(f)


#-------------------------------------------------------------------------
def U_f(x,y):    #上边界函数
    f=np.exp(x)*np.sin(np.pi)                                          #方程 1
    #f=0                                          #方程2
    #f=np.exp(x)*(np.sin(1)+np.cos(1))            #方程3
    #f=np.exp(x)                                  #变系数方程1
    return(f)
def D_f(x,y):    #下边界函数
    f=0                                          #方程 1
    #f=0                                          # 方程2
    #f=np.exp(x)                                  #方程3
    #f=1                                   #变系数方程1
    return(f)
def L_f(x,y):    #左边界函数
    f=np.sin(np.pi*y)                            #方程 1
    #f=0                                         #方程2
    #f=np.sin(y)+np.cos(y)                        #方程3
    #f=1                             #变系数方程1
    return(f)
def R_f(x,y):    #右边界函数
    f=np.exp(1)*np.sin(np.pi*y)                  #方程 1
    #f=0                                         #方程2
    #f=np.exp(1)*(np.sin(y)+np.cos(y))           #方程 3
    #f=np.exp(y)                                #变系数方程1
    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 项
    y = np.linspace ( m_a, m_b, m +2)
    #print(x)
    #系数矩阵A的赋值处理方法如下：
    T=np.zeros((m,m))
    I1=np.zeros((m,m))
    I2=np.zeros((m,m))
    I3=np.zeros((m,m))

    I=np.eye(m,k=0,dtype=int)  #https://blog.csdn.net/u013066730/article/details/56701436
    b=np.zeros((m*m))
    #print(b)
    #print('原始矩阵T如下\n')
    #print(T)
        
    print(T)
    
        #for i in range(0,)
    
    print('矩阵I如下\n')
    #print(I.shape)       
    print(I)
    #A1=np.zeros((m,m))
    #A2=np.tile(A1, (m, m))
    A2=np.zeros((m*m, m*m))
    #print('分块矩阵如下A2：')
    #print(A2.shape)
    #print(A2)
    
    # 类似分块矩阵       
    for j in range(0,m-1):
        #for k in range(1,m+1):
        for l in range(0,m-1):
    #for i in range(0,2):
            T[l][l] = T[l][l] -E_f(x[l+1],y[j+1])
            T[l][l+1]=T[l][l+1]+B_f(x[l+1],y[j+1])
            T[l+1][l]=T[l+1][l]+A_f(x[l+2],y[j+1])
            T[l+1][l+1]=T[l+1][l+1]+0
        
            I1[l][l]=I1[l][l]+D_ff(x[l+1],y[j+1])
        
            I2[l][l]=I2[l][l]+C_f(x[l+1],y[j+2])
        
                
        T[m-1][m-1]=-E_f(x[m],y[j+1])
        I1[m-1][m-1]=D_ff(x[m],y[j+1])
        I2[m-1][m-1]=C_f(x[m],y[j+2])
        
        #结束 子矩阵  的赋值   接下来是大矩阵的  赋值 
        A2[j*m:(j+1)*m,j*m:(j+1)*m]=A2[j*m:(j+1)*m,j*m:(j+1)*m]+T
            
        A2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]=A2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]+I1
            
        A2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=A2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+I2
            
        A2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]=A2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]+I3     
        
        # 将 T I1 I2 清零  不然 会累加
        T=np.zeros((m,m))
        I1=np.zeros((m,m))
        I2=np.zeros((m,m))
    
    #处理 几行 几列
    for l in range(0,m-1):
    #for i in range(0,2):
            T[l][l] = T[l][l] -E_f(x[l+1],y[m])
            T[l][l+1]=T[l][l+1]+B_f(x[l+1],y[m])
            T[l+1][l]=T[l+1][l]+A_f(x[l+2],y[m])
            T[l+1][l+1]=T[l+1][l+1]+0
        
    T[m-1][m-1]=-E_f(x[m],y[m])
    
    A2[m*m-m:m*m,m*m-m:m*m]=T
    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,m):
        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])*h*h-L_f(x[0],y[k+1])*A_f(x[j-(k*m-1)],y[k+1])-D_f(x[1],y[0])*C_f(x[j-(k*m-1)],y[k+1])
                elif(j==m*k+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-R_f(x[m+1],y[k+1])*B_f(x[j-(k*m-1)],y[k+1])-D_f(x[m],y[0])
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-D_f(x[j-m*k],y[0])*C_f(x[j-(k*m-1)],y[k+1]) #这个地方有问题 上下边值是否一样 不一样需要修改
            elif(k==m-1): #最后一行的b
                if(j==k*m):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-L_f(x[0],y[k+1])*A_f(x[j-(k*m-1)],y[k+1])-U_f(x[1],y[m+1])*D_ff(x[j-(k*m-1)],y[k+1])
                elif(j==m*k+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-R_f(x[m+1],y[k+1])*B_f(x[j-(k*m-1)],y[k+1])-U_f(x[m],y[m+1])*D_ff(x[j-(k*m-1)],y[k+1])
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-U_f(x[j-m*k],y[m+1])*D_ff(x[j-(k*m-1)],y[k+1]) 
            
            else:     #中间几行的b
                if(j==k*m):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-L_f(x[0],y[k+1])*A_f(x[j-(k*m-1)],y[k+1])
                elif(j==k*m+m-1):
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h-R_f(x[m+1],y[k+1])*B_f(x[j-(k*m-1)],y[k+1])
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*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([m,m,])
    U=U.T        # 需要转置 求解结果顺序为 U11 U21 U31...,U12 U22 U32...,U13 U23 U33..
    print('------------------------------------------------------------')
    print('即所计算的(%dx%d)个数值解：'%(m,m))
    print(U)
    
    #------------------------------------------------------------------------------------------  
    #精确解的3D图形程序如下
    fig=plt.figure()
    ax=Axes3D(fig)
    x1,y1=np.meshgrid(x,y)
    H=e_fn(x1,y1)
    ax.plot_surface(x1,y1,H.T,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')
    
    #--------------------------------------------------------------------------
    #近似解的3D图形程序如下
    print ( '' )
    print('------------------------------------------------------------')
    print('精确解结果如下：')
    print(e_fn(x1,y1))
    # 为了画出 数值解的 三维形状图  将程序中所求出的 数值mXm的结果替换  精确解中 中心位置mxm的结果 因为边值条件是一样的
    P_U= np.zeros((m+2,m+2))  # 定义完整的数值解
    #P_U1= np.zeros((m+2,m+2)) #完整的数值解
    P_U=H.T          #先把精确解 赋值 给他
    P_U[1:m+1,1:m+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,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')
    
    
    
    #------------------------------------------------------------------------------
    #误差分析 结果以及  三维图形 程序如下
    err=np.zeros((m+2,m+2))
    err=abs(e_fn(x1,y1)-P_U)/(m*m)
    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,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)
    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          数值解           解析解          Error' )
    print ( '' )
    
    for i in range ( 0, m + 2 ):
        print ( '  (%d,%d)  %14.6g  %14.6g  %14.6g' % ( 2,i, P_U[2][i], e_fn(x[2],y[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('')

只打印第三行的结果，步长为:h=0.0909091
  Node          数值解           解析解          Error

  (2,0)               0               0      0.00540641
  (2,1)        0.338971        0.337909      0.00253122
  (2,2)         0.65048        0.648442     2.03794e-05
  (2,3)        0.909292        0.906443      0.00199137
  (2,4)         1.09444         1.09101      0.00316699
  (2,5)         1.19092         1.18719      0.00339163
  (2,6)         1.19092         1.18719      0.00258102
  (2,7)         1.09444         1.09101     0.000728446
  (2,8)        0.909292        0.906443      0.00209526
  (2,9)         0.65048        0.648442      0.00574815
  (2,10)        0.338971        0.337909       0.0100293
  (2,11)     1.46884e-16     1.46884e-16       0.0146961
计算程序运行的时间

0.06643400000007205