Matlab code Collocation Method of the Resolution of the Poisson's equation with the WRM using global Shape Functions

From KratosWiki
Jump to: navigation, search
% Resolution of the Poisson's equation with the WRM 
% using global Shape Functions
% Collocation Method 
% domain between x0 - xL
% l= xL-x0
% Q constant over a part of the domain (x0-xf)
% k constant
% fi(x_0)=fi0
% fi(x_L)=fiL

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% analytical solution %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% parameters
x0=0;
xL=1;
xf=xL*0.5;
fi0=0;
fiL=0;
k=1;
Q=1;
Nx=101;
sep=(xL-x0)/(Nx-1);

xa=x0:sep:xf;
xb=xf+sep:sep:xL;
x=[xa xb];

Qa=Q*ones(size(xa));
Qb=zeros(size(xb));
Qt=[Qa Qb];

cte1a=(k/(xL-x0))*((fiL-fi0)-(Q/(2*k))*(x0^2-2*xL*xf+xf^2));
cte1b=-Q*xf+cte1a;
c1a=cte1a.*ones(size(xa));
c1b=cte1b.*ones(size(xb));
c1=[c1a c1b];

cte2a=((fi0*xL-fiL*x0)/(xL-x0))-(Q/(2*k))*(x0/(xL-x0))*(2*xL*xf+x0*xL-xf^2);
cte2b=cte2a+Q*(xf^2)/(2*k);
c2a=cte2a*ones(size(xa));
c2b=cte2b*ones(size(xb));
c2=[c2a c2b];

fi=-(Qt./(2*k)).*(x.^2)+(c1./k).*x+c2;
dfi=-(Qt./k).*x+(c1./k);
gfi=gradient(fi)./gradient(x);

kdfi=k.*dfi;

subplot(131),dibuja(x,fi,'x','fi','unknow'); 
subplot(132),dibuja(x,dfi,'x','dfi','gradient') 
subplot(133),dibuja(x,kdfi,'x','kdfi','flow')   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%  end of analytical solution %%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

l=xL-x0;
fit=fi;



% case 2: Fourier Series
% Kij= k (pi j / l)^2 sin(pi xi j / l)
% fi = Q(xi)

npoints=1;
xn(1)=x(round(length(x)/2));
Qn(1)=Q*0.5;
K(1,1)= k*(pi/l)^2*sin(pi*xn(1)/l);
f(1)=Qn(1);

a=inv(K)*f;

fi2a=a*sin(pi*x/l);
fit=[fit;fi2a];

figure(3)
mdibuja(x,fit,'x','fi','analytical fi vs numerical fi');

npoints=2;
xn(1)=x(round(length(x)/4));
xn(2)=x(round(3*length(x)/4));

Qn(1)=Q;
Qn(2)=0;
for i=1:npoints
    for j=1:npoints
        K(i,j)= k*(pi*j/l)^2*sin(pi*xn(i)*j/l);
    end
    f(i,1)=Qn(i);
end

a=inv(K)*f;

fi2b=0;
for i=1:npoints
    fi2b=fi2b+a(i)*sin(pi*x/l);
end

fit=[fit;fi2b];

figure(3)
mdibuja(x,fit,'x','fi','fi vs fi');


npoints=15;
for i=1:npoints
    xn(i)=x(round(i*length(x)/(npoints+1)));
    Qn(i)=Qt(round(i*length(x)/(npoints+1)));
end

for i=1:npoints
    for j=1:npoints
        K(i,j)= k*(pi*j/l)^2*sin(pi*xn(i)*j/l);
    end
    f(i,1)=Qn(i);
end

a=inv(K)*f;

fi2c=0;
for i=1:npoints
    fi2c=fi2c+a(i)*sin(pi*x/l);
end

fit=[fit;fi2c];

figure(3)
mdibuja(x,fit,'x','fi','fi vs fi');


% case 1: polynomials
% Nj(x)=[(xL/2)^2 - (x - xL/2)^2]^i
% Kij= k [(xL/2)^2 - (x-xL/2)^2]^(i-2) [-2 i [(xL/2)^2 - (x-xL/2)^2]- 2 (x-
% xL^2) i (i-1)]
% fi = Q(xi)

npoints=1;

xn(1)=x(round(length(x)/2));

Qn(1)=Q/2;

K(1,1)= -2*k;
f(1)=Qn(1);

a=inv(K)*f;

fi1a=a*((xL/2)^2 - (x-xL/2).^2);

fit=[fit;fi1a];

figure(2)
mdibuja(x,fit,'x','fi','fi vs fi');


npoints=2;
for i=1:npoints
    xn(i)=x(round(i*length(x)/(npoints+1)));
    Qn(i)=Qt(round(i*length(x)/(npoints+1)));
end

for i=1:npoints
    K(i,1)= -2*k;
    K(i,2)=k*(12*(xn(i)-(xL/2))^2-4*((xL/2)^2));
    cte=(xL/2)^2;
    y=xn(i)-(xL/2);
    for j=2:npoints
        K(i,j)=-2*j*((cte-y^2)^(j-2))*(cte-(2*j-1)*(y^2));
    end
    f(i,1)=Qn(i);
end

a=inv(K)*f;

fi1b=0;
for i=1:npoints
    y=(xL/2)^2 - (x - xL/2).^2;
    fi1b=fi1b+a(i)*(y.^i);
end

fit=[fit;fi1b];

figure(3)
mdibuja(x,fit,'x','fi','fi vs fi');


npoints=5;
for i=1:npoints
    xn(i)=x(round(i*length(x)/(npoints+1)));
    Qn(i)=Qt(round(i*length(x)/(npoints+1)));
end

for i=1:npoints
    K(i,1)= -2*k;
    K(i,2)=k*(12*(xn(i)-(xL/2))^2-4*((xL/2)^2));
    cte=(xL/2)^2;
    y=xn(i)-(xL/2);
    for j=2:npoints
        K(i,j)=-2*j*((cte-y^2)^(j-2))*(cte-(2*j-1)*(y^2));
    end
    f(i,1)=Qn(i);
end

a=inv(K)*f;

fi1c=0;
for i=1:npoints
    y=(xL/2)^2 - (x - xL/2).^2;
    fi1c=fi1c+a(i)*(y.^i);
end

fit=[fit;fi1c];

figure(3)
mdibuja(x,fit,'x','fi','fi vs fi');
 


% plot including labels, grid and ajusting axis...
function mdibuja(x,y,labx,laby,tit)

plot(x,y);
miny=min(1.1*min(min(y)),0.9*min(min(y)));
maxy=max(1.1*max(max(y)),0.9*max(max(y)));
axis([min(x) max(x) miny maxy]);
grid;
xlabel(labx);
ylabel(laby);
title(tit);



Download the Source code
Back to the example
Personal tools
Categories