2D热传导方程的FDM-齐次边界和右端项

自我o救赎

已于 2023-03-05 09:19:04 修改

阅读量89

点赞数

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

于 2023-03-05 09:14:35 首次发布

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

版权

Numerical Method For PDE 专栏收录该内容

12 篇文章 0 订阅

订阅专栏

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 28 08:07:25 2020

@author: Mr
"""
""" 

函数形如       狄利克雷条件
                Ut=a**2*Uxx  (0<x<l, t>0)  
                
                U(0,t)=0       t>0
                U(l,t)=0       t>0
                U(x,0)=φ(x)    0<x<l

例子：
                Ut=a**2*Uxx  (0<x<pi, t>0)  
                
                U(0,t)=0       t>0
                U(l,t)=0       t>0
                U(x,0)=100    0<x<pi

"""
import numpy as np
import matplotlib.pyplot as plt
from  matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
#from scipy.integrate import quad #,dblquad,nquad
from sympy import symbols,integrate
import seaborn as sns
#from CAL.PyCAL import *
#font.set_size(20)https://wenku.baidu.com/view/e645feafc9d376eeaeaad1f34693daef5ef71393.html   分离变量法求解

#%%
def Main_f(X,T,N,m,InitialCondition,Boundary_left,Boundary_right,D):
    
    if(D==1):
        
        n=N
        x1=np.linspace(0,1,50)
        y1=InitialCondition(x1)#map(InitialCondition,x)
        plt.figure(figsize=(10,6))
        plt.plot(x1,y1)
        plt.xlabel('$x$',fontsize=15)
        plt.ylabel('$f(x)$',fontsize=15)
        plt.title('1DHeat Equation Initial Condition')
        x=X
        t=T
        h=(max(x)-min(x))/n  #  步长 x轴方向
        k=(max(t)-min(t))/m # 时间步长
    
        kappa=1
        rho=kappa*k/h/h
        print('rho')
        '''
        显式格式不能任意取时间和空间的网格点数，即M与N不能随意取值。
        我们称显式格式为条件稳定。特别地，需       要满足所谓CFL条件（Courant–Friedrichs–Lewy）：
        '''
        if(rho<=0.5):
            print('显示格式此时满足条件稳定：')
        print('rh0=%f' % (rho))
        U=np.zeros((n+1,m+1)) #
        #print(U)
        
        for i in range(1,n):  #表示列数
            U[i][0]=InitialCondition(x[i])
        for j in range(0,m+1):   #表示 行数
            U[0][j]=Boundary_left(t[j])
            U[n][j]=Boundary_right(t[j])
        #print('初始化的U：')
        #print(U)
    
        for j in range(0,m):#  (1,  3)
            for i in range(1,n):  # (0,  5)
                #rho*U[j-1][k]+(1-2*rho)*U[j][k]+rho*U[j+1][k]
                #U[i+1][j]=rho*U[i][j-1]+(1-2*rho)*U[i][j]+rho*U[i][j+1]
                U[i][j+1]=rho*U[i+1][j]+(1-2*rho)*U[i][j]+rho*U[i-1][j]
        #print(U.shape)
        
        plt.figure(figsize=(10,6))
        plt.plot(x, U.T[0,:])
        plt.plot(x, U.T[int(0.10/ k),:])
        plt.plot(x, U.T[int(0.15/ k),:])
        plt.plot(x, U.T[int(0.20/ k),:])
        plt.xlabel('$x$', fontsize = 15)
        plt.ylabel(r'$U(dot, tau)$', fontsize = 15)
        plt.title(u'1D Heat equation', fontsize = 15)
        plt.legend([r'$t = 0.$', r'$t = 0.10$', r'$t = 0.15$', r'$t = 0.20$'], fontsize = 15)
        
        
        
        fig=plt.figure(figsize=(10,6))
        ax=Axes3D(fig)
        x1,y1=np.meshgrid(t,x)
        surf=ax.plot_surface(x1,y1,U,rstride=1,cstride=1,cmap=cm.coolwarm,alpha=0.9)   #绘面'rainbow'
        ax.set_xlabel('t')
        ax.set_ylabel('x')
        ax.set_zlabel('U')
        ax.set_xticks([0,0.05,0.1,0.15,0.2])
        ax.set_yticks([0,0.2,0.4,0.6,0.8,1])
        ax.set_zticks([0,0.2,0.4,0.6,0.8,1])
        ax.set_title('Hot equation$ u_t=A u_{xx} $')     
        ax.grid(False)
        
        cbar=plt.colorbar(surf,shrink=0.5,aspect=5)
        plt.show()
    #%%
        fig1=plt.figure(figsize=(10,6))
        #ax1=Axes3D(fig1)
        ax1=fig1.add_subplot(1,1,1,projection='3d')
        x1,y1=np.meshgrid(x,t)
        surf1=ax1.plot_surface(x1,y1,U.T,cmap=cm.coolwarm)   #绘面'rainbow' rstride=1,cstride=1,cmap=cm.coolwarm,alpha=0.9
        
        ax1.set_xlabel('x',fontdict={"size":16})
        ax1.set_ylabel('t',fontdict={"size":16})
        ax1.set_zlabel('U',fontdict={"size":16})
        ax1.set_xticks([0,0.2,0.4,0.6,0.8,1])
        ax1.set_yticks([0,0.05,0.1,0.15,0.2])
        ax1.set_xticks([0,0.2,0.4,0.6,0.8,1])
        ax1.set_title('Hot equation$ u_t=A u_{xx} ,$')     
        ax1.grid(False)
        cbar=fig1.colorbar(surf1,shrink=0.75,aspect=5)
        plt.show()
        
    
#%%    
    if(D==2):
        n = N[0] # 定义x方向节点数
        o = N[1] # 定义y方向节点数
        m = m # 定义时间步长
        nu = 0.05 # 扩散系数
        h = (max(X[0])-min(X[0])) / (n)
        h1 = (max(X[1])-min(X[1])) / (o)
        sigma = 0.25 # 库朗数
        k =T/(m)    #k=sigma * h * h1 / nu  该系数可以用现用替换 
                    # 非常不理解，nu是扩散系数，不是传播速度？还是类比扩散这个传播速度
        
        x = X[0]#np.linspace(0,2,n)
        y = X[1]#np.linspace(0,2,o)
        # 边界条件
        u = Boundary_left
        un = Boundary_right
        # 初始条件
        u[int(0.5/h):int(1/h+1),int(0.5/h1):int(1/h1+1)] = 2
        
        #m= int(input("Input:"))
        
        # 定义函数进行求解
        for n in range(m + 1):
        	un = u.copy()
        	u[1:-1,1:-1] = (un[1:-1,1:-1] +
        		             nu * k / h**2 * (un[2: ,1:-1] - 2 * un[1:-1,1:-1] + un[0:-2,1:-1]) +
        		             nu * k / h1**2 * (un[1:-1,2:] - 2 * un[1:-1,1:-1] + un[1:-1, 0:-2]) )
        	u[0, :] = 1
        	u[-1, :] = 1
        	u[:, 0] = 1
        	u[:, -1] = 1
        
        fig = plt.figure(figsize=(8,4))
        ax = Axes3D(fig)
        x, y = np.meshgrid(x,y)
        surf1=ax.plot_surface(x,y,u[:],rstride=1,cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False)  #cmap='viridis'
        ax.set_xlim(0,2)
        ax.set_ylim(0,2)
        ax.set_zlim(1,2.5)
        ax.set_xlabel('$x$',fontdict={"size":14})
        ax.set_ylabel('$y$',fontdict={"size":14})
        ax.set_zlabel('U',fontdict={"size":14})
        ax.set_xticks([0,0.4,0.8,1.2,1.6,2])
        ax.set_yticks([0,0.4,0.8,1.2,1.6,2])
        ax.set_zticks([1,1.5,2,2.5])
        ax.set_title('Hot equation$ u_t=A u_{xx}+A u_{yy} $')     
        ax.grid(False)
        cbar=fig.colorbar(surf1,shrink=0.75,aspect=5)
        plt.show()
        
#%%   
print('===================')
'''
X一维变量的取值范围为
T 时间变量的取值范围
初始条件的形状
    N X方向的网格数
    M 时间T方向的网格数
    T 最大时间期限
    X 最大空间范围
    U 用来存储差分网格点上值得矩阵

'''
#%%
D=2  #热传导方程的维数
if(D==1):
    N,m=25,1500
    T,x=np.linspace(0,0.2,m+1),np.linspace(0,1,N+1)  #m+1项数列 
    Initia=lambda x :4.0*(1-x)*x                     #初值条件
    Boundary_left=lambda t:0                         #左边值
    Boundary_right=lambda t:0                        #右边值
    Main_f(x,T,N,m,Initia,Boundary_left,Boundary_right,D=1)

#%%
elif(D==2):
    n=30   # X轴方向的网格数目
    o=30   # Y方向的网格数目
    m=1000   # Z时间方向的网格数目   一般 考虑 X 和Y轴网格数目固定 时  时间方向的网格数目变化 的情况
    InitialCondition=0
    T=0.2             #取一个 所研究温度的最大值 
    N=np.array([n,o])
    Boundary_left= np.ones((n,o))
    Boundary_right= np.ones((n,o))
    X=np.array([np.linspace(0,2,n),np.linspace(0,2,o)])
    
    
    
    Main_f(X,T,N,m,InitialCondition,Boundary_left,Boundary_right,D=2)
else:
    print('3D及以上还在努力中,请期待!!!!!!!!!!!!!!!')

在这里插入图片描述

# -*- coding: utf-8 -*-
"""
Created on Tue Apr 28 21:17:39 2020

