1D-Possion方程的FDM-Dirichlet-Neumann边界

# -*- coding: utf-8 -*-
"""
Created on Tue Mar 31 18:34:41 2020

@author: Mr.Yang
"""
""" 

函数形如
    u''(x)=f(x)     
    u(a)=A
    u(b)=B
例题1：u''(x)=sin(x)
      u(0)=
      u(1)=
      
例题2：u''(x)=x*x*exp(x)                       一个纽曼边值情况
      u'(0)=                                   一个纽曼边值情况
      u(1)=
      f(x)=x*x*exp(x)
      
例题3：
The PDE is defined for 0 < x < 1:
#      -u'' = f
#    with right hand side
#      f(x) = -(exact(x)'') 
#    and boundary conditions
#      u(0) = exact(0),
#      u(1) = exact(1).
#
#    The exact solution is:解析解
#      exact(x) = x * ( 1 - x ) * exp ( x )
#    The boundary conditions are
#      u(0) = 0.0 = exact(0.0),
#      u(1) = 0.0 = exact(1.0).
#    The right hand side is:
#      f(x) = x * ( x + 3 ) * exp ( x ) = - ( exact''(x) )

"""
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg as la
import time
#---------------------------------------------------------------------
#定义函数   先定义后调用
def f_x(x):
    #f=-np.sin(x)  #  u''(x)=sin(x)  u'=-cos(x)  u=-sin(x)   #方程 1
    #f=-x * ( x + 3 ) * np.exp ( x )    #注意这里的正负号和有限元不一样    #方程 2
    #f=np.exp(x)*(np.cos(x)-2*np.sin(x)-x*np.cos(x)-x*np.sin(x)) #注意这里的正负号和有限元不一样  #方程 3
    #f=np.sin(np.pi*x)      #方程 4
    f=np.exp(x)
    return(f)
    
def e_fn(x):
    #f1=np.sin(x)                                 #方程 1
    #f1=x * ( 1 - x ) * np.exp ( x )               #方程 2
    #f1=x*np.cos(x)                               #方程 3
    #f1=-1/np.pi**2*np.sin(np.pi*x)               #方程 4
    f1=np.exp(x)
    return(f1)
    
def de_fn(x):   
    #f1=np.sin(x)                                 #对应方程 1
    #f1=x * ( 1 - x ) * np.exp ( x )              #对应方程 2
    #f1=x*np.cos(x)                               #对应方程 3
    #f1=-1/np.pi**2*np.sin(np.pi*x)               #对应方程 4
    f1=np.exp(x)
    return(f1)    

