有限元课程实验

Gypsophila
  2024-10-07 18:48:57 2024-10-07 18:48:57 Created   2024-10-07 21:24:11 2024-10-07 21:24:11 Updated    

问题一：分片线性有限元方法求解一维Poisson方程

使用等距网格剖分的分片线性有限元方法求解如下带有Dirichlet边界条件的一维Poisson方程：

令右端项为

分别在网格内点数为时进行数值求解。此问题的真实解为

计算和范数误差：

最后，使用

计算该方法在和范数下的收敛阶数。

数值求解

Step 1:将微分问题转化为变分问题
将泛定方程两侧同时与某测试函数进行内积，假设具有足够好的光滑性，分部积分可得

由于，应当选取的解空间为，要求测试函数，于是，故上式对应的变分问题为：

其中

Step 2:由变分问题构造有限元问题
选取解空间的一个有限维子空间，在其中寻找有限元解，此时有限元问题为：

这里要求使用分片线性有限元方法，所以令

其中为网格单元，是的一个等距剖分，网格剖分的步长为，是至多一阶多项式函数空间。

Step 3:组装全局刚度矩阵和全局载荷向量
由于是有限维空间，并且不难发现是的一组基，其中

于是有限元解可表示为

其中是待定系数。此时有限元问题等价于

上式可以写作矩阵形式：

其中是刚度矩阵，是载荷向量。为了求解得到，首先需要装配刚度矩阵和载荷向量。注意到，因此当时，，故是一个三对角矩阵。计算可知

因此

另一方面，由于，直接计算可知

对和积分可得

所以

于是负荷向量为

求解可以得到，最后代入的表达式即可得到数值解。

Step 4:计算误差
为了计算和范数需要使用数值积分，这里使用复合梯形公式进行数值积分：

通过使用足够小的可以保证数值积分误差远小于有限元解本身误差，因此该方法可以得到可接受的和范数误差近似。根据范数的定义，有

在每一个单元上进行数值积分：

其中是一个线性函数

于是的弱导数在上为常数

因此可以类似地计算为

其中每个单元上的值近似为

将这些单元上的误差项求和并开方即得和范数误差。

数值结果

Something is Wrong
准确解以及 得到的有限元解
误差和收敛阶数
N L2 error order H1 error order
10 1.41e-03 - 4.90e-02 -
20 3.86e-04 1.8659 2.57e-02 0.9336
40 1.01e-04 1.9305 1.31e-02 0.9650
80 2.60e-05 1.9646 6.65e-03 0.9833
160 6.57e-06 1.9821 3.35e-03 0.9905
320 1.65e-06 1.9910 1.67e-03 0.9987

程序代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Homework of FEM, Week 1                                       %
% -u'' = f, u(0)=u(1)=0;                                        %
% equal-size grid, linear interpolation                         %
% f(x) = -2*cos(x)+(x-1)*sin(x)                                 %
% Exact solution: u(x) = (x-1)*sin(x)                           %
% Input: N (number of interior points)                          %
% Output:L2, H1 (L2 error and H1 error), error order            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fem_hw1
N = 10*2.^(0:5);
L2 = zeros(length(N), 1); L2order = L2;
H1 = L2; H1order = H1;
fprintf('N\t L^2 error\t order\t\t H1 error\t order\n')
for i = 1:length(N)
    [L2(i), H1(i)] = femsol(N(i));
    if i>1
        L2order(i) = abs(log(L2(i)/L2(i-1))/log(2));
        H1order(i) = abs(log(H1(i)/H1(i-1))/log(2));
        fprintf('%d\t %.2e\t %.4f\t\t %.2e\t %.4f\n',...
            N(i),L2(i),L2order(i),H1(i),H1order(i))
    else
        fprintf('%d\t %.2e\t -\t\t %.2e\t -\n', N(i),L2(i),H1(i))
    end
