Bohrium
robot
新建

空间站广场

论文
Notebooks
比赛
课程
Apps
我的主页
我的Notebooks
我的论文库
我的足迹

我的工作空间

任务
节点
文件
数据集
镜像
项目
数据库
公开
大作业 曹一帆 202105755309 第一题
python
python
bohre001cd
更新于 2024-07-14
推荐镜像 :FEALPy :2.0.1.1
推荐机型 :c2_m4_cpu
赞
1
定义模型数据
Du_FortFrankel
求解

定义模型数据

代码
文本
[23]
import numpy as np
from sympy import *
import sympy as sp
#利用sympy求出右端项f(x,y,t)
x=symbols('x')
y=symbols('y')
t=symbols('t')

u=exp(-t)*sin(pi*x)*sin(pi*y)
f= diff(u,t,1)-(diff(u,x,2)+diff(u,y,2))
u_0= u.subs(t, 0)

print("u=", u)
print("f=", f)
print("u(x,y,0)=", u_0)

u=sp.Array(u)
u = sp.lambdify((x,y,t), u, 'numpy')
f=sp.Array(f)
f= sp.lambdify((x,y,t),f , 'numpy')
u_0=sp.Array(u_0)
u_0 = sp.lambdify((x,y), u_0, 'numpy')

class SinExpPDEData:
def __init__(self, D=[0, 1, 0, 1], T=[0, 1]):
"""
@brief 模型初始化函数
@param[in] D 模型空间定义域
@param[in] T 模型时间定义域
"""
self._domain = D
self._duration = T

def domain(self):
"""
@brief 空间区间
"""
return self._domain

def duration(self):
"""
@brief 时间区间
"""
return self._duration

def solution(self, p, t):
"""
@brief 真解函数

@param[in] p numpy.ndarray, 空间点
@param[in] t float, 时间点

@return 真解函数值
"""
pi = np.pi
x = p[..., 0]
y = p[..., 1]
return u(x,y,t)

def init_solution(self, p):
"""
@brief 初始解

@param[in] p numpy.ndarray, 空间点
@param[in] t float, 时间点

@return 初始解函数值
"""
pi = np.pi
x = p[..., 0]
y = p[..., 1]
return u_0(x,y)
def source(self, p, t):
"""
@brief 方程右端项

@param[in] p numpy.ndarray, 空间点
@param[in] t float, 时间点

@return 方程右端函数值
"""
pi = np.pi
x = p[..., 0]
y = p[..., 1]
return f(x,y,t)
def gradient(self, p, t):
"""
@brief 真解导数

@param[in] p numpy.ndarray, 空间点
@param[in] t float, 时间点

@return 真解导函数值
"""
x = p[..., 0]
y = p[..., 1]
pi = np.pi
val = np.zeros(p.shape, dtype=np.float64)
u_x=diff(u, x, 1)
u_y=diff(u, y, 1)
u_t=diff(u, t, 1)
u_x=sp.Array(u_x)
u_x = sp.lambdify((x,y, t), u_x, 'numpy')
u_y=sp.Array(u_y)
u_y = sp.lambdify((x,y, t), u_y, 'numpy')
u_t=sp.Array(u_t)
u_t = sp.lambdify((x,y, t), u_t, 'numpy')
val[..., 0] = u_x(x,y,t)
val[..., 1] = u_y(x,y,t)
val[..., 2] = u_t(x,y,t)
return val
def dirichlet(self, p, t):
"""
@brief Dirichlet 边界条件

@param[in] p numpy.ndarray, 空间点
@param[in] t float, 时间点
"""
return self.solution(p, t)


pde = SinExpPDEData()
domain = pde.domain()
# 测试
print('\nquestion:')
print('domain :', domain)
print(np.abs(pde.solution(np.array([0, 0]),0) < 1e-12))
print(np.abs(pde.solution(np.array([0.5, 0.5]),1)-exp(-1))< 1e-12)
print(np.abs(pde.source(np.array([0.5, 0.5]),0)-(-1+2*pi*pi))< 1e-12)
print(np.abs(pde.init_solution(np.array([0.5, 0.5]))-1)< 1e-12)

u= exp(-t)*sin(pi*x)*sin(pi*y)
f= -exp(-t)*sin(pi*x)*sin(pi*y) + 2*pi**2*exp(-t)*sin(pi*x)*sin(pi*y)
u(x,y,0)= sin(pi*x)*sin(pi*y)

question:
domain : [0, 1, 0, 1]
True
True
True
True
代码
文本

Du_FortFrankel

代码
文本
[26]
from fealpy.mesh.uniform_mesh_2d import UniformMesh2d
import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse.linalg import spsolve
from scipy.sparse import diags
from sympy import *
import sympy as sp
from scipy.sparse import coo_matrix, csr_matrix, diags, spdiags