#-------------------------------------------------------------------------

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 项

    print(x)
    #系数矩阵A的赋值处理方法如下：
    A=np.zeros((m+2,m+2))
    b=np.zeros((m+2))
    #print(b)

    print('原始矩阵A如下\n')
    print(A)
    #----------------求解纽曼边值 的方法一A的赋值  但不是 试用中心差商
    
    for l in range(0,m+1):
    #for i in range(0,2):
            A[l][l]=A[l][l]+0
            A[l][l+1]=A[l][l+1]+1
            A[l+1][l]=A[l+1][l]+1
            A[l+1][l+1]=A[l+1][l+1]-2
    print('赋值后的矩阵A如下\n')        
    print(A)
    
    A[0][0]=-h
    A[0][1]=h
    A[m+1][m]=0
    A[m+1][m+1]=h*h
    '''
    #----------------求解纽曼边值 的方法三 A的赋值  用中心差商
    
    for l in range(0,m+1):
    #for i in range(0,2):
            A[l][l]=A[l][l]+0
            A[l][l+1]=A[l][l+1]+1
            A[l+1][l]=A[l+1][l]+1
            A[l+1][l+1]=A[l+1][l+1]-2
    
    print('赋值后的矩阵A如下\n')        
    print(A)
    
    A[0][0]=0
    A[0][1]=0
    A[m+1][m]=0
    A[m+1][m+1]=h*h
            
    print('赋值后的矩阵A如下\n')        
    print(A)
    '''
    
    
    #print(A)
    
    A=A/(h*h)
    print('最终所得的矩阵A如下:\n')        
    print(A)
    
    #方程右端系数矩阵b的赋值处理方法如下


    for i in range(0,m+2):
        if(i==0):
            b[i]=b[i]+de_fn(m_a)
        elif(i==m+1):
            b[i]=b[i]+e_fn(m_b)
        else:
            b[i]=b[i]+f_x(x[i])
        
    #print(b)       
    #b[0]=b[0]-e_fn(m_a)/(h*h)
    #b[m-1]=b[m-1]-e_fn(m_b)/(h*h)

    print('最终所得的矩阵b如下：\n')
    print(b)
    
    b=b.T
    u = la.solve ( A, b)
    #b1=b.T   转置与否，结果 都一样 

    print('数值解计算结果u如下：')
    print(u)
    #u1 = la.solve ( A, b1)
    #print(u1)
    
    npp = 51
    xp = np.linspace (m_a,m_b,npp)
    up = np.zeros ( npp )
    for i in range ( 0, npp):
        up[i] = e_fn ( xp[i] )

    uex = np.zeros ( m + 2 )
    for i in range ( 0, m + 2 ):
        uex[i] = e_fn(x[i] )
    #
    #  Compare the solution and the error at the nodes.
    #
    print ( '' )
    print ( '  Node          数值解           解析解          Error' )
    print ( '' )
    err=np.zeros(m+2)  #是在使用范数时才零为数组的   不然 使用常数也可以解决
    for i in range ( 0, m + 2 ):
        err[i] = abs ( uex[i] - u[i] )
        print ( '  %4d  %14.6g  %14.6g  %14.6g' % ( i, u[i], uex[i], err[i] ) )
    

    #plt.plot ( x, U, 'bo-',xp, up, 'r.')
    plt.plot( x, u, 'bo-',lw=1,label="Numerical ")
    plt.plot( xp, up, 'r.',lw=6,label="Real")
    plt.savefig ( 'fem1d.png' )
    plt.xlabel('$x$', fontsize=14,color='blue')
    plt.ylabel('$F(x)$', fontsize=14,color='blue')

    plt.legend(fontsize=14)
    #
    #  Terminate.
    # 
    plt.show ( )
    print ( '<--------------------------------------------------->' )
    print('使用误差证明该方法是二阶的即 O(h^2)')
    print('在无穷范数情况下：')
    norm1=max(abs(err))
    print('无穷范数下的误差为：%14.6g' % (norm1))
    print ( '<--------------------------------------------------->' )
    print('在1范数情况下：')
    norm2=sum(abs(err))*h
    print('1范数下的误差为：%14.6g' % (norm2))
    print ( '<--------------------------------------------------->' )
    print('在2范数情况下：')
    norm3=np.sqrt(sum(abs(err)*(abs(err)).T)*h)
    print('2范数下的误差为：%14.6g' % (norm3))
    
    print ( 'PDE——有限差分：中心差商求解：：1D 纽曼边值 问题' )
    print ( '  Normal end of execution.' )
    
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('')

