Page 433 - Applied Numerical Methods Using MATLAB

422 PARTIAL DIFFERENTIAL EQUATIONS
T
which, with u(x, y) = c ϕ(x, y), can be written as
∂ ∂ T T ∂ ∂ T

T
I = c ϕ ϕ c + c ϕ ϕ c
R ∂x ∂x ∂y ∂y

T
T
T
− g(x, y)c ϕϕ c + 2f (x, y)c ϕ dx dy (9.4.12)
The condition for this functional to be minimized with respect to c is
d ∂ ∂ T ∂ ∂ T
I = ϕ 2 ϕ c + ϕ 2 ϕ c
dc 2 R ∂x ∂x ∂y ∂y

T
− g(x, y)ϕ 2 ϕ c + f (x, y)ϕ 2 dx dy = 0 (9.4.13)

N s

≈ A 1 c 1 + A 2 c 2 + f(x s ,y s )ϕ 2,s S s = 0 (9.4.14)
s=1
See [R-1] for details.
The objectives of the MATLAB routines “fem_basis_ftn()”and
“fem_coef()” are to construct the basis function φ n,s (x, y)’s for each node
n = 1,...,N n and each subregion s = 1,...,N s and to get the coefficient vector
c of the solution (9.4.4) via Eq. (9.4.7) and the solution polynomial φ s (x, y)’s
via Eq. (9.4.6) for each subregion s = 1,...,N s , respectively.
Before going into a specific example of applying the FEM method to solve
a PDE, let us take a look at the basis (shape) function φ n (x, y) for each node
n = 1,...,N n , which is defined collectively for all of the (triangular) subregions
so that φ n is 1 only at node n, and 0 at all other nodes and can be generated by
the routine “fem_basis_ftn()”.

function p = fem_basis_ftn(N,S)
%p(i,s,1:3): coefficients of each basis ftn phi_i
% for s-th subregion(triangle)
%N(n,1:2):x&y coordinates of the n-th node
%S(s,1:3) : the node #s of the s-th subregion(triangle)
N_n = size(N,1); % the total number of nodes
N_s = size(S,1); % the total number of subregions(triangles)
for n = 1:N_n
for s = 1:N_s
for i = 1:3
A(i,1:3) = [1 N(S(s,i),1:2)];
b(i) = (S(s,i) == n); %The nth basis ftn is 1 only at node n.
end
pnt=A\b’;
for i=1:3, p(n,s,i) = pnt(i); end
end
end