21.3 The Finite Element Method
249
for ni=1:n
for nj=1:ni
zwi = polyint(conv(Cli(ni,:),Cli(nj,:)));
% x^0 *
Clint0(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
Clint0(nj,ni) = Clint0(ni,nj);
...
% x^1 *
...
% x^2 *
%
Derivation
zwi = polyint(conv(Cliabl(ni,:),Cliabl(nj,:)));%x^0
Clabint0(ni,nj) = polyval(zwi,1)-polyval(zwi,0);
Clabint0(nj,ni) = Clabint0(ni,nj);
...
% x^1
...
% x^2
end
end
%% FE
quadratic spacing
ml = 1:m;
% index FE
h0 = rmax/(m-1)^2;
% stepsize
h = (2 * ml-1) * h0;
r0 = (ml-1).^2 * h0;
%% Computation of the integrals
S = h(1) * (r0(1)^2 * Clint0 + 2 * h(1) * r0(1) * Clint1 + ...
h(1)^2 * Clint2);
S = S(:);
%
T = (r0(1)^2 * Clabint0 + 2 * h(1) * r0(1) * Clabint1 + ...
h(1)^2 * Clabint2)/h(1);
T = T(:);
%
V = -2 * h(1) * (r0(1) * Clint0 + h(1) * Clint1);
V = V(:);
%
for k=2:m
Sneu = h(k) * (r0(k)^2 * Clint0+2 * h(k) * r0(k) * Clint1+ ...
h(k)^2 * Clint2);
Sneu = Sneu(:);
S(end) = S(end)+Sneu(1);
Sneu(1) = [];
S = [S;Sneu];
...
T
...
V
end
%% Eigenvalue equation
H = T + V;
%% Indexgymnastics
...
% creating sparse matrices
S = sparse(iz,is,S);
H = sparse(iz,is,H);
figure, spy(H), shg
%% solving the generalized eigenvalue problem
E = eigs(H,S,25,-1);%’sa’)% 25 eigenvalues
E = E/2;
% use format rat