# -*- 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 分离变量法求解
'''a=1
l=np.pi
def g(x1):
f=100
return f
def Fenli_f(x,t):
s=0
for i in range(1,100):#np.inf
x1=symbols('x1')
f=x1#g(x1)*np.sin((i*x1*np.pi)/l)
A=integrate(f,(x1,0,l))
s=s+A*(np.sin(i*np.pi*x/l)*np.exp(-((i*a*np.pi)/l)**2*t))
print(s)
#s=np.sin(x*t)
return s
x=np.linspace(0,np.pi,10)
t=np.linspace(0,1,5)
x1,t1=np.meshgrid(x,t)
z=Fenli_f(x1,t1)
fig1=plt.figure()
ax1=Axes3D(fig1)
ax1.plot_surface(x1,t1,z,rstride=1,cstride=1,cmap='rainbow',alpha=0.9) #绘面
ax1.set_xlabel('x')
ax1.set_ylabel('t')
ax1.set_zlabel('U')
ax1.set_title('Frnli variability method')
ax1.set_xlim(0,np.pi)
ax1.set_ylim(0,1)
ax1.set_zlim(0,100)
'''
#%%
def E_f(x,t):
B=150
k=0.02
f=np.exp(-(x-0.4)**2/(4*k*t+1/B))/(np.sqrt(4*B*k*t+1))
return f
def Main_f(x,T,N,m,InitialCondition,Boundary_left,Boundary_right):
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')
print(x)
x=x
t=T
h=(max(x)-min(x))/n # 步长 x轴方向
k=(max(t)-min(t))/m # 时间步长
kappa=0.02 #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)
U[i][j+1]=rho*U[i+1][j]+(1-2*rho)*U[i][j]+rho*U[i-1][j]
#%%
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_yticks([0,0.2,0.4,0.6,0.8,1])
#ax1.set_zticks([0,0.2,0.4,0.6,0.8,1])
ax1.set_title('Numeri Hot equation$ u_t=A u_{xx} ,$')
ax1.grid(False)
cbar=fig1.colorbar(surf1,shrink=0.75,aspect=5)
plt.show()
#%%
fig2=plt.figure(figsize=(10,6))
#ax1=Axes3D(fig1)
ax2=fig2.add_subplot(1,1,1,projection='3d')
x2,y2=np.meshgrid(x,t)
surf2=ax2.plot_surface(x2,y2,E_f(x2,y2),cmap=cm.coolwarm) #绘面'rainbow' rstride=1,cstride=1,cmap=cm.coolwarm,alpha=0.9
ax2.set_xlabel('x',fontdict={"size":16})
ax2.set_ylabel('t',fontdict={"size":16})
ax2.set_zlabel('U',fontdict={"size":16})
ax2.set_xticks([0,0.2,0.4,0.6,0.8,1])
#ax2.set_yticks([0,0.05,0.1,0.15,0.2])
ax2.set_yticks([0,0.2,0.4,0.6,0.8,1])
#ax2.set_zticks([0,0.2,0.4,0.6,0.8,1])
ax2.set_title('Real Hot equation$ u_t=A u_{xx} ,$')
ax2.grid(False)
cbar=fig2.colorbar(surf2,shrink=0.75,aspect=5)
plt.show()
#%%
fig3=plt.figure(figsize=(10,6))
#ax1=Axes3D(fig1)
ax3=fig3.add_subplot(1,1,1,projection='3d')
x3,y3=np.meshgrid(x,t)
surf3=ax3.plot_surface(x3,y3,abs(E_f(x2,y2)-U.T),cmap=cm.coolwarm) #绘面'rainbow' rstride=1,cstride=1,cmap=cm.coolwarm,alpha=0.9
ax3.set_xlabel('x',fontdict={"size":16})
ax3.set_ylabel('t',fontdict={"size":16})
ax3.set_zlabel('U',fontdict={"size":16})
ax3.set_xticks([0,0.2,0.4,0.6,0.8,1])
#ax3.set_yticks([0,0.05,0.1,0.15,0.2])
ax3.set_yticks([0,0.2,0.4,0.6,0.8,1])
#ax3.set_zticks([0,0.2,0.4,0.6,0.8,1])
ax3.set_title('Error Hot equation$ u_t=A u_{xx} ,$')
ax3.grid(False)
cbar=fig3.colorbar(surf3,shrink=0.75,aspect=5)
plt.show()
return U
print('===================')
'''
X一维变量的取值范围为
T 时间变量的取值范围
初始条件的形状
N X方向的网格数
M 时间T方向的网格数
T 最大时间期限
X 最大空间范围
U 用来存储差分网格点上值得矩阵
'''
# 热传导方程 案例1
N,m=25,1500 #6000,1500#
T,x=np.linspace(0,1,m+1),np.linspace(0,1,N+1) #m+1项数列 T 0.2
Initia=lambda x:E_f(x,0 )#lambda x :4.0*(1-x)*x #初值条件
Boundary_left=lambda t:E_f(0,t) #lambda t:0 #左边值
Boundary_right=lambda t:E_f(1,t) #lambda t:0 #右边值
print(Main_f(x,T,N,m,Initia,Boundary_left,Boundary_right))
# 热传导方程 案例2
输出结果
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
