maple_FEM_solution_basic.mws

> restart;
phi:=proc(x,L)
      if x>L then
         0
      else
         if x>=0 then 1-x else
            phi(-x,L)
         end if;
       end if;
end proc;

phiLocal:=proc(i,x,h,L)
    phi(x/h-(i-1),L)
end proc;

yapprox:=proc(x,c,h,n,L)
   local i;
   add(c[i]*phiLocal(i,x,h,L),i=1..n)
end proc;

L:=1; nPoints:=100; nElements:=nPoints-1;
q:=4; f:=4;
nShapeFunctions:=nPoints;
h:=evalf(L/nElements);
#plot([''phiLocal(i,x,h,L)'' $i=1..nPoints],x=0..L);
with(LinearAlgebra):

makeMatrix:=proc(q,h,f,nPoints)
local A,i,j;
A:=Matrix(nPoints);
for i from 1 to nPoints do
  for j from 1 to nPoints do
    if i=j then
       A[i,j]:=2+2/3 * q*h^2;
    else
       if j=i-1 or j=i+1 then
          A[i,j]:= -1+1/6 *q *h^2;
       end if;
    end if;
  end do;
end do;
A
end proc:


A:=makeMatrix(q,h,f,nPoints):
b:=Vector(nPoints):
b[1..-1]:=h^2*f:
c:=LinearSolve(A,b):
with(plots):

listplot(['yapprox(i*h,c,h,nPoints,L)' $i=0..nPoints-1]);

 

 

 

 

 

 

 

 

 

 

 

proc (x, L) if `<`(L, x) then 0 else if `<=`(0, x) then `+`(1, `-`(x)) else phi(`+`(`-`(x)), L) end if end if end proc
proc (i, x, h, L) phi(`+`(`/`(`*`(x), `*`(h)), `-`(i), 1), L) end proc
proc (x, c, h, n, L) local i; add(`*`(c[i], `*`(phiLocal(i, x, h, L))), i = 1 .. n) end proc
1
100
99
4
4
100
0.1010101010e-1
Plot_2d
 

>