1 - 问题域为 (Omega)
2 - 方程及其变分形式
[egin{aligned}
-
abla^2 u(oldsymbol{x}) &= f(oldsymbol{x}), quad &oldsymbol{x} in Omega \
u(oldsymbol{x}) &= u_D(oldsymbol{x}), quad & oldsymbol{x} in partialOmega
end{aligned}]
试函数及其积分：
[-int_{Omega}(
abla^2 u) v dx = int_{Omega} f v dx
]
变分形式：
由 ( abla(fcdot g) = abla f cdot abla g) 可得：（设(f= abla u) and (g=v)）
[
abla(
abla u cdot v) = (
abla^2u)cdot v +
abla u cdot
abla v
]
那么
[-int_{Omega}(
abla^2 u)vdx = int_{Omega}
abla u cdot
abla v dx - int_{Omega}
abla(
abla u cdot v) dx
]
其中，上式右侧第二项，可通过高斯定理进行化简为：
[int_{Omega}
abla(
abla u cdot v) dx = int_{partialOmega}
abla u cdot v cdot n ds
]
由于 (v) 在边界(partialOmega) 上为 0，因此，上式的值为零。那么，前一个积分公式为：
[-int_{Omega}(
abla^2 u)vdx = int_{Omega}
abla u cdot
abla v dx
]
最终得到变分形式：
[int_{Omega}
abla u cdot
abla v dx = int_{Omega} f v dx
]
从编程实现的角度讲，分四个步骤：
1 - 定义问题域
# create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)
2 - 定义边界条件
# define boundary condition
u_D = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]", 2)
def boundary(x, on_boundary):
return on_boundary
bc = DirchletBC(V, u_D, boundary)
3 - 定义问题
# define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
4 - 求解问题
# compute solution
u_s = Function(V)
solve(a == L, u_s, bc)
# plot solution
plto(u_s)