目标方程
我们需要解以下方程：
u
t
−
∇
⋅
(
c
∇
u
)
=
f
在
Ω
×
[
0
,
T
]
,
u_t - \nabla \cdot (c \nabla u) = f \quad \text{在} \; \Omega \times [0, T],
ut−∇⋅(c∇u)=f在Ω×[0,T],
其中，
∇
\nabla
∇ 表示梯度运算符，
u
t
u_t
ut 表示
u
u
u 关于时间
t
t
t 的偏导数。
u = g 在 ∂ Ω × [ 0 , T ] , u = g \quad \text{在} \; \partial \Omega \times [0, T], u=g在∂Ω×[0,T],
u = u 0 在 t = 0 时和 u ∈ Ω . u = u_0 \quad \text{在} \; t = 0 \; \text{时和} \; u \in \Omega. u=u0在t=0时和u∈Ω.
这里 Ω \Omega Ω 是2D区域， f ( x , y , t ) f(x,y,t) f(x,y,t) 和 c ( x , y , t ) c(x,y,t) c(x,y,t) 是已知的, u 0 ( x , y ) u_0 (x,y) u0(x,y) 给已经给出的初值函数， u ( x , y , t ) u(x,y,t) u(x,y,t)是未知函数。
弱格式
方程两边同时乘
v
(
x
,
y
)
v(x,y)
v(x,y)，
u
t
v
−
∇
⋅
(
c
∇
u
)
v
=
f
v
u_tv - \nabla \cdot (c \nabla u)v = fv
utv−∇⋅(c∇u)v=fv
∫
Ω
u
t
v
d
x
d
y
−
∫
Ω
∇
⋅
(
c
∇
u
)
v
d
x
d
y
=
∫
Ω
f
v
d
x
d
y
.
\int_{\Omega} u_t v \, dx \, dy - \int_{\Omega} \nabla \cdot (c \nabla u) v \, dx \, dy = \int_{\Omega} f v \, dx \, dy.
∫Ωutvdxdy−∫Ω∇⋅(c∇u)vdxdy=∫Ωfvdxdy.
u u u 叫做trial function， v v v 叫做test function。
然后使用格林公式
∫
Ω
∇
⋅
(
c
∇
u
)
v
d
x
d
y
=
∫
∂
Ω
(
c
∇
u
⋅
n
)
v
d
s
−
∫
Ω
c
∇
u
⋅
∇
v
d
x
d
y
,
\int_{\Omega} \nabla \cdot (c \nabla u) v \, dx \, dy = \int_{\partial \Omega} (c \nabla u \cdot \mathbf{n}) v \, ds - \int_{\Omega} c \nabla u \cdot \nabla v \, dx \, dy,
∫Ω∇⋅(c∇u)vdxdy=∫∂Ω(c∇u⋅n)vds−∫Ωc∇u⋅∇vdxdy,
因为这里是狄利克雷边界条件，我们选择
v
(
x
,
y
)
=
0
v (x, y)=0
v(x,y)=0 在边界上。这是因为，在边界上的
u
u
u值已经给出来了，所以不需要再去求解边界上的
u
u
u了。 这里可以去看我的b站讲解。 积分酱-有限元学习 1D-边界处理
我也有手写的笔记，在我的资源里可以下载，手写笔记-有限元学习 1D
现在我们可以得到简化后的形式：
∫
Ω
u
t
v
d
x
d
y
+
∫
Ω
c
∇
u
⋅
∇
v
d
x
d
y
=
∫
Ω
f
v
d
x
d
y
.
\int_{\Omega} u_t v \, dx \, dy + \int_{\Omega} c \nabla u \cdot \nabla v \, dx \, dy = \int_{\Omega} f v \, dx \, dy.
∫Ωutvdxdy+∫Ωc∇u⋅∇vdxdy=∫Ωfvdxdy.
这里
u
u
u 和
v
v
v 都属于 Sobolev空间。
定义：
H
1
(
0
,
T
;
H
1
(
Ω
)
)
=
{
v
(
t
,
⋅
)
∣
∂
v
∂
t
(
t
,
⋅
)
∈
H
1
(
Ω
)
,
∀
t
∈
[
0
,
T
]
}
H^1(0, T; H^1(\Omega)) = \{ v(t, \cdot) \;|\; \frac{\partial v}{\partial t}(t, \cdot) \in H^1(\Omega), \; \forall t \in [0, T] \}
H1(0,T;H1(Ω))={v(t,⋅)∣∂t∂v(t,⋅)∈H1(Ω),∀t∈[0,T]}
弱格式：
找到
u
∈
H
1
(
0
,
T
;
H
1
(
Ω
)
)
u \in H^1(0, T; H^1(\Omega))
u∈H1(0,T;H1(Ω)),使得：
∫
Ω
u
t
v
d
x
d
y
+
∫
Ω
c
∇
u
⋅
∇
v
d
x
d
y
=
∫
Ω
f
v
d
x
d
y
.
\int_{\Omega} u_t v \, dx \, dy + \int_{\Omega} c \nabla u \cdot \nabla v \, dx \, dy = \int_{\Omega} f v \, dx \, dy.
∫Ωutvdxdy+∫Ωc∇u⋅∇vdxdy=∫Ωfvdxdy.
成立，对于任意的
v
∈
H
0
1
(
Ω
)
v \in H_0^1(\Omega)
v∈H01(Ω),
H
0
1
(
Ω
)
H_0^1(\Omega)
H01(Ω)是指边界上为0的
H
1
(
Ω
)
H^1(\Omega)
H1(Ω)。
Galerkin格式
Galerkin格式的基本思路是：用有限维空间去近似无限维空间。
半离散
考虑有限元空间 设
U
h
=
s
p
a
n
{
φ
j
}
j
=
1
N
b
⊂
H
1
(
Ω
)
U_h = span \{\varphi_j \}_{j=1}^{N_b} \subset H^1(\Omega)
Uh=span{φj}j=1Nb⊂H1(Ω),其中，
{
φ
j
}
j
=
1
N
b
\{ \varphi_j \}_{j=1}^{N_b}
{φj}j=1Nb是全局有限元基函数的集合，
N
b
N_b
Nb是基函数的个数。
(这里基函数的定义和形式 在何晓明老师的课件中给出 何老师的基函数定义和形式)
定义
U
h
0
U_h^0
Uh0 为包含
U
h
U_h
Uh 中在 Dirichlet 边界上取值为 0 的函数的空间：
U
h
0
=
{
v
∈
U
h
∣
v
=
0
on the Dirichlet boundary
}
.
U_h^0 = \{ v \in U_h \mid v = 0 \text{ on the Dirichlet boundary} \}.
Uh0={v∈Uh∣v=0 on the Dirichlet boundary}.
我们需要找到
u
h
∈
H
1
(
0
,
T
;
U
h
)
u_h \in H^1(0, T; U_h)
uh∈H1(0,T;Uh)使得
∫
Ω
u
h
t
v
h
d
x
d
y
+
∫
Ω
c
∇
u
h
⋅
∇
v
h
d
x
d
y
=
∫
Ω
f
v
h
d
x
d
y
\int_{\Omega} u_{ht} v_h \, dx \, dy + \int_{\Omega} c \nabla u_h \cdot \nabla v_h \, dx \, dy = \int_{\Omega} f v_h \, dx \, dy
∫Ωuhtvhdxdy+∫Ωc∇uh⋅∇vhdxdy=∫Ωfvhdxdy
对所有
v
h
∈
U
h
0
v_h \in U_h^0
vh∈Uh0 都成立。
考虑
u
h
∈
H
1
(
0
,
T
;
U
h
)
u_h \in H^1(0, T; U_h)
uh∈H1(0,T;Uh)，且
U
h
=
s
p
a
n
{
φ
j
}
j
=
1
N
b
U_h = \mathrm{span} \{ \varphi_j \}_{j=1}^{N_b}
Uh=span{φj}j=1Nb，我们可以写成下面的形式
u
h
(
x
,
y
,
t
)
=
∑
j
=
1
N
b
u
j
(
t
)
φ
j
(
x
,
y
)
,
u_h(x, y, t) = \sum_{j=1}^{N_b} u_j(t) \varphi_j(x, y),
uh(x,y,t)=j=1∑Nbuj(t)φj(x,y),
其中
u
j
(
t
)
u_j(t)
uj(t)是待求的时间依赖系数，
j
=
1
,
2
,
…
,
N
b
j = 1, 2, \ldots, N_b
j=1,2,…,Nb。
如果我们能够为 u j ( t ) u_j(t) uj(t)建立一个线性代数系统，并求解这个系统，那么我们可以得到有限元解 u h u_h uh。
因此，我们选择测试函数$v_h = \varphi_i $（其中
i
=
1
,
…
,
N
b
i = 1, \ldots, N_b
i=1,…,Nb。那么我们有：
∫
Ω
(
∑
j
=
1
N
b
u
j
(
t
)
φ
j
)
t
φ
i
d
x
d
y
+
∫
Ω
c
∇
(
∑
j
=
1
N
b
u
j
(
t
)
φ
j
)
⋅
∇
φ
i
d
x
d
y
=
∫
Ω
f
φ
i
d
x
d
y
,
i
=
1
,
…
,
N
b
.
\int_{\Omega} \left( \sum_{j=1}^{N_b} u_j(t) \varphi_j \right)_t \varphi_i \, dx \, dy + \int_{\Omega} c \nabla \left( \sum_{j=1}^{N_b} u_j(t) \varphi_j \right) \cdot \nabla \varphi_i \, dx \, dy = \int_{\Omega} f \varphi_i \, dx \, dy, \quad i = 1, \ldots, N_b.
∫Ω(j=1∑Nbuj(t)φj)tφidxdy+∫Ωc∇(j=1∑Nbuj(t)φj)⋅∇φidxdy=∫Ωfφidxdy,i=1,…,Nb.
这等价于：
∑
j
=
1
N
b
u
j
′
(
t
)
[
∫
Ω
φ
j
φ
i
d
x
d
y
]
+
∑
j
=
1
N
b
u
j
(
t
)
[
∫
Ω
c
∇
φ
j
⋅
∇
φ
i
d
x
d
y
]
=
∫
Ω
f
φ
i
d
x
d
y
,
i
=
1
,
…
,
N
b
.
\sum_{j=1}^{N_b} u_j'(t) \left[ \int_{\Omega} \varphi_j \varphi_i \, dx \, dy \right] + \sum_{j=1}^{N_b} u_j(t) \left[ \int_{\Omega} c \nabla \varphi_j \cdot \nabla \varphi_i \, dx \, dy \right] = \int_{\Omega} f \varphi_i \, dx \, dy, \quad i = 1, \ldots, N_b.
j=1∑Nbuj′(t)[∫Ωφjφidxdy]+j=1∑Nbuj(t)[∫Ωc∇φj⋅∇φidxdy]=∫Ωfφidxdy,i=1,…,Nb.
这里
φ
i
\varphi_i
φi 不依赖于时间
t
t
t.
定义 the stiffness matrix
A
(
t
)
=
[
a
i
j
]
i
,
j
=
1
N
b
=
[
∫
Ω
c
∇
φ
j
⋅
∇
φ
i
d
x
d
y
]
i
,
j
=
1
N
b
A(t) = [a_{ij}]_{i,j=1}^{N_b}= \left[ \int_{\Omega} c \nabla \varphi_j \cdot \nabla \varphi_i \, dx \, dy \right]_{i,j=1}^{N_b}
A(t)=[aij]i,j=1Nb=[∫Ωc∇φj⋅∇φidxdy]i,j=1Nb
定义 the mass matrix
M
=
[
m
i
j
]
i
,
j
=
1
N
b
=
[
∫
Ω
φ
j
φ
i
d
x
d
y
]
i
,
j
=
1
N
b
M = [m_{ij}]_{i,j=1}^{N_b}=\left[ \int_{\Omega} \varphi_j \varphi_i \, dx \, dy \right]_{i,j=1}^{N_b}
M=[mij]i,j=1Nb=[∫Ωφjφidxdy]i,j=1Nb
定义 the load vector
b
~
(
t
)
=
[
b
i
]
i
=
1
N
b
=
[
∫
Ω
f
φ
i
d
x
d
y
]
i
=
1
N
b
\tilde{b}(t) = [b_{i}]_{i=1}^{N_b} = \left[ \int_{\Omega} f \varphi_i \, dx \, dy \right]_{i=1}^{N_b}
b~(t)=[bi]i=1Nb=[∫Ωfφidxdy]i=1Nb
定义 the unknown vector
X
~
(
t
)
=
[
u
j
(
t
)
]
j
=
1
N
b
\tilde{X}(t) = [u_{j}(t)]_{j=1}^{N_b}
X~(t)=[uj(t)]j=1Nb
现在我们就设置好了线性系统：
M
X
′
~
(
t
)
+
A
(
t
)
X
~
(
t
)
=
b
~
(
t
)
M \tilde{X'}(t) + A(t)\tilde{X}(t) = \tilde{b}(t)
MX′~(t)+A(t)X~(t)=b~(t)
全离散
回顾一下有限差分法： IVP （Initial value problem）
y
′
(
t
)
=
f
(
t
,
y
(
t
)
)
(
a
≤
t
≤
b
)
,
y
(
a
)
=
g
a
y ′(t) = f (t, y (t)) (a ≤ t ≤ b), y (a) = g_a
y′(t)=f(t,y(t))(a≤t≤b),y(a)=ga
将时间区域
[
a
,
b
]
[a,b]
[a,b]等分，为
M
m
M_m
Mm个，步长为
δ
t
\delta t
δt,
t
m
=
a
+
m
δ
t
tm=a+m\delta t
tm=a+mδt。
θ
θ
θ格式:
y
m
+
1
−
y
m
δ
t
=
θ
f
(
t
m
+
1
,
y
m
+
1
)
+
(
1
−
θ
)
f
(
t
m
,
y
m
)
,
y
m
=
y
(
t
m
)
\frac{y^{m+1} - y^m}{\delta t} = \theta f(t_{m+1}, y^{m+1}) + (1 - \theta) f(t_m, y^m), y^{m} = y(t_m)
δtym+1−ym=θf(tm+1,ym+1)+(1−θ)f(tm,ym),ym=y(tm)
θ
=
0
θ=0
θ=0: 向前欧拉格式
θ
=
1
θ=1
θ=1: 向后欧拉格式
θ
=
1
/
2
θ=1/2
θ=1/2: Crank-Nicolson格式
求解
M
X
′
~
(
t
)
+
A
(
t
)
X
~
(
t
)
=
b
~
(
t
)
M \tilde{X'}(t) + A(t)\tilde{X}(t) = \tilde{b}(t)
MX′~(t)+A(t)X~(t)=b~(t) 我们用有限差分法求解这个ODE（ordinary differential equatio常微分方程）。
对应的
θ
θ
θ格式就是：
M
X
~
m
+
1
−
X
~
m
δ
t
+
θ
A
(
t
m
+
1
)
X
~
m
+
1
+
(
1
−
θ
)
A
(
t
m
)
X
~
m
=
θ
b
~
(
t
m
+
1
)
+
(
1
−
θ
)
b
~
(
t
m
)
,
m
=
0
,
…
,
M
m
−
1.
M \frac{\tilde{X}^{m+1} - \mathbf{\tilde{X}}^m}{\delta t} + \theta A(t_{m+1}) \tilde{X}^{m+1} +(1- \theta)A(t_{m}) \tilde{X}^{m} =\theta \tilde{b}(t_{m+1}) + (1 - \theta) \mathbf{\tilde{b}}(t_m), \quad m = 0, \ldots, M_m - 1.
MδtX~m+1−X~m+θA(tm+1)X~m+1+(1−θ)A(tm)X~m=θb~(tm+1)+(1−θ)b~(tm),m=0,…,Mm−1.
这里
y
=
X
y=X
y=X
重新整理一下：
[
M
δ
t
+
θ
A
(
t
m
+
1
)
]
X
~
m
+
1
=
M
δ
t
X
~
m
+
θ
b
~
(
t
m
+
1
)
+
(
1
−
θ
)
b
~
(
t
m
)
−
(
1
−
θ
)
A
(
t
m
)
X
~
m
[\frac{M}{\delta t} + \theta A(t_{m+1}) ]\tilde{X}^{m+1}= \frac{M}{\delta t} \tilde{X}^{m}+ \theta \tilde{b}(t_{m+1}) + (1 - \theta) \mathbf{\tilde{b}}(t_m) -(1- \theta)A(t_{m}) \tilde{X}^{m}
[δtM+θA(tm+1)]X~m+1=δtMX~m+θb~(tm+1)+(1−θ)b~(tm)−(1−θ)A(tm)X~m
化简一下：
A
^
m
+
1
=
M
δ
t
+
θ
A
(
t
m
+
1
)
\hat{A}^{m+1} = \frac{M}{\delta t} + \theta A(t_{m+1})
A^m+1=δtM+θA(tm+1)
b
^
m
+
1
=
M
δ
t
X
~
m
+
θ
b
~
(
t
m
+
1
)
+
(
1
−
θ
)
b
~
(
t
m
)
−
(
1
−
θ
)
A
(
t
m
)
X
~
m
\hat{b}^{m+1} = \frac{M}{\delta t} \tilde{X}^{m}+ \theta \tilde{b}(t_{m+1}) + (1 - \theta) \mathbf{\tilde{b}}(t_m) -(1- \theta)A(t_{m}) \tilde{X}^{m}
b^m+1=δtMX~m+θb~(tm+1)+(1−θ)b~(tm)−(1−θ)A(tm)X~m
我们就可以得到迭代格式：
A
^
m
+
1
X
~
m
+
1
=
b
^
m
+
1
\hat{A}^{m+1} \tilde{X}^{m+1} = \hat{b}^{m+1}
A^m+1X~m+1=b^m+1
算法:
形成网格信息 P and T .
组装M矩阵 .
初始化 X_0
时间迭代:
FOR m = 0, · · · , Mm − 1:
tmp1 = (m + 1)dt;
tm = mdt;
组装矩阵A(tmp1) 和 A(tm)
组装向量b(tm+1) 和 b(tm)
得到A_hat和 b_hat
处理边界条件（A_hat，b_hat）;
求解迭代格式 X_mp1 = A_hat\b_hat
END
如果时间无关的组装部分，也可以在时间迭代前组装起来，这样可以节约计算时间。
function [err] = FE_solver_2D_second_order_parabolic_equation(element_type,start_t,end_t,dt,theta,domain,h,basis_type,basis_type_trial, ...
basis_type_test,Gauss_type)
format shorte
%% ----------------------------------------
% basis_type = 201 : 2D linear
% basis_type = 202 : 2D Qudratic
% P : information matrix consisting of the coordinates of all mesh nodes.
% T: information matrix consisting of the global node indices of the mesh nodes of all the mesh elements.
% Pb : an information matrix consisting of the coordinates of all finite element nodes
% Tb : an information matrix consisting of the global node indices of the finite element nodes of all the mesh elements.
% der_x_trail_A %r
% der_y_trail_A %s
% der_x_test_A %p
% der_y_test_A %q
% der_x_test_b %p
% der_y_test_b %q
%% ================= obtain Global info P T N Gauss_type element_type ===============
% Get P, T, N for global info
PT = generate_2D_PT(domain, h, element_type); % (domain, h,element_type)
N = size(PT.T,2);
Global.T = PT.T;
Global.P = PT.P;
Global.N = N;
Global.element_type=element_type;
Global.Gauss_type = Gauss_type;
%% ================= obtain Basis info Pb Tb Nb ===============
PbTb_Trail = generate_2D_PbTb(domain, h, basis_type_trial,element_type);
N_lb_trial =size(PbTb_Trail.Tb,1);
PbTb_Test = generate_2D_PbTb(domain, h, basis_type_test,element_type);
N_lb_test =size(PbTb_Test.Tb,1);
Nb_trail = size(PbTb_Trail.Pb,2);
Nb_test =size(PbTb_Test.Pb,2);
%% ------------------Obtain all functions ---------------------
all_funs = All_functions;
%% ================= Pre der info for matrix and vector =================
[der_A,der_b]=get_all_der_A_b;
%% ------------------Obtain A stiff matrix ---------------------
% trial and test are same; For A1,A2,A3,A4,A5,A6,A7,A8,A9
Basis.Nb_test= Nb_test;
Basis.Nb_trail= Nb_trail;
Basis.matrix_size= [Basis.Nb_test, Basis.Nb_trail];
Basis.Pb_trail = PbTb_Trail.Pb;
Basis.Pb_test = PbTb_Test.Pb;
Basis.Tb_trail = PbTb_Trail.Tb;
Basis.Tb_test = PbTb_Test.Tb;
Basis.N_lb_trial = N_lb_trial;
Basis.N_lb_test = N_lb_test;
Basis.basis_type_trial = basis_type_trial;
Basis.basis_type_test = basis_type_test;
A1 = assemble_matrix_2D(all_funs.fun_c,Global,Basis,der_A.der_A1);
A2 = assemble_matrix_2D(all_funs.fun_c,Global,Basis,der_A.der_A2);
A =A1+A2;
%% -----------------Obtain M matrix-------------------
M = assemble_matrix_2D(all_funs.fun_c_M,Global,Basis,der_A.der_M);
X_old = generate_initial_vector(all_funs.fun_initial,Basis.Pb_trail);
%A = sparse (matrix_size(1), matrix_size(2));
%A_tlide\b_tlide
number_of_time_step=(end_t-start_t)/dt;
A_tlide= M./dt +theta.*A;
temp = M./dt - (1-theta) .* A;
for m = 0:number_of_time_step-1
% X_old: Xm;
tm = start_t + m*dt;
tmp1= start_t + (m+1)*dt;
%% ------------------Obtain b load vector ---------------------
b_tm= assmble_vector_2D_time(all_funs.fun_f,Global,Basis,der_b.der_b1,tm);
b_tmp1= assmble_vector_2D_time(all_funs.fun_f,Global,Basis,der_b.der_b1,tmp1);
b_tlide= theta * b_tmp1 + (1-theta)*b_tm+temp*X_old;
%% ------------------ Handle Boundary conditions ---------------------
[boundaryedges,boundarynodes] = generate_boundary_info(element_type,basis_type,domain,h);
[A_tlide,b_tlide] = treat_Dirichlet_boundary_1_variable_time(all_funs.fun_Diri,tmp1, A_tlide, b_tlide, boundarynodes, Basis.Pb_test);
%% ------------------ Get Solution --------------
solution = A_tlide\b_tlide;
X_old = solution;
end
%% --------------------- compute error--------------------------%%
t = start_t:dt:end_t;
err = get_all_errors_time(solution,all_funs,Global,Basis,t);
结果：
h
1
h
2
L
∞
error
Convergence
L
∞
L
2
error
Convergence
L
2
H
1
semi-error
Convergence
H
1
0.1250
0.1250
9.8704
e
−
02
5.0853
e
−
02
1.2865
e
+
00
0.1250
0.0625
2.5483
e
−
02
1.9536
e
+
00
1.2871
e
−
02
1.9822
e
+
00
6.4214
e
−
01
1.0025
e
+
00
0.1250
0.0312
6.4745
e
−
03
1.9767
e
+
00
3.2279
e
−
03
1.9954
e
+
00
3.2092
e
−
01
1.0007
e
+
00
\begin{array}{cccccccc} \hline h_1 & h_2 &L_\infty \text{ error} & \text{Convergence } L_\infty & L_2 \text{ error} & \text{Convergence } L_2 & H^1 \text{ semi-error} & \text{Convergence } H^1 \\ \hline 0.1250 & 0.1250 & 9.8704e-02 & & 5.0853e-02 & & 1.2865e+00& \\ 0.1250 & 0.0625 & 2.5483e-02 & 1.9536e+00 & 1.2871e-02 & 1.9822e+00 & 6.4214e-01 & 1.0025e+00\\ 0.1250 & 0.0312 & 6.4745e-03 & 1.9767e+00 & 3.2279e-03 & 1.9954e+00 & 3.2092e-01 &1.0007e+00\\ \hline \end{array}
h10.12500.12500.1250h20.12500.06250.0312L∞ error9.8704e−022.5483e−026.4745e−03Convergence L∞1.9536e+001.9767e+00L2 error5.0853e−021.2871e−023.2279e−03Convergence L21.9822e+001.9954e+00H1 semi-error1.2865e+006.4214e−013.2092e−01Convergence H11.0025e+001.0007e+00
收敛阶的计算公式：
Convergence = log(error(1:end-1) ./ error(2:end)) ./ log(h(1:end-1) ./ h(2:end))
扫一扫
686
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
