代码bug #1703
Description
关键在
errors = mesh.error(pde.solution, uh, errortype='all')
print(errors)
发现输出的三种误差的数据类型不一样，有时候不太抗压方便使用。
输出示例：
(np.float64(0.0130138678719538), array([0.01017449]), array([0.02505717]))
具体完整代码如下，方便调试
import matplotlib.pyplot as plt#用于绘制网格和结果可视化
from fealpy.backend import backend_manager as bm#FEALPy的计算后端管理模块，这里设置为NumPy
bm.set_backend('numpy')
from fealpy.mesh import UniformMesh#生成任意维维均匀网格
from fealpy.sparse import csr_matrix, spdiags#用于生成稀疏矩阵
from fealpy.fdm import LaplaceOperator
from fealpy.fdm import DirichletBC
from fealpy.solver import spsolve
from fealpy.model import PDEDataManager
class SinPDEData:
def domain(self):
""" 方程的定义域
"""
return [-1, 1]
def solution(self, p):
x = p[..., 0]
return bm.exp(-x**2) * (1 - x**2)#方程的真解
def geo_dimension(self):
return 1
def gradient(self, p):
x = p[..., 0]
return 2 * bm.exp(-x**2)* ( x**3-2*x)#方程的一阶导数
def source(self, p):
x = p[..., 0]
return 2 * bm.exp(-x**2)* (2* x**4-8*x**2+3)#方程的右端项
def dirichlet(self, p):#方程的边界条件函数，用于处理边界
x = p[..., 0]
return bm.exp(-x**2) * (1 - x**2)
pde = SinPDEData()#创建该微分方程
domain = pde.domain()#定义该方程的定义域
nx = 10
extent=[0, nx]*pde.geo_dimension()
mesh = UniformMesh(domain, extent)#生成网格
node = mesh.entity('node')#网格节点
NN=len(node)
uI = pde.solution(node)# 网格函数
A = LaplaceOperator(mesh=mesh).assembly()
data = bm.ones(NN)
I = spdiags(data, 0, NN, NN, format='csr')
A=A+2*I
f = pde.source(node)
bc = DirichletBC(mesh=mesh, gd=pde.dirichlet)
A, f = bc.apply(A, f)
uh = spsolve(A, f, "scipy")
errors = mesh.error(pde.solution, uh, errortype='all')
print(errors)