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

最新推荐文章于 2024-07-29 14:37:54 发布

自我o救赎

最新推荐文章于 2024-07-29 14:37:54 发布

阅读量120

点赞数

分类专栏： Numerical Method For PDE 文章标签： python

本文链接：https://blog.csdn.net/y15520833229/article/details/129342649

版权

Numerical Method For PDE 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

# -*- coding: utf-8 -*-
"""
Created on Tue Mar 31 22:11:52 2020

@author: Mr.Yang
"""


""" 

函数形如
                a(x)u''(x)+b(x)u'(x)+c(x)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
#---------------------------------------------------------------------
#定义函数   先定义后调用
es=0.01  #奇异摄动二阶非齐次
def f_x(x):
    #f=-np.sin(x)  #  u''(x)=sin(x)  u'=-cos(x)  u=-sin(x) 
    #f=-x * ( x + 3 ) * np.exp ( x )    #注意这里的正负号和有限元不一样
    #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)
    f=-1/es       #  作业 
    return(f)
    
def e_fn(x):
    #f1=np.sin(x)
    #f1=x * ( 1 - x ) * np.exp ( x )
    #f=x*np.cos(x)  #方程3
    #f1=-1/np.pi**2*np.sin(np.pi*x)
    f=(np.exp(1/es)-2)/(np.exp(1/es)-1)+(np.exp(1/es*x))/(np.exp(1/es)-1)+x  #  作业 https://wenku.baidu.com/view/cde993ba4028915f814dc24e.html  解析解 求解网址
    return(f)
#-------------------------------------------------------------------------
def a_f(x):
    
    #f=1  #方程3
    f=1       #  作业 
    return(f)
def b_f(x):
    #f=0  #方程3
    f=-1/es     #  作业 
    return(f)
def c_f(x):
    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 项

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

    print('原始矩阵A如下\n')
    print(A)
    for l in range(0,m-1):
    #for i in range(0,2):
            A[l][l] = A[l][l] + h*h*c_f(x[l+1])-2*a_f(x[l+1])
            A[l][l+1]=A[l][l+1]+a_f(x[l+1])+h*b_f(x[l+1])/2
            A[l+1][l]=A[l+1][l]+a_f(x[l+2])-h*b_f(x[l+2])/2
            A[l+1][l+1]=A[l+1][l+1]+0
    print('赋值后的矩阵A如下\n')        
    print(A)

    A[m-1][m-1]=h*h*c_f(x[m])-2*a_f(x[m])
    A=A/(h*h)
    print('最终所得的矩阵A如下:\n')        
    print(A)

#方程右端系数矩阵b的赋值处理方法如下


    for i in range(0,m):
    #if(i==0):
    #    b[i]=b[i]+f_x(x[i+1])-(e_fn(m_a)/(h*h))
    #elif(i==m-1):
    #    b[i]=b[i]+f_x(x[i+1])-(e_fn(m_b)/(h*h))
    #else:
          b[i]=b[i]+f_x(x[i+1])
        
    print(b)       
    b[0]=b[0]-(a_f(x[1])/(h*h)-b_f(x[1])/(2*h))*e_fn(m_a)   #     边值的处理
    b[m-1]=b[m-1]-(a_f(x[m])/(h*h)+b_f(x[m])/(2*h))*e_fn(m_b)

    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 = n+2#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] )

    U=np.zeros((m+2))
    for j in range(0,m+2):
        if(j==0):
            U[j]=e_fn(m_a)
        elif(j==m+1):
            U[j]=e_fn(m_b)
        else:
            U[j]=u[j-1]
        
    print('')        
    print('已知有m个U的值，再加上左端点值，右端点值')
    print(U)

    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           Numer           True          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.fill_between(x,U ,up, color = 'green')
    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('Under inf Norm：')
    #print(abs(err))
    norm1=max(abs(err))
    print('The error of under inf Norm：%14.6g' % (norm1))
    print ( '<--------------------------------------------------->' )
    print('Under 1 Norm：')
    norm2=sum(abs(err))*h
    print('The error of under 1 Norm：%14.6g' % (norm2))
    print ( '<--------------------------------------------------->' )
    print('Under 2 Norm：')
    norm3=np.sqrt(sum(abs(err)*(abs(err)).T)*h)
    print('The error of under 2 Norm：%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           Numer           True          Error

     0               1               1               0
     1          1.0099          1.0099               0
     2          1.0198          1.0198     2.22045e-16
     3          1.0297          1.0297     2.22045e-16
     4          1.0396          1.0396     4.44089e-16
     5          1.0495          1.0495     4.44089e-16
     6         1.05941         1.05941     6.66134e-16
     7         1.06931         1.06931     8.88178e-16
     8         1.07921         1.07921     1.11022e-15
     9         1.08911         1.08911     1.77636e-15
    10         1.09901         1.09901      1.9984e-15
    11         1.10891         1.10891     2.44249e-15
    12         1.11881         1.11881     2.88658e-15
    13         1.12871         1.12871     2.88658e-15
    14         1.13861         1.13861     3.10862e-15
    15         1.14851         1.14851     3.10862e-15
    16         1.15842         1.15842     3.33067e-15
    17         1.16832         1.16832     3.77476e-15
    18         1.17822         1.17822      3.9968e-15
    19         1.18812         1.18812      3.9968e-15
    20         1.19802         1.19802     4.44089e-15
    21         1.20792         1.20792     4.66294e-15
    22         1.21782         1.21782     4.88498e-15
    23         1.22772         1.22772     5.32907e-15
    24         1.23762         1.23762     5.32907e-15
    25         1.24752         1.24752     5.55112e-15
    26         1.25743         1.25743     5.77316e-15
    27         1.26733         1.26733      5.9952e-15
    28         1.27723         1.27723     6.43929e-15
    29         1.28713         1.28713     6.88338e-15
    30         1.29703         1.29703     7.32747e-15
    31         1.30693         1.30693     7.77156e-15
    32         1.31683         1.31683     7.99361e-15
    33         1.32673         1.32673     8.21565e-15
    34         1.33663         1.33663     8.43769e-15
    35         1.34653         1.34653     8.88178e-15
    36         1.35644         1.35644     8.88178e-15
    37         1.36634         1.36634     9.32587e-15
    38         1.37624         1.37624     9.99201e-15
    39         1.38614         1.38614     1.06581e-14
    40         1.39604         1.39604     1.11022e-14
    41         1.40594         1.40594     1.15463e-14
    42         1.41584         1.41584     1.22125e-14
    43         1.42574         1.42574     1.26565e-14
    44         1.43564         1.43564     1.33227e-14
    45         1.44554         1.44554     1.37668e-14
    46         1.45545         1.45545     1.44329e-14
    47         1.46535         1.46535      1.4877e-14
    48         1.47525         1.47525     1.53211e-14
    49         1.48515         1.48515     1.62093e-14
    50         1.49505         1.49505     1.66533e-14
    51         1.50495         1.50495     1.70974e-14
    52         1.51485         1.51485     1.77636e-14
    53         1.52475         1.52475     1.84297e-14
    54         1.53465         1.53465     1.88738e-14
    55         1.54455         1.54455     1.93179e-14
    56         1.55446         1.55446      1.9762e-14
    57         1.56436         1.56436     2.04281e-14
    58         1.57426         1.57426     2.10942e-14
    59         1.58416         1.58416     2.15383e-14
    60         1.59406         1.59406     2.19824e-14
    61         1.60396         1.60396     2.26485e-14
    62         1.61386         1.61386     2.30926e-14
    63         1.62376         1.62376     2.37588e-14
    64         1.63366         1.63366     2.42029e-14
    65         1.64356         1.64356      2.4869e-14
    66         1.65347         1.65347     2.64233e-14
    67         1.66337         1.66337     2.86438e-14
    68         1.67327         1.67327     3.30846e-14
    69         1.68317         1.68317     4.44089e-14
    70         1.69307         1.69307     7.28306e-14
    71         1.70297         1.70297      1.4766e-13
    72         1.71287         1.71287       3.475e-13
    73         1.72277         1.72277     8.79741e-13
    74         1.73267         1.73267     2.29972e-12
    75         1.74257         1.74257     6.08713e-12
    76         1.75248         1.75248     1.61811e-11
    77         1.76238         1.76238     4.30553e-11
    78         1.77228         1.77228      1.1452e-10
    79         1.78218         1.78218     3.04303e-10
    80         1.79208         1.79208     8.07511e-10
    81         1.80198         1.80198     2.13942e-09
    82         1.81188         1.81188     5.65771e-09
    83         1.82178         1.82178     1.49302e-08
    84         1.83168         1.83168     3.93037e-08
    85         1.84158         1.84158     1.03178e-07
    86         1.85149         1.85149     2.69982e-07
    87         1.86139         1.86139     7.03796e-07
    88         1.87129         1.87129     1.82658e-06
    89         1.88119          1.8812     4.71583e-06
    90          1.8911         1.89111     1.20993e-05
    91         1.90101         1.90104     3.08084e-05
    92         1.91095         1.91103     7.77191e-05
    93         1.92096         1.92116      0.00019378
    94         1.93119         1.93167      0.00047596
    95         1.94208         1.94322      0.00114604
    96         1.95489         1.95757      0.00268483
    97         1.97341         1.97945       0.0060427
    98         2.00883         2.02159       0.0127598
    99         2.09427         2.11824        0.023968
   100         2.32785         2.36164       0.0337916
   101               3               3               0

<--------------------------------------------------->
使用误差证明该方法是二阶的即 O(h^2)
Under inf Norm：
The error of under inf Norm：     0.0337916
<--------------------------------------------------->
Under 1 Norm：
The error of under 1 Norm：    0.00080387
<--------------------------------------------------->
Under 2 Norm：
The error of under 2 Norm：    0.00436509
PDE——有限差分：中心差商求解：：1D 狄利克雷边值 问题
  Normal end of execution.
计算程序运行的时间

0.33772579999958907