Node          数值解           解析解          Error

     0         1.00497               1       0.0049727
     1         1.01487         1.00995      0.00492352
     2         1.02487            1.02      0.00487435
     3         1.03497         1.03015      0.00482517
     4         1.04517          1.0404      0.00477599
     5         1.05548         1.05075      0.00472681
     6         1.06588         1.06121      0.00467763
     7         1.07639         1.07177      0.00462844
     8         1.08701         1.08243      0.00457926
     9         1.09773          1.0932      0.00453008
    10         1.10856         1.10408      0.00448089
    11         1.11949         1.11506      0.00443171
    12         1.13054         1.12616      0.00438252
    13          1.1417         1.13736      0.00433333
    14         1.15296         1.14868      0.00428415
    15         1.16434         1.16011      0.00423496
    16         1.17584         1.17165      0.00418577
    17         1.18745         1.18331      0.00413658
    18         1.19917         1.19509      0.00408738
    19         1.21102         1.20698      0.00403819
    20         1.22298         1.21899        0.003989
    21         1.23506         1.23112       0.0039398
    22         1.24726         1.24337      0.00389061
    23         1.25958         1.25574      0.00384141
    24         1.27202         1.26823      0.00379221
    25         1.28459         1.28085      0.00374301
    26         1.29729          1.2936      0.00369381
    27         1.31011         1.30647      0.00364461
    28         1.32306         1.31947      0.00359541
    29         1.33614          1.3326      0.00354621
    30         1.34935         1.34586      0.00349701
    31         1.36269         1.35925       0.0034478
    32         1.37617         1.37277      0.00339859
    33         1.38978         1.38643      0.00334939
    34         1.40353         1.40023      0.00330018
    35         1.41741         1.41416      0.00325097
    36         1.43143         1.42823      0.00320176
    37         1.44559         1.44244      0.00315255
    38          1.4599         1.45679      0.00310333
    39         1.47434         1.47129      0.00305412
    40         1.48893         1.48593       0.0030049
    41         1.50367         1.50071      0.00295569
    42         1.51855         1.51565      0.00290647
    43         1.53358         1.53073      0.00285725
    44         1.54877         1.54596      0.00280803
    45          1.5641         1.56134      0.00275881
    46         1.57959         1.57688      0.00270959
    47         1.59523         1.59257      0.00266036
    48         1.61102         1.60841      0.00261114
    49         1.62698         1.62442      0.00256191
    50         1.64309         1.64058      0.00251268
    51         1.65937          1.6569      0.00246345
    52          1.6758         1.67339      0.00241422
    53         1.69241         1.69004      0.00236499
    54         1.70917         1.70686      0.00231576
    55         1.72611         1.72384      0.00226653
    56         1.74321         1.74099      0.00221729
    57         1.76048         1.75832      0.00216805
    58         1.77793         1.77581      0.00211881
    59         1.79555         1.79348      0.00206957
    60         1.81335         1.81133      0.00202033
    61         1.83132         1.82935      0.00197109
    62         1.84947         1.84755      0.00192184
    63         1.86781         1.86594       0.0018726
    64         1.88632          1.8845      0.00182335
    65         1.90503         1.90325       0.0017741
    66         1.92392         1.92219      0.00172485
    67         1.94299         1.94132       0.0016756
    68         1.96226         1.96063      0.00162635
    69         1.98172         1.98014      0.00157709
    70         2.00137         1.99984      0.00152784
    71         2.02122         2.01974      0.00147858
    72         2.04127         2.03984      0.00142932
    73         2.06152         2.06014      0.00138006
    74         2.08197         2.08064      0.00133079
    75         2.10262         2.10134      0.00128153
    76         2.12348         2.12225      0.00123226
    77         2.14455         2.14336      0.00118299
    78         2.16582         2.16469      0.00113372
    79         2.18731         2.18623      0.00108445
    80         2.20902         2.20798      0.00103518
    81         2.23094         2.22995     0.000985903
    82         2.25308         2.25214     0.000936626
    83         2.27544         2.27455     0.000887347
    84         2.29802         2.29718     0.000838067
    85         2.32083         2.32004     0.000788784
    86         2.34386         2.34312       0.0007395
    87         2.36713         2.36644     0.000690214
    88         2.39063         2.38999     0.000640926
    89         2.41436         2.41377     0.000591636
    90         2.43833         2.43778     0.000542344
    91         2.46253         2.46204      0.00049305
    92         2.48698         2.48654     0.000443754
    93         2.51167         2.51128     0.000394457
    94         2.53661         2.53627     0.000345157
    95          2.5618          2.5615     0.000295855
    96         2.58724         2.58699     0.000246551
    97         2.61293         2.61273     0.000197245
    98         2.63888         2.63873     0.000147937
    99         2.66508         2.66498     9.86268e-05
   100         2.69155          2.6915     4.93145e-05
   101         2.71828         2.71828               0


<--------------------------------------------------->
使用误差证明该方法是二阶的即 O(h^2)
在无穷范数情况下：
无穷范数下的误差为：     0.0049727
<--------------------------------------------------->
在1范数情况下：
1范数下的误差为：    0.00251212
<--------------------------------------------------->
在2范数情况下：
2范数下的误差为：    0.00289326
PDE——有限差分：中心差商求解：：1D 纽曼边值 问题
  Normal end of execution.
计算程序运行的时间

0.27608379999992394