@author: Mr.Yang
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from  matplotlib import cm
def Main_f(X,T,N,m,InitialCondition,Boundary_left,Boundary_right):
    
    n = N[0] # 定义x方向节点数
    o = N[1] # 定义y方向节点数
    m = m # 定义时间步长
    nu = 0.05 # 扩散系数
    h = (max(X[0])-min(X[0])) / (n)
    h1 = (max(X[1])-min(X[1])) / (o)
    sigma = 0.25 # 库朗数
    k =T/(m)    #k=sigma * h * h1 / nu  该系数可以用现用替换 
                # 非常不理解，nu是扩散系数，不是传播速度？还是类比扩散这个传播速度
    
    x = X[0]#np.linspace(0,2,n)
    y = X[1]#np.linspace(0,2,o)
    # 边界条件
    u = Boundary_left
    un = Boundary_right
    # 初始条件
    u[int(0.5/h):int(1/h+1),int(0.5/h1):int(1/h1+1)] = 2
    
    #m= int(input("Input:"))
    
    # 定义函数进行求解
    for n in range(m + 1):
    	un = u.copy()
    	u[1:-1,1:-1] = (un[1:-1,1:-1] +
    		             nu * k / h**2 * (un[2: ,1:-1] - 2 * un[1:-1,1:-1] + un[0:-2,1:-1]) +
    		             nu * k / h1**2 * (un[1:-1,2:] - 2 * un[1:-1,1:-1] + un[1:-1, 0:-2]) )
    	u[0, :] = 1
    	u[-1, :] = 1
    	u[:, 0] = 1
    	u[:, -1] = 1
    
    fig = plt.figure(figsize=(8,4))
    ax = Axes3D(fig)
    x, y = np.meshgrid(x,y)
    surf1=ax.plot_surface(x,y,u[:],rstride=1,cstride=1,cmap=cm.coolwarm,linewidth=0,antialiased=False)  #cmap='viridis'
    ax.set_xlim(0,2)
    ax.set_ylim(0,2)
    ax.set_zlim(1,2.5)
    ax.set_xlabel('$x$',fontdict={"size":14})
    ax.set_ylabel('$y$',fontdict={"size":14})
    ax.set_zlabel('U',fontdict={"size":14})
    
    ax.set_xticks([0,0.4,0.8,1.2,1.6,2])
    ax.set_yticks([0,0.4,0.8,1.2,1.6,2])
    ax.set_zticks([1,1.5,2,2.5])
    ax.set_title('Hot equation$ u_t=A u_{xx}+A u_{yy} $')     
    ax.grid(False)
    cbar=fig.colorbar(surf1,shrink=0.75,aspect=5)
    plt.show()
    return u

n=30   # X轴方向的网格数目
o=30   # Y方向的网格数目
m=1000   # Z时间方向的网格数目   一般 考虑 X 和Y轴网格数目固定 时  时间方向的网格数目变化 的情况
InitialCondition=0
Boundary_left=0
Boundary_right=0
T=0.2             #取一个 所研究温度的最大值 
N=np.array([n,o])

u_xy0=lambda x,y: 0
u_xyt=lambda x,y:x*np.sin(np.pi*y)-y*np.sin(np.pi*x)#x**2*np.sin(y)-y**2*np.sin(x)

Boundary_left= np.ones((n,o))
Boundary_right= np.ones((n,o))
X=np.array([np.linspace(0,2,n),np.linspace(0,2,o)])
Main_f(X,T,N,m,InitialCondition,Boundary_left,Boundary_right)
#print(N.)
#print(L,R)