Figure 2.30:

Approximate solutions to \eqref {eq:MWRexample} using the finite element method with hat functions for N=6 and N=12. The exact solution also is shown.

Code for Figure 2.30

Text of the GNU GPL.

main.m

 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
% finite element solution of u''+u=-x,u(0)=u(1)=0
% linear elements
% MDG 1/4/13
N=5; %no basis functions phi_0 or phi_N so here N only includes interior elements
h=1/(N+1); %N+1 here corresponds to N in the book
A=zeros(N,N);
u=zeros(N,1);
x1=(0:h:1)';
b=zeros(N,1);
A(1,1)=2/3 *h-2/h;
A(2,1)=h/6+1/h;
A(N,N)=2/3 *h-2/h;
A(N-1,N)=h/6+1/h;
for i=2:N-1
    A(i,i)=2/3 *h-2/h;
    A(i-1,i)=h/6+1/h;
    A(i+1,i)=h/6+1/h;
end
for i=1:N
    b(i)=-i*h^2;
end
u= A \ b;
us=zeros(N+2,1); %vector of solution values including boundaries
us(1)=0;
us(N+2)=0;
for i=1:N
    us(i+1)=u(i);
end
us1 = us;

% exact solution
Nex=100;
hex=1/(Nex+1);
uex=zeros(Nex+2,1);
xex=0:hex:1;
uex=-xex+sin(xex)*csc(1);

N=11; %repeat but with larger N
h=1/(N+1);
A=zeros(N,N);
u=zeros(N,1);
x2=(0:h:1)';
b=zeros(N,1);
A(1,1)=2/3 *h-2/h;
A(2,1)=h/6+1/h;
A(N,N)=2/3 *h-2/h;
A(N-1,N)=h/6+1/h;
for i=2:N-1
    A(i,i)=2/3 *h-2/h;
    A(i-1,i)=h/6+1/h;
    A(i+1,i)=h/6+1/h;
end
for i=1:N
    b(i)=-i*h^2;
end
u = A \ b;
us=zeros(N+2,1); %vector of solution values including boundaries
us(1)=0;
us(N+2)=0;
for i=1:N
    us(i+1)=u(i);
end
us2 = us;
plot(xex, uex, '-', x1, us1, '-.', x2, us2,'--')
legend('Exact','N=6','N=12')

t1 = [xex', uex']; t2 = [x1, us1]; t3=[x2, us2];
save "FEGalerkin1.dat" t1 t2 t3