def Du_FortFrankel(self, tau):
rx = tau/self.h[0]**2
ry = tau/self.h[1]**2

NN = self.number_of_nodes()
n0 = self.nx + 1
n1 = self.ny + 1
K = np.arange(NN).reshape(n0, n1)

val1 = np.broadcast_to(-rx, (NN-n1, ))
val2 = np.broadcast_to(-ry, (NN-n0, ))
I = K[1:, :].flat
J = K[0:-1, :].flat
A = diags([1 + rx + ry], [0], shape=(NN, NN), format='csr')
A += csr_matrix((val1, (I, J)), shape=(NN, NN), dtype=self.ftype)
A += csr_matrix((val2, (J, I)), shape=(NN, NN), dtype=self.ftype)

B = diags([1], [0], shape=(NN, NN), format='csr')
B += csr_matrix((-val1, (J, I)), shape=(NN, NN), dtype=self.ftype)
B += csr_matrix((-val2, (I, J)), shape=(NN, NN), dtype=self.ftype)
C = diags([-rx-ry], [0], shape=(NN, NN), format='csr')
return A, B, C

def DFF(n):
global uh0
global uh1
"""
@brief 时间步进格式为 CN 方法
@param[in] n int, 表示第 n 个时间步（当前时间步）
"""
t = duration[0] + n*tau
if n == 0:
return uh0, t
elif n == 1:
A1, B1 = mesh.parabolic_operator_crank_nicholson(tau)
source = lambda p: pde.source(p, t )
f = mesh.interpolate(source, intertype='node') # f.shape = (nx+1,ny+1)
f *= tau
f.flat[:] += B1@uh0.flat[:]
gD = lambda p: pde.dirichlet(p, t )
A1, f = mesh.apply_dirichlet_bc(gD, A1, f)
uh1.flat = spsolve(A1, f)
return uh1,t
else:
A, B, C = Du_FortFrankel(mesh,tau)
source1=lambda p: pde.source(p, t+tau)
source2=lambda p: pde.source(p, t)
source3=lambda p: pde.source(p, t-tau)
f1 = mesh.interpolate(source1, intertype='node')
f2 = mesh.interpolate(source2, intertype='node')
f3 = mesh.interpolate(source3, intertype='node')
# print("f1",f1)
f=f1+2*f2+f3
# f=ff(f1)
f *= tau/4
# print("f",f)
f.flat[:]+= B@uh1.flat[:]+C@uh0.flat[:]
gD = lambda p: pde.dirichlet(p, t+tau)
A, f = mesh.apply_dirichlet_bc(gD, A, f)

uh0=uh1;
uh1.flat = spsolve(A, f)
solution = lambda p: pde.solution(p, t+tau)
e = mesh.error(solution, uh1, errortype='max')
print(f"the max error is {e}")

return uh1, t
代码
文本

求解

代码
文本
[29]
from mpl_toolkits.mplot3d import Axes3D
nx = 20
ny = 20
hx = (domain[1] - domain[0])/nx
hy = (domain[3] - domain[2])/ny
mesh = UniformMesh2d([0, nx, 0, ny], h=(hx, hy), origin=(domain[0], domain[2]))

# 时间离散
duration = pde.duration()
nt =800
tau = (duration[1] - duration[0])/nt
global uh1
global uh0


from mpl_toolkits.mplot3d import Axes3D
maxit = 3

print(1)
for i in range(maxit):
uh1 = mesh.interpolate(pde.init_solution, intertype='node')
uh0 = np.zeros(mesh.number_of_nodes()).reshape(mesh.nx+1,mesh.ny+1)
for n in range(nt):
uh1,_=DFF(n)
x=np.arange(0,1+mesh.h[0],mesh.h[0])
y=np.arange(0,1+mesh.h[1],mesh.h[1])
x,y=np.meshgrid(x,y)
z=pde.solution(mesh.node,1)#真解
fig=plt.figure()
ax=fig.add_subplot(1,1,1,projection='3d')
surf=ax.plot_surface(x,y,z, cmap='rainbow')#真解图像
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5) # 添加颜色标识
plt.show()
fig=plt.figure()
ax1=fig.add_subplot(1,1,1,projection='3d')
surf1=ax1.plot_surface(x,y,uh1, cmap='rainbow')#拟合图像
fig.colorbar(surf1, ax=ax1, shrink=0.5, aspect=5) # 添加颜色标识

