Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
heat equation in 1D and finite element methods

Subject: heat equation in 1D and finite element methods

From: ritchie31

Date: 10 Dec, 2012 21:49:54

Message: 1 of 1

Hello forum.
I am solving the heat equation i 1d using FEM. It works fine for initial condition. However when I tried to run it for t>0 it does not work properly.

Here is my code, you can run it with this inputs and see
[Uh]=ParaFEM1D(1,0,1,10,7)


function [Uh, xnod]=ParaFEM1D(Tf,a, b, Nt, xnel)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This function shows how to set up a finite element solution
%% of a two-point BVP
%% -(ku')' =f, u(a)=ua, u(b)=ub
%% using Galerkin linear Finite Elements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% function sol = solfem (a,b,xnel)
%% k,f,ua,ub are provided as functions in code
%% numerical integration is used in computation of stiffness matrix and rhs
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% INPUT:
%% <a,b> = endpoints of interval
%% <xnel> = the number or location of elements
%% if xnel is a scalar, then xnel = number of elements,
%% if xnel is a vector, then it contains position of first node of each element
%% OUTPUT:
%% <sol> is the computed numerical solution at nodes
%% <xnod> position of all nodes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clf;
set(0,'DefaultLineLineWidth',3);
%%%%%%%%%%%%%%%%%%%%%%%%% set up grid
dt = Tf/Nt;
t = 0:dt:Tf;

if size(xnel) == 1 %% uniform grid
    
   nels = xnel;
   hel = (b-a)./nels;
   xel = a : hel : b-hel;

elseif size(xnel,1) == 1

   nels = size(xnel,2);
   xel=xnel;
   
   for i=1:nels-1
       
       hel(i) = xnel(i+1)-xnel(i);
   
   end;

else

   error('wrong number of elements/grid');

end;

%% set up uniform order of elements = 1 + degree of polynomial =
%% = number of degrees of freedom

ord = zeros(nels,1) + 2; %% type of elements: linear
maxord = max(ord);

%% number of nodes
nnodes = sum(ord-1)+1;

%% derive global indexing of nodes:
%% nod(i,j) is the global number of j'th node in element i
nod = zeros(nels,maxord); myel = zeros(nnodes,2);
n = 1;
for i = 1:nels
    for j = 1:ord(i)
        nod(i,j) = n;
        if j == 1
            myel(n,2) = i;
        elseif j == ord(i)
            myel(n,1) = i;
        else
            myel(n,1) = i;
            myel(n,2) = i;
        end;
        
        if j ~= ord(i)
            n = n+1;
        end
    end;
end;

%% xnod (i=1..nnodes): coordinates of node i
xnod = zeros(nnodes,1);
for i=1:nels-1
    h = xel(i+1)-xel(i);
    hi = h/(ord(i)-1);
    for j=1:ord(i)
        xnod (nod(i,j)) = xel(i) + hi*(j-1);
    end;
end;

i = nels;
h = b-xel(i);
hi=h/(ord(i)-1);

for j=1:ord(i)
    xnod (nod(i,j)) = xel(i) + hi*(j-1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%% set up numerical integration

%% set up quadrature parameters on the reference element (-1,1)
%% set up number of integration points nw, nodes xw, and weights w

if maxord == 1 %% exact for linears
      nw = 1;
      xw(1) = 0.;
      w(1) = 2.;
elseif maxord == 2 %% exact for cubics
      nw = 2;
      xw(1) = -1/sqrt(3); xw(2) = -xw(1);
      w(1) = 1; w(2) = 1;
elseif maxord == 3 %% exact for polynomials of degree 5
      nw = 3;
      xw(1) = -sqrt(3./5.); xw(2)=0.; xw(3) =- xw(1);
      w(1) = 5./9.; w(2)=8./9.; w(3)=w(1);
end;

%%%%%%%%%%%%%%%%%%%%%%%%% matrix and rhs of linear system
b_global = sparse(nnodes,nnodes); m_global = sparse(nnodes,nnodes);
gs_global = sparse(nnodes,nnodes); rhsf = zeros(nnodes,1);
for el = 1:nels
    x1 = xnod(nod(el,1)); %% left endpoint
    x2 = xnod(nod(el,2)); %% right endpoint
    dx = (x2-x1)/2.; %% Jacobian of transformation
    %% compute element stiffness matrix and load vector
    b_local = zeros(ord(el),ord(el)); %% element stiffness matrix
% gs_local = zeros(ord(el),ord(el));
    m_local = zeros(ord(el),ord(el));
% gs_local = zeros(ord(el),ord(el)); %% GS = B+dt*A
    f = zeros(ord(el),1); %% element load vector
    %nw is the order of Gaussian Integration
    for i = 1:nw
        x = x1 + (1 + xw(i))*dx; %% x runs in true element,
                                                %% xw runs in reference element
        [psi,dpsi] = shape(xw(i),ord); %% calculations on ref.element
% kval = feval(@kfun,x);
% fval = feval(@rhsfun,0,x);
        fval = feval(@exfun,0,x);
        f = f + fval * psi * w(i)*dx;
        
        b_local = b_local + (psi*psi')* w(i)*dx; %mass vector
        m_local = m_local + w(i)*(dt*(dpsi*dpsi')/dx/dx) * dx;%stiffness matrix
% gs_local = gs_local + w(i)*(dt*(dpsi*dpsi')/dx/dx) * dx + (psi*psi')* w(i)*dx ;
     end
    %% add the computed element stiffness matrix and load vector to
    %% the global matrix and vector
    rhsf(nod(el,:)) = rhsf(nod(el,:)) + f(:); % F = load vector
    b_global(nod(el,:),nod(el,:)) = b_global(nod(el,:),nod(el,:)) + b_local; % B = mass matrix
    m_global(nod(el,:),nod(el,:)) = m_global(nod(el,:),nod(el,:)) + m_local; % A = stiffness matrix
% gs_global(nod(el,:),nod(el,:)) = gs_global(nod(el,:),nod(el,:)) + gs_local;
end

%% impose Dirichlet boundary conditions
xa = exfun(0,a); xb = exfun(0,b) ;
Uh = zeros(Nt,nnodes);
for i=1:nnodes
% rhsf(i) = rhsf(i) - m_global(i,1)*xa - m_global(i,nnodes)*xb ;
     rhsf(i) = rhsf(i) - (dt*m_global(i,1)-b_global(i,1))*xa - ...
         (dt*m_global(i,nnodes)-b_global(i,nnodes))*xb ;
    m_global(i,1) =0; m_global(1,i)=0;
    b_global(i,1) =0; b_global(1,i)=0;
    m_global(i,nnodes)=0; m_global(nnodes,i)=0;
    b_global(i,nnodes)=0; b_global(nnodes,i)=0;
end

m_global(1,1) = 1; b_global(1,1) = 1; rhsf(1,1) = xa;
m_global(nnodes,nnodes)=1;b_global(nnodes,nnodes)=1; rhsf(nnodes,1) = xb;

for n=1:Nt
   fprintf('itetation %d of %d\n',n,Nt);
% k = t(n+1) - t(n); % time step;
   if n == 1
      Uh(n,:) = b_global\rhsf;
   else
% Uh(ntime,:) = Uh(ntime-1,:)*(gs_global\(dt*b_global));
      Uh(n,:) = (b_global+dt*m_global)\(b_global*Uh(n,:)'+dt*rhsf);
   end
  figure(1)
  exfun(dt*(n-1),xnod)'
  plot(xnod,exfun(dt*(n-1),xnod),'k',xnod,Uh(n,:),'r--');
  pause(0.5),legend('Exact solution','Numerical solution (linear FE)');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% end of algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [y,dy] = shape(x,n)

%% shape function on reference element (-1,1)
%% n = 2: linear
%% n = 3: quadratic (must be coded)

if n == 2

   y (1,:) = .5.*(1-x);
   y (2,:) = .5.*(1+x);
   dy (1,:) = -.5;
   dy (2,:) = .5;

end


function y = exfun(t, x)
     
y= sin(pi*x).*exp(-pi^2*t) + sin(2*pi*x).*exp(-4*pi^2*t);
% dy = pi*cos(pi*x).*exp(-pi^2*t) + 2*pi*cos(2*pi*x).*exp(-4*pi^2*t);

Please help

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us