end
% disp([N', L2, L2order, H1, H1order])

function [L2,H1] = femsol(N)        % N: number of interior points
% Parameters
h = 1/(N+1);        % step size
x = 0:h:1;          % gird

% Stiff Matrix
d = 1/h * ones(1,N);
A = diag(2*d);
A(2:N+1:N^2-N) = -1/h;              % lower subdiagonal elements
A(N+1:N+1:N^2-1) = -1/h;            % upper subdiagonal elements

% Load vector
F1 = @(x) -2*cos(x)+x*cos(x)-(x^2-2)*cos(x)-sin(x);
F2 = @(x) (1-x)*cos(x)-sin(x);
L = zeros(N, 1);
for j = 2:N+1
    L(j-1) = 1/h * (F1(x(j))-F1(x(j-1)) - x(j-1)*(F2(x(j))-F2(x(j-1))))...
        - 1/h * (F1(x(j+1))-F1(x(j)) - x(j+1)*(F2(x(j+1))-F2(x(j))));
    % L(j-1) = 1/h*(x(j-1)-1)*(sin(x(j))-sin(x(j-1)))...
    %     + 1/h*(x(j+1)-1)*(sin(x(j+1))-sin(x(j)));
end

% Solve linear system AU=L
U = A\L;

% Finite element solution
U = [0; U; 0];

%Plot
if N==20
u = @(x) (x-1).*sin(x) ;      % Exact solution
plot(x,u(x),'LineWidth',1)
hold on
plot(x, U, 'o', 'LineWidth',1)               % Finite element solution
hold off
legend('Exact Solution', 'FEM Solution')
xlabel('x')
ylabel('u')
end

% Error
L2 = L2error(U,x);
H1 = H1error(U,x);

function I = trapezoid(v, s)    % Numerically integration by Trapezoid's rule 
n = length(v);
v(2:n-1) = 2*v(2:n-1);
I = sum(v)*s/2;

function L2 = L2error(U,x)
L2 = 0;
s = 1e-5;                       % Intergration step
u = @(y) (y-1).*sin(y);
for i = 2:length(U)
    phi = @(y) U(i)*(y-x(i-1))/(x(i)-x(i-1)) + U(i-1)*(y-x(i))/(x(i-1)-x(i));
    E = @(y) (phi(y) - u(y)).^2;
    L2 = L2 + trapezoid(E(x(i-1):s:x(i)), s);
end
L2 = sqrt(L2);

function H1 = H1error(U,x)
H1 = L2error(U,x).^2;
s = 1e-5;                       % Intergration step
du = @(y) sin(y)+(y-1).*cos(y);
for i = 2:length(U)
    dphi = (U(i)-U(i-1))/(x(i)-x(i-1));
    E = @(y) (dphi - du(y)).^2;
    H1 = H1 + trapezoid(E(x(i-1):s:x(i)), s);
end
H1 = sqrt(H1);

问题二：局部刚度矩阵和局部载荷向量

使用等距网格剖分的分片线性有限元方法，构造局部刚度矩阵和局部负载向量求解与问题一相同的带有Dirichlet边界条件的一维Poisson方程：

令右端项为

分别在网格内点数为时进行数值求解。此问题的真实解为

计算和范数误差：

最后，使用

计算该方法在和范数下的收敛阶数。

数值求解

基本流程与问题一相同，只是在计算刚度矩阵和载荷向量时，使用局部刚度矩阵和局部载荷向量进行计算，因此在第三步中，组装刚度矩阵和载荷向量的过程应该按照如下方式修改。

首先需要指明各单元上的局部基函数为

使用如下线性变换

可得到标准单元上的基函数

同时，由于可以被紧致嵌入，因此全局基函数应当至少连续，所以这里选取相应的全局基函数为

现在我们可以使用这些基函数构造刚度矩阵和载荷向量。借助之前定义的线性变换，将各单元上的问题拉回到标准单元上可得

而有限元问题为

其中

现在将带入有限元问题可得

其中第一行只与第单元有关，而第二行仅与第单元有关。类似地，令可得

上式第一行只与第单元有关，第二行仅与第单元有关。将以上两式中与有关的项抽离出来可得

其中左侧的系数矩阵为需要的局部刚度矩阵，右侧的载荷向量为局部载荷向量。使用之前的线性变换，可知局部刚度矩阵为

直接积分可知局部载荷向量为

其中

使用这些局部刚度矩阵和局部载荷向量，可以组装全局刚度矩阵和全局载荷向量


以及

其中 。最后求解线性方程组 即可得到数值解。

数值结果

由于只是组装全局刚度矩阵和全局载荷向量时使用了局部刚度矩阵和局部载荷向量，全局刚度矩阵与负载向量与问题一相同，因此数值结果与问题一相同。

程序代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Homework of FEM, Week 2                                       %
% -u'' = f, u(0)=u(1)=0;                                        %
% equal-size grid, linear interpolation                         %
% f(x) = -2*cos(x)+(x-1)*sin(x)                                 %
% Exact solution: u(x) = (x-1)*sin(x)                           %
% Remark: Construct local stiff matrix and local load vecter    %
% Input: N (number of interior points)                          %
% Output:L2, H1 (L2 error and H1 error), error order            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fem_hw2
N = 10*2.^(0:5);
L2 = zeros(length(N), 1); L2order = L2;
H1 = L2; H1order = H1;
fprintf('N\t L^2 error\t order\t\t H1 error\t order\n')
for i = 1:length(N)
    [L2(i), H1(i)] = femsol(N(i));
    if i>1
        L2order(i) = abs(log(L2(i)/L2(i-1))/log(2));
        H1order(i) = abs(log(H1(i)/H1(i-1))/log(2));
        fprintf('%d\t %.2e\t %.4f\t\t %.2e\t %.4f\n',...
            N(i),L2(i),L2order(i),H1(i),H1order(i))
    else
        fprintf('%d\t %.2e\t -\t\t %.2e\t -\n', N(i),L2(i),H1(i))
    end
end
% disp([N', L2, L2order, H1, H1order])

function [L2,H1] = femsol(N)        % N: number of interior points
% Parameters
h = 1/(N+1);                        % step size
x = 0:h:1;                          % gird

% Stiff Matrix
A = zeros(N+2);                     % Initialization
K = 1/h * [1,-1;-1,1];              % Local Stiff Matrix
for j=2:N+2
    A(j-1:j,j-1:j)=A(j-1:j,j-1:j)+K;
end
A = A(2:N+1,2:N+1);                 % Gobal Stiff Matrix


% Load vector
F1 = @(x) -2*cos(x)+x*cos(x)-(x^2-2)*cos(x)-sin(x);
F2 = @(x) (1-x)*cos(x)-sin(x);
F = zeros(N+2, 1);                  % Initialization
for j = 2:N+2
    % Local load vector
    L = 1/h * [-(F1(x(j))-F1(x(j-1)) - x(j)*(F2(x(j))-F2(x(j-1))));
               (F1(x(j))-F1(x(j-1)) - x(j-1)*(F2(x(j))-F2(x(j-1))))];
    F(j-1:j) = F(j-1:j) + L;
    % L(j-1) = 1/h * (F1(x(j))-F1(x(j-1)) - x(j-1)*(F2(x(j))-F2(x(j-1))))...
    %     - 1/h * (F1(x(j+1))-F1(x(j)) - x(j+1)*(F2(x(j+1))-F2(x(j))));
end
F = F(2:N+1);                       % Gobal load vector

% Solve linear system AU=L
U = A\F;

% Finite element solution
U = [0; U; 0];

%Plot
if N==20
u = @(x) (x-1).*sin(x) ;                    % Exact solution
plot(x,u(x),'LineWidth',1)
hold on
plot(x, U, 'o', 'LineWidth',1)              % Finite element solution
hold off
legend('Exact Solution', 'FEM Solution')
xlabel('x')
ylabel('u')
end

% Error
L2 = L2error(U,x);
H1 = H1error(U,x);

function I = trapezoid(v, s)    % Numerically integration
n = length(v);
v(2:n-1) = 2*v(2:n-1);
I = sum(v)*s/2;

function L2 = L2error(U,x)
L2 = 0;
s = 1e-5;                       % Intergration step
u = @(y) (y-1).*sin(y);
for i = 2:length(U)
    phi = @(y) U(i)*(y-x(i-1))/(x(i)-x(i-1)) + U(i-1)*(y-x(i))/(x(i-1)-x(i));
    E = @(y) (phi(y) - u(y)).^2;
    L2 = L2 + trapezoid(E(x(i-1):s:x(i)), s);
end
L2 = sqrt(L2);

function H1 = H1error(U,x)
H1 = L2error(U,x).^2;
s = 1e-5;                       % Intergration step
du = @(y) sin(y)+(y-1).*cos(y);
for i = 2:length(U)
    dphi = (U(i)-U(i-1))/(x(i)-x(i-1));
    E = @(y) (dphi - du(y)).^2;
    H1 = H1 + trapezoid(E(x(i-1):s:x(i)), s);
end
H1 = sqrt(H1);

问题三：分片二次有限元方法求解一维Poisson方程

数值求解

按照与问题一相同的思路，只是在分片线性有限元方法的基础上，将插值函数改为二次多项式，这将导致局部刚度矩阵和局部载荷向量由二阶变为三阶，本质没有发生改变。

Step 1:将微分问题转化为变分问题
将泛定方程两侧同时与某测试函数进行内积，假设具有足够好的光滑性，分部积分可得

由于，应当选取的解空间为，要求测试函数，于是，故上式对应的变分问题为：

其中

Step 2:由变分问题构造有限元问题
选取解空间的一个有限维子空间，在其中寻找有限元解，此时有限元问题为：

这里要求使用分段连续二阶有限元方法，所以令

其中为网格单元，是的一个等距剖分，网格剖分的步长为，是至多二阶多项式函数空间。

Step 3:组装全局刚度矩阵和全局载荷向量
为了在每个单元上都是用二次多项式插值，需要一个额外的插值节点，这里选取，在每个单元上使用二次多项式插值，于是每个单元上的局部基函数为

使用如下线性变换

可得到标准单元上的基函数

同时，由于可以被紧致嵌入，因此全局基函数应当至少连续，所以这里选取相应的全局基函数为

现在我们可以使用这些基函数构造刚度矩阵和载荷向量。

首先，使用上述线性变换可以得到

而有限元问题为

其中

现在将带入有限元问题可得

其中第一行只与第单元有关，而第二行仅与第单元有关。类似地，令可得

最后令可得

现在将上面三个等式写成矩阵形式

其中左侧的系数矩阵为需要的局部刚度矩阵，右侧的载荷向量为局部载荷向量。使用之前的线性变换，可知局部刚度矩阵为

直接积分可知局部载荷向量为

其中

之后使用上述局部刚度矩阵和局部载荷向量组装全局刚度矩阵和全局载荷向量可得


以及

其中 。最后解线性方程组 即可得到数值解。

Step 4:误差计算
为了计算和范数需要使用数值积分，这里使用复合辛普森公式进行数值积分：

通过使用足够小的可以保证数值积分误差远小于有限元解本身误差，因此该方法可以得到可接受的和范数误差近似。根据范数的定义，有

在每一个单元上进行数值积分，之后可以类似地计算为

将这些单元上的误差项求和并开方即得和范数误差。

数值结果

Something is Wrong
准确解以及 得到的有限元解
误差和收敛阶数
N L2 error order H1 error order
10 4.25e-06 - 3.03e-04 -
20 6.11e-07 2.7967 8.32e-05 1.8641
40 8.21e-08 2.8952 2.18e-05 1.9300
80 1.07e-08 2.9468 5.59e-06 1.9648

程序代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Homework of FEM, Week 3                                       %
% -u'' = f, u(0)=u(1)=0;                                        %
% equal-size grid, quadratic interpolation                      %
% f(x) = -2*cos(x)+(x-1)*sin(x)                                 %
% Exact solution: u(x) = (x-1)*sin(x)                           %
% Remark: Construct local stiff matrix and local load vecter    %
% Input: N (number of interior points)                          %
% Output:L2, H1 (L2 error and H1 error), error order            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fem_hw3
N = 10*2.^(0:3);
L2 = zeros(length(N), 1); L2order = L2;
H1 = L2; H1order = H1;
fprintf('N\t L^2 error\t order\t\t H1 error\t order\n')
for i = 1:length(N)
    [L2(i), H1(i)] = femsol(N(i));
    if i>1
        L2order(i) = abs(log(L2(i)/L2(i-1))/log(2));
        H1order(i) = abs(log(H1(i)/H1(i-1))/log(2));
        fprintf('%d\t %.2e\t %.4f\t\t %.2e\t %.4f\n',...
            N(i),L2(i),L2order(i),H1(i),H1order(i))
    else
        fprintf('%d\t %.2e\t -\t\t %.2e\t -\n', N(i),L2(i),H1(i))
    end
end
% disp([N', L2, L2order, H1, H1order])

function [L2,H1] = femsol(N)        % N: number of interior points
% Parameters
h = 1/(N+1);                        % step size
x = 0:h:1;                          % gird (N+2 points)

A = zeros(2*N+3);
F = zeros(2*N+3,1);
F0 = @(y) (1-y).*cos(y)-sin(y);
F1 = @(y) (1-y).*y.*cos(y)-sin(y);
F2 = @(y) (-2+2*y+y.^2-y.^3).*cos(y)+(-2-2*y+y.^2).*sin(y);

for j = 1:N+1
    K = 1/3/h*[7, -8, 1; -8, 16, -8; 1, -8, 7];                     % Local
    A(2*j-1:2*j+1, 2*j-1:2*j+1) = A(2*j-1:2*j+1, 2*j-1:2*j+1) + K;

    b1 = 2/h^2*(...
        F2(x(j+1))-(2*x(j+1)-h/2)*F1(x(j+1))+x(j+1)*(x(j+1)-h/2)*F0(x(j+1))-...
        (F2(x(j))-(2*x(j+1)-h/2)*F1(x(j))+x(j+1)*(x(j+1)-h/2)*F0(x(j)))...
        );
    b2 = -4/h^2*(...
        F2(x(j+1))-(2*x(j+1)-h)*F1(x(j+1))+x(j+1)*(x(j+1)-h)*F0(x(j+1))-...
        (F2(x(j))-(2*x(j+1)-h)*F1(x(j))+x(j+1)*(x(j+1)-h)*F0(x(j)))...
        );
    b3 = 2/h^2*(...
        F2(x(j+1))-(2*x(j)+h/2)*F1(x(j+1))+x(j)*(x(j)+h/2)*F0(x(j+1))-...
        (F2(x(j))-(2*x(j)+h/2)*F1(x(j))+x(j)*(x(j)+h/2)*F0(x(j)))...
        );
    b = [b1;b2;b3];                                                 % Local
    F(2*j-1:2*j+1) = F(2*j-1:2*j+1) + b;
end
A = A(2:end-1,2:end-1);
F = F(2:end-1);

% Solve linear system AU=L
U = A\F;

% Finite element solution
U = [0; U; 0];

% Plot
x = 0:1/(N+1)/2:1;
if N==20
u = @(x) (x-1).*sin(x) ;                    % Exact solution
plot(x,u(x),'LineWidth',2)
hold on
plot(x, U, '-o', 'LineWidth',1)             % Finite element solution
hold off
legend('Exact Solution', 'FEM Solution')
xlabel('x')
ylabel('u')
end

% Error
L2 = L2error(U,x);
H1 = H1error(U,x);

function I = Simpson(v, s)      % Numerically integration by Simpson's rule 
v(2:2:end-1) = 4*v(2:2:end-1);
v(3:2:end-2) = 2*v(3:2:end-2);
I = sum(v)*s/6;

function L2 = L2error(U,x)
L2 = 0;
s = 1e-6;                       % Intergration step
u = @(y) (y-1).*sin(y);
N = (length(U)-3)/2;
for i = 1:N+1
    % I_i=[x_{i-1},x_i]: x(2*i-1),x(2*i),x(2*i+1)
    % U(2*i-1)-x_{i-1}, U(2*i)-x_{i-1/2},   U(2*i+1)-x_i
    phi = @(y) U(2*i+1)*(y-x(2*i)).*(y-x(2*i-1))...
        / (x(2*i+1)-x(2*i))/(x(2*i+1)-x(2*i-1))...
        + U(2*i)*(y-x(2*i-1)).*(y-x(2*i+1))...
        / (x(2*i)-x(2*i-1))/(x(2*i)-x(2*i+1))...
        + U(2*i-1)*(y-x(2*i)).*(y-x(2*i+1))...
        / (x(2*i-1)-x(2*i))/(x(2*i-1)-x(2*i+1));
    E = @(y) (phi(y) - u(y)).^2;
    L2 = L2 + Simpson(E(x(2*i-1):s:x(2*i+1)), s);
end
L2 = sqrt(L2);

function H1 = H1error(U,x)
H1 = L2error(U,x).^2;
s = 1e-6;                       % Intergration step
du = @(y) sin(y)+(y-1).*cos(y);
N = (length(U)-3)/2;
for i = 1:N+1
    dphi =@(y) U(2*i+1)*(y-x(2*i))/(x(2*i+1)-x(2*i))/(x(2*i+1)-x(2*i-1))...
        + U(2*i+1)*(y-x(2*i-1))/(x(2*i+1)-x(2*i))/(x(2*i+1)-x(2*i-1))...
        + U(2*i)*(y-x(2*i-1))/(x(2*i)-x(2*i-1))/(x(2*i)-x(2*i+1))...
        + U(2*i)*(y-x(2*i+1))/(x(2*i)-x(2*i-1))/(x(2*i)-x(2*i+1))...
        + U(2*i-1)*(y-x(2*i))/(x(2*i-1)-x(2*i))/(x(2*i-1)-x(2*i+1))...
        + U(2*i-1)*(y-x(2*i+1))/(x(2*i-1)-x(2*i))/(x(2*i-1)-x(2*i+1));
    E = @(y) (dphi(y) - du(y)).^2;
    H1 = H1 + Simpson(E(x(2*i-1):s:x(2*i+1)), s);
end
H1 = sqrt(H1);

问题四：分片线性有限元方法求解一维Poisson方程（Neumann边界条件）

数值求解

将泛定方程两侧同时与某测试函数进行内积，假设具有足够好的光滑性，分部积分可得

由于，所以上式变为

为了保证解的唯一性，需要令解空间为

于是为了从上述空间中寻找解，经过验证有限元解与真实解的关系可知这样得到的解满足

为了将解修正为真实解，需要添加一个常数修正项，于是最终的变分问题为

于是相应的有限元问题为

其中为分段连续线性有限元空间。于是与之前类似地，我们可以通过局部基函数得到局部刚度矩阵和局部载荷向量，之后组装全局刚度矩阵和全局载荷向量，最后解线性方程组即可得到数值解。此时全局刚度矩阵形如

全局载荷向量形如

计算得到


最终的数值解为

数值结果

Something is Wrong
准确解以及 得到的有限元解
误差和收敛阶数
N L2 error order H1 error order
10 6.61e-04 - 2.31e-02 -
20 1.84e-04 1.8448 1.22e-02 0.9158
40 4.85e-05 1.9247 6.29e-03 0.9605
80 1.24e-05 1.9631 3.19e-03 0.9811

程序代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Homework Two of FEM, Week 3                                   %
% -u'' = f, u'(0)=u'(1)=0;                                      %
% equal-size grid, Neumann Boundary Condition, linear           %
% f(x) = -12*x^2+12*x-2                                         %
% Exact solution: u(x) = x^2*(x-1)^2                            %
% Remark: Construct local stiff matrix and local load vecter    %
% Input: N (number of interior points)                          %
% Output:L2, H1 (L2 error and H1 error), error order            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fem_hw3_2
N = 10*2.^(0:3);
L2 = zeros(length(N), 1); L2order = L2;
H1 = L2; H1order = H1;
fprintf('N\t L^2 error\t order\t\t H1 error\t order\n')
for i = 1:length(N)
    [L2(i), H1(i)] = femsol(N(i));
    if i>1
        L2order(i) = abs(log(L2(i)/L2(i-1))/log(2));
        H1order(i) = abs(log(H1(i)/H1(i-1))/log(2));
        fprintf('%d\t %.2e\t %.4f\t\t %.2e\t %.4f\n',...
            N(i),L2(i),L2order(i),H1(i),H1order(i))
    else
        fprintf('%d\t %.2e\t -\t\t %.2e\t -\n', N(i),L2(i),H1(i))
    end
end
% disp([N', L2, L2order, H1, H1order])

function [L2,H1] = femsol(N)        % N: number of interior points
% Parameters
h = 1/(N+1);                        % step size
x = 0:h:1;                          % gird (N+2 points)

% Stiff Matrix
A = zeros(N+2);                     % Initialization
K = 1/h * [1,-1;-1,1];              % Local Stiff Matrix
for j=2:N+2
    A(j-1:j,j-1:j)=A(j-1:j,j-1:j)+K;
end
A = A(2:N+1,2:N+1);                 % Gobal Stiff Matrix
A = [A,h*ones(N,1);h*ones(1,N),0];

% Load vector
% f(x) = -12*x^2+12*x-2
F = zeros(N+2, 1);                  % Initialization
F1 = @(y) 3*y^4-4*y^3+y^2;
F2 = @(y) 4*y^3-6*y^2+2*y;
for j = 2:N+2
    % Local load vector
    K1 = 1/h*(F1(x(j))-F1(x(j-1))) - x(j)/h*(F2(x(j))-F2(x(j-1)));
    K2 = -1/h*(F1(x(j))-F1(x(j-1))) + x(j-1)/h*(F2(x(j))-F2(x(j-1)));
    F(j-1:j) = F(j-1:j) + [K1;K2];
end
F = F(2:N+1);                       
F = [F;0];                          % Gobal load vector

% Solve linear system AU=L
U = A\F;

% Finite element solution
U = [0; U(1:end-1); 0] + U(end)/2*(x.*(1-x))';

% Plot
x = 0:1/(N+1):1;
if N==20
u = @(x) (x-1).^2.*x.^2;                        % Exact solution
plot(x,u(x),'LineWidth',2)
hold on
plot(x, U, '-o', 'LineWidth',1)                 % Finite element solution
hold off
legend('Exact Solution', 'FEM Solution')
xlabel('x')
ylabel('u')
end

% Error
L2 = L2error(U,x);
H1 = H1error(U,x);

function I = trapezoid(v, s)    % Numerically integration by Trapezoid's rule 
n = length(v);
v(2:n-1) = 2*v(2:n-1);
I = sum(v)*s/2;

function L2 = L2error(U,x)
L2 = 0;
s = 1e-6;                       % Intergration step
u = @(y) ((y-1).*y).^2;
for i = 2:length(U)
    phi = @(y) U(i)*(y-x(i-1))/(x(i)-x(i-1)) + U(i-1)*(y-x(i))/(x(i-1)-x(i));
    E = @(y) (phi(y) - u(y)).^2;
    L2 = L2 + trapezoid(E(x(i-1):s:x(i)), s);
end
L2 = sqrt(L2);

function H1 = H1error(U,x)
H1 = L2error(U,x).^2;
s = 1e-6;                       % Intergration step
du = @(y) 2*y.*(1-3*y+2*y.^2);
for i = 2:length(U)
    dphi = (U(i)-U(i-1))/(x(i)-x(i-1));
    E = @(y) (dphi - du(y)).^2;
    H1 = H1 + trapezoid(E(x(i-1):s:x(i)), s);
end
H1 = sqrt(H1);
  • Title: 有限元课程实验
  • Author: Gypsophila
  • Created at : 2024-10-07 18:48:57
  • Updated at : 2024-10-07 21:24:11
  • Link: https://chenx.space/2024/10/07/FEMprojs/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments
On this page
有限元课程实验
  1. 问题一：分片线性有限元方法求解一维Poisson方程
    1. 数值求解
    2. 数值结果
    3. 程序代码
  2. 问题二：局部刚度矩阵和局部载荷向量
    1. 数值求解
    2. 数值结果
    3. 程序代码
  3. 问题三：分片二次有限元方法求解一维Poisson方程
    1. 数值求解
    2. 数值结果
    3. 程序代码
  4. 问题四：分片线性有限元方法求解一维Poisson方程（Neumann边界条件）
    1. 数值求解
    2. 数值结果
    3. 程序代码