plt.show()
fig = plt.figure()
ax2 = fig.add_subplot(1, 1, 1, projection='3d')
surf2 = ax2.plot_surface(x, y, z-uh1, cmap='rainbow') # 误差图像，颜色映射
fig.colorbar(surf2, ax=ax2, shrink=0.5, aspect=5) # 添加颜色标识
ax2.set_xlim(0, 1)
ax2.set_ylim(0, 1)
ax2.set_zlim(0, 1)
plt.show()
if i < maxit:
mesh.uniform_refine()
1
the max error is 0.9506275647043294
the max error is 0.9274403759198202
the max error is 0.90481726916043
the max error is 0.8827445236782716
the max error is 0.8612087524670211
the max error is 0.8401968941440385
the max error is 0.8196962050299444
the max error is 0.7996942514208535
the max error is 0.780178902048576
the max error is 0.7611383207242166
the max error is 0.7425609591607093
the max error is 0.7244355499699358
the max error is 0.7067510998301808
the max error is 0.6894968828197817
the max error is 0.6726624339129286
the max error is 0.6562375426336726
the max error is 0.640212246864292
the max error is 0.6245768268042642
the max error is 0.609321799076177
the max error is 0.5944379109750086
the max error is 0.5799161348572853
the max error is 0.5657476626667193
the max error is 0.5519239005930004
the max error is 0.5384364638605101
the max error is 0.5252771716437914
the max error is 0.5124380421066957
the max error is 0.4999112875621961
the max error is 0.48768930974993296
the max error is 0.4757646952286273
the max error is 0.4641302108805686
the max error is 0.4527787995254512
the max error is 0.4417035756408989
the max error is 0.430897821187082
the max error is 0.42035498153289763
the max error is 0.4100686614812389
the max error is 0.4000326213909464
the max error is 0.39024077339308816
the max error is 0.3806871776992743
the max error is 0.37136603899976794
the max error is 0.3622717029492071
the max error is 0.35339865273780846
the max error is 0.34474150574597207
the max error is 0.3362950102802592
the max error is 0.32805404238876323
the max error is 0.3200136027539442
the max error is 0.31216881366103966
the max error is 0.3045149160402175
the max error is 0.2970472665806735
the max error is 0.2897613349149265
the max error is 0.2826527008716009
the max error is 0.2757170517950329
the max error is 0.2689501799300742
the max error is 0.2623479798705073
the max error is 0.2559064460695252
the max error is 0.24962167041076522
the max error is 0.2434898398384261
the max error is 0.23750723404502794
the max error is 0.23167022321541653
the max error is 0.22597526582563976
the max error is 0.2204189064953671
the max error is 0.21499777389254415
the max error is 0.2097085786890136
the max error is 0.2045481115658665
the max error is 0.1995132412673074
the max error is 0.19460091270185909
the max error is 0.18980814508975252
the max error is 0.1851320301553807
the max error is 0.1805697303637166
the max error is 0.17611847719963125
the max error is 0.17177556948906259
the max error is 0.16753837176102138
the max error is 0.16340431264943944
the max error is 0.1593708833338917
the max error is 0.15543563601824306
the max error is 0.15159618244630468
the max error is 0.14785019245359166
the max error is 0.14419539255430913
the max error is 0.1406295645627087
the max error is 0.13715054424797712
the max error is 0.13375622002184484
the max error is 0.13044453165811687
the max error is 0.12721346904334851
the max error is 0.1240610709579103
the max error is 0.12098542388670175
the max error is 0.11798466085879378
the max error is 0.11505696031529489
the max error is 0.11220054500475485
the max error is 0.10941368090543724
the max error is 0.1066946761738069
the max error is 0.10404188011859195
the max error is 0.10145368219980377
the max error is 0.09892851105210154
the max error is 0.0964648335319146
the max error is 0.09406115378774127
the max error is 0.09171601235306281
the max error is 0.08942798526131834
the max error is 0.08719568318241167
the max error is 0.08501775058021654
the max error is 0.08289286489057879
the max error is 0.08081973571930923
the max error is 0.07879710405968698
the max error is 0.07682374152899596
the max error is 0.0748984496236289
the max error is 0.07302005899231434
the max error is 0.07118742872701767
the max error is 0.06939944567109357
the max error is 0.06765502374426324
the max error is 0.06595310328401471
the max error is 0.06429265040301957
the max error is 0.06267265636218144
the max error is 0.06109213695893445
the max error is 0.059550131930418804
the max error is 0.05804570437117629
the max error is 0.05657794016500772
the max error is 0.05514594743064949
the max error is 0.053748855980935595
the max error is 0.052385816795111784
the max error is 0.05105600150398826
the max error is 0.04975860188761272
the max error is 0.04849282938516408
the max error is 0.047257914616765984
the max error is 0.046053106916933095
the max error is 0.04487767387936337
the max error is 0.04373090091280574
the max error is 0.042612090807727876
the max error is 0.04152056331352738
the max error is 0.04045565472602408
代码
文本
[36]
print(range(3))
range(0, 3)
代码
文本
python
python
点个赞吧
推荐阅读
公开
5/10 曹一帆 202105755309
python
python
bohre001cd
更新于 2024-07-07
公开
5/10
python
python
bohr7c365e
发布于 2024-05-18