文章目录
- 引言
- 题目
- 补全方程
- 刚度矩阵
- 构造基底
- 边值条件非齐次
- 左边值条件非齐次
- 右边值条件非齐次
- 非齐次边值条件有限元方程
- 求数值解
- 直接求总刚度矩阵
- 先求单元刚度矩阵
引言
本文参考李荣华教授的《偏微分方程数值解法》一书
题目
对于非齐次第二边值问题
{
−
d
d
x
(
p
d
u
d
x
)
+
q
u
=
f
,
a
<
x
<
b
u
(
a
)
=
α
,
u
′
(
b
)
=
β
\begin{cases} -\dfrac{d}{dx}(p\dfrac{du}{dx})+qu=f,a<x<b\\ u(a)=\alpha,u'(b)=\beta \end{cases}
⎩
⎨
⎧−dxd(pdxdu)+qu=f,a<x<bu(a)=α,u′(b)=β
已知精确解
u
(
x
)
=
sin
5
x
+
x
3
(
2
−
x
)
+
2
u(x)=\sin5x+x^3(2-x)+2
u(x)=sin5x+x3(2−x)+2,并且
{
p
(
x
)
=
cos
x
q
(
x
)
=
x
\begin{cases} p(x)=\cos x\\ q(x)=x \end{cases}
{p(x)=cosxq(x)=x
求剖分区间数分别为16和32时的精确解并画图
补全方程
使用 S a g e M a t h SageMath SageMath进行符号运算,求出 u ′ ( x ) u'(x) u′(x)和 f ( x ) f(x) f(x),代码如下
x = var('x')
u = sin(5*x)+x^3*(2-x)+2
p = cos(x)
q = x
f = -derivative(p * derivative(u, x), x) + q * u
print(f)
print(derivative(u, x))
结果为:
-((x - 2)*x^3 - sin(5*x) - 2)*x + (6*(x - 2)*x + 6*x^2 + 25*sin(5*x))*cos(x) - (3*(x - 2)*x^2 + x^3 - 5*cos(5*x))*sin(x)
-3*(x - 2)*x^2 - x^3 + 5*cos(5*x)
于是原方程为
{
−
d
d
x
(
p
d
u
d
x
)
+
q
u
=
−
(
(
x
−
2
)
x
3
−
sin
5
x
−
2
)
x
+
(
6
(
x
−
2
)
x
+
6
x
2
+
25
sin
5
x
)
cos
(
x
)
−
(
3
(
x
−
2
)
x
2
+
x
3
−
5
cos
5
x
)
sin
x
,
a
<
x
<
b
u
(
a
)
=
sin
5
a
+
a
3
(
2
−
a
)
+
2
u
′
(
b
)
=
−
3
(
x
−
2
)
x
2
−
x
3
+
5
cos
5
x
\begin{cases} -\dfrac{d}{dx}(p\dfrac{du}{dx})+qu=-((x - 2)x^3 - \sin5x - 2)x +\\ (6(x - 2)x + 6x^2 + 25\sin 5x)\cos(x) - (3(x - 2)x^2 + x^3 - 5\cos 5x)\sin x,\\&a<x<b\\ u(a)=\sin5a+a^3(2-a)+2\\ u'(b)=-3(x - 2)x^2 - x^3 + 5\cos5x \end{cases}
⎩
⎨
⎧−dxd(pdxdu)+qu=−((x−2)x3−sin5x−2)x+(6(x−2)x+6x2+25sin5x)cos(x)−(3(x−2)x2+x3−5cos5x)sinx,u(a)=sin5a+a3(2−a)+2u′(b)=−3(x−2)x2−x3+5cos5xa<x<b
刚度矩阵
对于网格剖分
a
=
x
0
<
x
1
<
⋯
<
x
n
=
b
a=x_0<x_1<\cdots<x_n=b
a=x0<x1<⋯<xn=b
对于每个单元
I
i
=
[
x
i
−
1
,
x
i
]
I_i = [x_{i-1}, x_i]
Ii=[xi−1,xi],按照线性插值公式,有数值解
u
h
=
x
i
−
x
h
i
u
i
−
1
+
x
−
x
i
−
1
h
i
u
i
,
x
∈
I
i
,
i
=
1
,
2
,
⋯
,
n
u_h = \frac{x_i-x}{h_i}u_{i-1} + \frac{x-x_{i-1}}{h_i}u_i,x\in I_i,i=1,2,\cdots,n
uh=hixi−xui−1+hix−xi−1ui,x∈Ii,i=1,2,⋯,n
其中
h
i
=
x
i
−
x
i
−
1
h_i=x_i-x_{i-1}
hi=xi−xi−1
为使按段插值标准化,通常用仿射变换
ξ
=
x
−
x
i
−
1
h
i
\xi =\frac{x-x_{i-1}}{h_i}
ξ=hix−xi−1
把
I
i
I_i
Ii变到
ξ
\xi
ξ轴上的参考单元
[
0
,
1
]
[0,1]
[0,1],令
N
0
(
ξ
)
=
1
−
ξ
,
N
i
(
ξ
)
=
ξ
N_0(\xi)=1-\xi,N_i(\xi)=\xi
N0(ξ)=1−ξ,Ni(ξ)=ξ
则
u
h
=
N
0
(
ξ
)
u
i
−
1
+
N
1
(
ξ
)
u
i
u_h=N_0(\xi)u_{i-1}+N_1(\xi)u_i
uh=N0(ξ)ui−1+N1(ξ)ui
代入泛函
J
(
u
)
=
1
2
∫
a
b
(
p
u
′
2
+
q
u
2
−
2
f
u
)
d
x
J(u)=\frac{1}{2}\int_a^b (pu'^2+qu^2-2fu)dx
J(u)=21∫ab(pu′2+qu2−2fu)dx
令
KaTeX parse error: Undefined control sequence: \part at position 8: \frac{\̲p̲a̲r̲t̲ ̲J(u_h)}{\part u…
其中,
u
j
=
u
h
(
x
j
)
,
j
=
1
,
2
,
⋯
,
n
,
u
0
=
0
u_j=u_h(x_j),j=1,2,\cdots,n,u_0=0
uj=uh(xj),j=1,2,⋯,n,u0=0
{
a
j
−
1
,
j
=
∫
0
1
[
−
p
(
x
j
−
1
+
h
j
ξ
)
/
h
j
+
h
j
q
(
x
j
−
1
+
h
j
ξ
)
ξ
(
1
−
ξ
)
]
d
ξ
a
j
+
1
,
j
=
∫
0
1
[
−
p
(
x
j
+
h
j
+
1
ξ
)
/
h
j
+
1
+
h
j
+
1
q
(
x
j
+
h
j
+
1
ξ
)
ξ
(
1
−
ξ
)
]
d
ξ
a
j
j
=
∫
0
1
[
p
(
x
j
−
1
+
h
j
ξ
)
/
h
j
+
h
j
q
(
x
j
−
1
+
h
j
ξ
)
ξ
2
]
d
ξ
+
∫
0
1
[
p
(
x
j
+
h
j
+
1
ξ
)
/
h
j
+
1
+
h
j
+
1
q
(
x
j
+
h
j
+
1
ξ
)
(
1
−
ξ
)
2
]
d
ξ
b
j
=
∫
0
1
[
h
j
f
(
x
j
−
1
+
h
j
ξ
)
ξ
]
d
ξ
+
∫
0
1
[
h
j
+
1
f
(
x
j
+
h
j
+
1
ξ
)
(
1
−
ξ
)
]
d
ξ
⋯
(
∗
)
\begin{cases} a_{j-1,j}&=\int_0^1[-p(x_{j-1}+h_j\xi)/h_j+h_jq(x_{j-1}+h_j\xi)\xi(1-\xi)]d\xi\\ a_{j+1,j}&=\int_0^1[-p(x_{j}+h_{j+1}\xi)/h_{j+1}+h_{j+1}q(x_{j}+h_{j+1}\xi)\xi(1-\xi)]d\xi\\ a_{jj}&=\int_0^1[p(x_{j-1}+h_j\xi)/h_j+h_jq(x_{j-1}+h_j\xi)\xi^2]d\xi+\\ &\int_0^1[p(x_{j}+h_{j+1}\xi)/h_{j+1}+h_{j+1}q(x_{j}+h_{j+1}\xi)(1-\xi)^2]d\xi\\ b_{j}&=\int_0^1[h_jf(x_{j-1}+h_j\xi)\xi]d\xi+\int_0^1[h_{j+1}f(x_{j}+h_{j+1}\xi)(1-\xi)]d\xi \end{cases}\cdots(*)
⎩
⎨
⎧aj−1,jaj+1,jajjbj=∫01[−p(xj−1+hjξ)/hj+hjq(xj−1+hjξ)ξ(1−ξ)]dξ=∫01[−p(xj+hj+1ξ)/hj+1+hj+1q(xj+hj+1ξ)ξ(1−ξ)]dξ=∫01[p(xj−1+hjξ)/hj+hjq(xj−1+hjξ)ξ2]dξ+∫01[p(xj+hj+1ξ)/hj+1+hj+1q(xj+hj+1ξ)(1−ξ)2]dξ=∫01[hjf(xj−1+hjξ)ξ]dξ+∫01[hj+1f(xj+hj+1ξ)(1−ξ)]dξ⋯(∗)
于是可以得到总刚度矩阵
K
K
K,满足
K
u
=
b
Ku=b
Ku=b
发现
a
j
j
a_{jj}
ajj和
b
j
b_j
bj由不同的两项相加,于是可以得到单元刚度矩阵
K
(
i
)
K^{(i)}
K(i)
K
(
i
)
=
(
a
i
−
1
,
i
−
1
(
i
)
a
i
−
1
,
i
(
i
)
a
i
,
i
−
1
(
i
)
a
i
i
(
i
)
)
K^{(i)}= \begin{pmatrix} a^{(i)}_{i-1,i-1}&a^{(i)}_{i-1,i}\\ a^{(i)}_{i,i-1}&a^{(i)}_{ii} \end{pmatrix}
K(i)=(ai−1,i−1(i)ai,i−1(i)ai−1,i(i)aii(i))
其中
{
a
i
−
1
,
i
−
1
(
i
)
=
a
i
−
1
,
i
−
1
(
i
)
=
∫
0
1
[
−
p
(
x
i
−
1
+
h
i
ξ
)
/
h
i
+
h
i
q
(
x
i
−
1
+
h
i
ξ
)
ξ
(
1
−
ξ
)
]
d
ξ
a
i
−
1
,
i
−
1
(
i
)
=
∫
0
1
[
p
(
x
i
−
1
+
h
i
ξ
)
/
h
i
+
h
i
q
(
x
i
−
1
+
h
i
ξ
)
(
1
−
ξ
)
2
]
d
ξ
a
i
i
(
i
)
=
∫
0
1
[
p
(
x
i
+
h
i
+
1
ξ
)
/
h
i
+
1
+
h
i
+
1
q
(
x
i
+
h
i
+
1
ξ
)
ξ
2
]
d
ξ
⋯
(
∗
∗
)
\begin{cases} a^{(i)}_{i-1,i-1}=a^{(i)}_{i-1,i-1}=\int_0^1[-p(x_{i-1}+h_i\xi)/h_i+h_iq(x_{i-1}+h_i\xi)\xi(1-\xi)]d\xi\\ a^{(i)}_{i-1,i-1}=\int_0^1[p(x_{i-1}+h_i\xi)/h_i+h_iq(x_{i-1}+h_i\xi)(1-\xi)^2]d\xi\\ a^{(i)}_{ii}=\int_0^1[p(x_{i}+h_{i+1}\xi)/h_{i+1}+h_{i+1}q(x_{i}+h_{i+1}\xi)\xi^2]d\xi\\ \end{cases}\cdots (**)
⎩
⎨
⎧ai−1,i−1(i)=ai−1,i−1(i)=∫01[−p(xi−1+hiξ)/hi+hiq(xi−1+hiξ)ξ(1−ξ)]dξai−1,i−1(i)=∫01[p(xi−1+hiξ)/hi+hiq(xi−1+hiξ)(1−ξ)2]dξaii(i)=∫01[p(xi+hi+1ξ)/hi+1+hi+1q(xi+hi+1ξ)ξ2]dξ⋯(∗∗)
不难发现,总刚度矩阵可以由单元刚度矩阵
K
(
i
)
K^{(i)}
K(i)“组装”而成,即对于刚度矩阵
K
K
K中的元素
a
i
j
a_{ij}
aij和单元刚度矩阵
K
(
i
)
K^{(i)}
K(i)中的元素
a
i
j
(
i
)
a_{ij}^{(i)}
aij(i),有如下关系
a
i
j
=
{
a
i
,
i
−
1
(
i
)
,
j
=
i
−
1
a
i
i
(
i
)
+
a
i
i
(
i
+
1
)
,
j
=
i
a
i
,
i
+
1
(
i
+
1
)
,
j
=
i
+
1
0
,
∣
j
−
i
∣
≥
2
a_{ij}= \begin{cases} a_{i,i-1}^{(i)}, &j=i-1\\ a_{ii}^{(i)}+a_{ii}^{(i+1)}, &j=i\\ a_{i,i+1}^{(i+1)}, &j=i+1\\ 0, &|j-i|\ge 2 \end{cases}
aij=⎩
⎨
⎧ai,i−1(i),aii(i)+aii(i+1),ai,i+1(i+1),0,j=i−1j=ij=i+1∣j−i∣≥2
不难发现,总刚度矩阵
K
K
K为三对角矩阵
对应的,向量
b
b
b也可以进行“拆分”,令
{
f
i
−
1
(
i
)
=
∫
0
1
[
h
j
f
(
x
j
−
1
+
h
j
ξ
)
(
1
−
ξ
)
]
d
ξ
f
i
(
i
)
=
∫
0
1
[
h
j
+
1
f
(
x
j
+
h
j
+
1
ξ
)
ξ
]
d
ξ
\begin{cases} f^{(i)}_{i-1}=\int_0^1[h_jf(x_{j-1}+h_j\xi)(1-\xi)]d\xi\\ f^{(i)}_{i}=\int_0^1[h_{j+1}f(x_{j}+h_{j+1}\xi)\xi]d\xi \end{cases}
{fi−1(i)=∫01[hjf(xj−1+hjξ)(1−ξ)]dξfi(i)=∫01[hj+1f(xj+hj+1ξ)ξ]dξ
则有
b
1
=
f
1
(
1
)
+
f
1
(
2
)
,
b
1
=
f
2
(
2
)
+
f
2
(
3
)
,
⋯
,
b
1
=
f
n
(
n
)
⋯
(
∗
∗
∗
)
b_1=f_1^{(1)}+f_1^{(2)},b_1=f_2^{(2)}+f_2^{(3)},\cdots,b_1=f_n^{(n)}\cdots (***)
b1=f1(1)+f1(2),b1=f2(2)+f2(3),⋯,b1=fn(n)⋯(∗∗∗)
构造基底
再单元
I
i
,
I
i
+
1
I_i,I_{i+1}
Ii,Ii+1考察线性插值公式及
u
i
u_i
ui的系数,我们自然对每一节点
x
i
x_i
xi构造山形函数:
{
φ
i
(
x
)
=
{
1
+
x
−
x
i
h
i
,
x
i
−
1
≤
x
≤
x
i
1
−
x
−
x
i
h
i
+
1
,
x
i
≤
x
≤
x
i
+
1
0
,
o
t
h
e
r
w
i
s
e
i
=
1
,
2
,
⋯
,
n
−
1
φ
n
(
x
)
=
{
1
+
x
−
x
n
h
n
,
x
n
−
1
≤
x
≤
x
n
0
,
o
t
h
e
r
w
i
s
e
\begin{cases} \varphi_i(x)&=\begin{cases} 1+\dfrac{x-x_i}{h_i},&x_{i-1}\le x\le x_i\\ 1-\dfrac{x-x_i}{h_{i+1}},&x_{i}\le x\le x_{i+1}\\ 0,&otherwise \end{cases}\\ &i=1,2,\cdots,n-1\\ \varphi_n(x)&=\begin{cases} 1+\dfrac{x-x_n}{h_n},&x_{n-1}\le x\le x_n\\ 0,&otherwise \end{cases}\\ \end{cases}
⎩
⎨
⎧φi(x)φn(x)=⎩
⎨
⎧1+hix−xi,1−hi+1x−xi,0,xi−1≤x≤xixi≤x≤xi+1otherwisei=1,2,⋯,n−1=⎩
⎨
⎧1+hnx−xn,0,xn−1≤x≤xnotherwise
于是
u
h
(
x
)
=
∑
i
=
1
n
u
i
φ
i
(
x
)
,
u
i
=
u
h
(
x
i
)
u_h(x)=\sum\limits_{i=1}^{n}u_i\varphi_i(x),u_i=u_h(x_i)
uh(x)=i=1∑nuiφi(x),ui=uh(xi)
边值条件非齐次
左边值条件非齐次
u ( α ) = a u(\alpha)=a u(α)=a
则增加一基函数
φ
0
(
x
)
=
{
a
−
x
−
x
0
h
1
,
x
0
≤
x
≤
x
1
0
,
o
t
h
e
r
w
i
s
e
\varphi_0(x)=\begin{cases} a-\dfrac{x-x_0}{h_1},&x_0\le x\le x_1\\ 0,&otherwise \end{cases}
φ0(x)=⎩
⎨
⎧a−h1x−x0,0,x0≤x≤x1otherwise
将
u
h
u_h
uh表为
u
h
=
α
φ
0
(
x
)
+
∑
i
=
1
n
u
i
φ
i
(
x
)
u_h=\alpha \varphi_0(x)+\sum_{i=1}^n u_i\varphi_i(x)
uh=αφ0(x)+i=1∑nuiφi(x)
有限元方程为
∑
i
=
1
n
a
(
φ
i
,
φ
j
)
u
i
=
(
f
,
φ
j
)
−
α
a
(
φ
0
,
φ
j
)
,
j
=
1
,
2
,
⋯
,
n
\sum_{i=1}^n a(\varphi_i,\varphi_j)u_i=(f,\varphi_j)-\alpha a(\varphi_0,\varphi_j),j=1,2,\cdots,n
i=1∑na(φi,φj)ui=(f,φj)−αa(φ0,φj),j=1,2,⋯,n
实际上只修改了第一个方程的右端,因为当
j
=
2
,
3
,
⋯
,
n
j=2,3,\cdots,n
j=2,3,⋯,n时,
a
(
φ
0
,
φ
j
)
=
0
a(\varphi_0,\varphi_j)=0
a(φ0,φj)=0
右边值条件非齐次
则右端修改为
(
f
,
φ
j
)
+
p
(
b
)
β
φ
j
(
b
)
(f,\varphi_j)+p(b)\beta\varphi_j(b)
(f,φj)+p(b)βφj(b),有限元方程变为
∑
i
=
1
n
u
i
∫
a
b
[
p
φ
i
′
φ
j
′
+
q
φ
i
φ
j
]
d
x
=
(
f
,
φ
j
)
+
p
(
b
)
β
φ
j
(
b
)
,
j
=
1
,
2
,
⋯
,
n
\sum_{i=1}^n u_i\int_a^b[p\varphi_i'\varphi_j'+q\varphi_i\varphi_j]dx=(f,\varphi_j)+p(b)\beta\varphi_j(b),j=1,2,\cdots,n
i=1∑nui∫ab[pφi′φj′+qφiφj]dx=(f,φj)+p(b)βφj(b),j=1,2,⋯,n
实际上只修改了第
n
n
n个方程的右端项,因为当
j
=
1
,
2
,
⋯
,
n
−
1
j=1,2,\cdots,n-1
j=1,2,⋯,n−1时,
φ
j
(
b
)
=
0
\varphi_j(b)=0
φj(b)=0
非齐次边值条件有限元方程
综上,我们得到非齐次边值条件有限元方程
(
a
(
φ
1
,
φ
1
)
⋯
a
(
φ
n
,
φ
1
)
a
(
φ
1
,
φ
2
)
⋯
a
(
φ
n
,
φ
2
)
⋮
⋱
⋮
a
(
φ
1
,
φ
n
−
1
)
⋯
a
(
φ
n
,
φ
n
−
1
)
a
(
φ
1
,
φ
n
)
⋯
a
(
φ
n
,
φ
n
)
)
(
u
1
u
2
⋮
u
n
−
1
u
n
)
=
(
(
f
,
φ
1
)
−
α
a
(
φ
0
,
φ
1
)
(
f
,
φ
2
)
⋮
(
f
,
φ
n
−
1
)
(
f
,
φ
n
)
+
p
(
b
)
β
φ
n
(
b
)
)
\begin{pmatrix} a(\varphi_1,\varphi_1)&\cdots&a(\varphi_n,\varphi_1)\\ a(\varphi_1,\varphi_2)&\cdots&a(\varphi_n,\varphi_2)\\ \vdots&\ddots&\vdots\\ a(\varphi_1,\varphi_{n-1})&\cdots&a(\varphi_n,\varphi_{n-1})\\ a(\varphi_1,\varphi_n)&\cdots&a(\varphi_n,\varphi_n)\\ \end{pmatrix} \begin{pmatrix} u_1\\ u_2\\ \vdots\\ u_{n-1}\\ u_n \end{pmatrix} =\begin{pmatrix} (f,\varphi_1)-\alpha a(\varphi_0,\varphi_1)\\ (f,\varphi_2)\\ \vdots\\ (f,\varphi_{n-1})\\ (f,\varphi_n)+p(b)\beta\varphi_n(b) \end{pmatrix}
a(φ1,φ1)a(φ1,φ2)⋮a(φ1,φn−1)a(φ1,φn)⋯⋯⋱⋯⋯a(φn,φ1)a(φn,φ2)⋮a(φn,φn−1)a(φn,φn)
u1u2⋮un−1un
=
(f,φ1)−αa(φ0,φ1)(f,φ2)⋮(f,φn−1)(f,φn)+p(b)βφn(b)
求数值解
我们以区间 [ − 1 , 2 ] [-1,2] [−1,2]为例做均匀剖分
在求解数值解 u = ( u 1 , u 2 , ⋯ , u n ) T u=(u_1,u_2,\cdots,u_n)^T u=(u1,u2,⋯,un)T时,我们可以直接求总刚度矩阵,也可以先求单元刚度矩阵再进行“组装”
直接求总刚度矩阵
M a t l a b Matlab Matlab代码如下
% 直接求刚度矩阵
n = ;
a = -1;
b_ = 2;
A = zeros(n, n);
b = zeros(n, 1);
p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
p(1,i) = a + (b_-a)/n.*(i-1);
if i<n+1
I(1,i) = i;
I(2,i) = i+1;
end
end
u0 = @(x) sin(5.*x)+x.^3.*(2-x)+2;
du = @(x) -3.*(x - 2).*x.^2 - x.^3 + 5.*cos(5.*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
aii = @(x1, x2, h) integral(@(ksi) (px(x1+h.*ksi)/h + h.*qx(x1+h.*ksi).*ksi.^2 + px(x2+h.*ksi)/h + h.*qx(x2+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
ann = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
bi = @(x1, x2, h) integral(@(ksi) (h.*f(x1+h.*ksi).*ksi + h.*f(x2+h.*ksi).*(1-ksi)), 0, 1);
bn = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);
for i=1:1:n-1
x1 = p(1, I(1,i));
x2 = p(1, I(2,i));
h = x2 - x1;
A(I(1,i), I(1,i)) = aii(x1, x2, h);
b(I(1, i), 1) = bi(x1, x2, h);
end
for i=1:1:n-1
x1 = p(1, I(1,i+1));
x2 = p(1, I(2,i+1));
h = x2 - x1;
A(I(1,i), I(2,i)) = a23(x1, h);
A(I(2,i), I(1,i)) = a23(x1, h);
end
x1 = p(1, I(1,n));
x2 = p(1, I(2,n));
h = x2 - x1;
A(n, n) = ann(x1, h);
% b(n, 1) = bi(x1, x2, h);
b(n, 1) = bn(p(1, I(1,n)), h);
% 边值条件
x1 = p(1, I(1,1));
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);
dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
for j=1:1:length(x)
if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
h = p(1, I(2, i)) - p(1, I(1, i));
u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
h = p(I(2, i+1)) - p(I(1, i+1));
u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
end
end
end
for j=1:1:length(x)
if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
h = p(1, I(2, n)) - p(1, I(1, n));
u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
end
end
for j=1:1:length(x)
if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
h = p(1, I(2, n)) - p(1, I(1, n));
u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
end
end
u_ = sin(5.*x)+x.^3.*(2-x)+2;
% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'b', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'g', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;
当n=16
时结果如下
当n=32
时结果如下
先求单元刚度矩阵
M a t l a b Matlab Matlab代码如下
% 先求单元刚度矩阵
n = ;
a = -1;
b_ = 2;
A = zeros(n, n);
b = zeros(n, 1);
p = zeros(1, n+1);
I = zeros(2, n);
for i=1:1:n+1
p(1,i) = a + (b_-a)/n.*(i-1);
if i<n+1
I(1,i) = i;
I(2,i) = i+1;
end
end
u0 = @(x) sin(5*x)+x.^3.*(2-x)+2;
du = @(x) -3*(x - 2).*x.^2 - x.^3 + 5*cos(5*x);
alpha = u0(a);
beta = du(b_);
px = @(x) cos(x);
qx = @(x) x;
f = @(x) -((x - 2).*x.^3 - sin(5*x) - 2).*x + (6*(x - 2).*x + 6*x.^2 + 25*sin(5*x)).*cos(x) - (3*(x - 2).*x.^2 + x.^3 - 5*cos(5*x)).*sin(x);
a1 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*(1-ksi).^2), 0, 1);
a23 = @(x, h) integral(@(ksi) (-px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.*(1-ksi)),0, 1);
a4 = @(x, h) integral(@(ksi) (px(x+h.*ksi)/h + h.*qx(x+h.*ksi).*ksi.^2), 0, 1);
b1 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*(1-ksi)), 0, 1);
b2 = @(x, h) integral(@(ksi) (h.*f(x+h.*ksi).*ksi), 0, 1);
for i=1:1:n-1
x1 = p(1, I(1,i+1));
x2 = p(1, I(2,i+1));
h = x2 - x1;
a1_ = a1(x1, h);
a2_ = a23(x1, h);
a3_ = a23(x1, h);
a4_ = a4(x1, h);
A(I(1,i), I(1,i)) = A(I(1,i), I(1,i)) + a1_;
A(I(1,i), I(2,i)) = A(I(1,i), I(2,i)) + a2_;
A(I(2,i), I(1,i)) = A(I(2,i), I(1,i)) + a3_;
A(I(2,i), I(2,i)) = A(I(2,i), I(2,i)) + a4_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
A(I(1,1), I(1,1)) = A(I(1,1), I(1,1)) + a4(x1, h);
for i=2:1:n
x1 = p(1, I(1,i));
x2 = p(1, I(2,i));
h = x2 - x1;
b1_ = b1(x1, h);
b2_ = b2(x1, h);
b(I(1,i-1), 1) = b(I(1,i-1), 1) + b1_;
b(I(2,i-1), 1) = b(I(2,i-1), 1) + b2_;
end
x1 = p(1, I(1,1));
x2 = p(1, I(2,1));
h = x2 - x1;
b(I(1,1), 1) = b(I(1,1), 1) + b2(x1, h);
% 边值条件
b(I(1,1), 1) = b(I(1,1), 1) - alpha.*a23(x1, h);
b(I(2,n-1), 1) = b(I(2,n-1), 1) + px(b_).*beta;
u = A\b;
c = cond(A);
dx = 1/1000;
x = a:dx:b_;
u1 = zeros(1, length(x));
for i=1:1:n-1
for j=1:1:length(x)
if p(1, I(1,i))<=x(1, j) && x(1, j)<=p(1, I(2,i)) % x_{i-1} <= x <= x_i
h = p(1, I(2, i)) - p(1, I(1, i));
u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, i))) / h) * u(I(1, i), 1);
elseif p(1, I(1,i+1))<=x(1, j) && x(1, j)<=p(1, I(2,i+1)) % x_i <= x <= x_{i+1}
h = p(I(2, i+1)) - p(I(1, i+1));
u1(1, j) = u1(1, j) + (1 - (x(1, j) - p(I(1, i+1))) / h) * u(I(1, i), 1);
end
end
end
for j=1:1:length(x)
if p(1, I(1,n))<=x(1, j) && x(1, j)<=p(1, I(2,n)) % x_{n-1} <= x <= x_n
h = p(1, I(2, n)) - p(1, I(1, n));
u1(1, j) = u1(1, j) + ((x(1, j) - p(1, I(1, n))) / h) * u(I(1, n), 1);
end
end
for j=1:1:length(x)
if p(1, I(1,1))<=x(1, j) && x(1, j)<=p(1, I(2,1)) % x_0 <= x <= x_1
h = p(1, I(2, n)) - p(1, I(1, n));
u1(1, j) = u1(1, j) + alpha * (1 - (x(1, j) - p(1, I(1, 1))) / h);
end
end
u_ = sin(5.*x)+x.^3.*(2-x)+2;
% 绘制数值解和精确解的图像
figure(1);
plot(x, u1, 'r', 'LineWidth', 1.5);
hold on;
plot(x, u_, 'b', 'LineWidth', 1.5);
xlabel('x');
ylabel('u(x)');
legend('数值解', '精确解');
title('两点边值问题的有限元方法');
grid on;
hold off;
box off;
当n=16
时结果如下
当n=32
时结果如下