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

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

@author: Mr.Yang
"""


""" 

函数形如
                Uxx+Uyy=f(x)   狄利克雷条件    步长一样的时候 h1 = h2 =h
                ▽•△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)   #作业
    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)             #作业
    return(f1)
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
def U_f(x,y):    #上边界
    #f=0                                          #方程 1
    f=0                                          #方程2
    #f=np.exp(x+1/2) #作业
    return(f)
def D_f(x,y):    #下边界
    #f=0                                          #方程 1
    f=0                                          # 方程2
    #f=np.exp(x+0)  #作业
    return(f)
def L_f(x,y):    #左边界
    #f=np.sin(np.pi*y)                            #方程 1
    f=0                                         #方程2
    #f=np.exp(0+y/2)              #作业
    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)    #作业
    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))
    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)
    for l in range(0,m-1):
    #for i in range(0,2):
            T[l][l] = T[l][l] -4
            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]=-4
    print(T)
    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):
            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]+I
            
            A2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=A2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+I
            
            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]+np.zeros((m,m))     
           
    #print(A2)
    A2[m*m-m:m*m,m*m-m:m*m]=T
    #A=A/(h*h)
    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])-D_f(x[1],y[0])
                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])-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]) #这个地方有问题 上下边值是否一样 不一样需要修改
            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])-U_f(x[1],y[m+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])-U_f(x[m],y[m+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]) #4月5号 由 D_f -.>U_F()
            
            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])
                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])
                else:
                    b[j]=F_f(x[j-(k*m-1)],y[k+1])*h*h
    
    #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([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)
    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,m+2))  # 定义完整的数值解
    #P_U1= np.zeros((m+2,m+2)) #完整的数值解
    P_U=e_fn(x1,y1).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.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)
    '''
    uex = np.zeros((m,m))
    err=np.zeros((m,m))  
    for i in range ( 0, m):
        for j in range(0,m):
            uex[i][j] = e_fn (x[i+1],y[j+1])
    err= abs (uex - U)   # 实际值 - 精确值 =误差
    print('------------------------------------------------------------')
    print('相同坐标下的精确值为：')
    print(uex)
    print('------------------------------------------------------------')
    print('误差如下：')
    print(err)
    '''
    #------------------------------------------------------------------------------
    #误差分析 结果以及  三维图形 程序如下
    err=np.zeros((m+2,m+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, m + 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('')