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

% Resolution of the Poisson's equation with the WRM 
% using global Shape Functions
% Least Squares 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;

figure(1)
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 3: 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;
K(1,1)=k*(pi/2)*(pi/l)^2;
f(1)=Qn(1);

a=inv(K)*f;

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

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

npoints=2;
xn(1)=x0;
xn(2)=x(round(length(x)/2));
xn(3)=xL;

Qn(1)=Q;
Qn(2)=0;

for i=1:npoints
    for j=1:npoints
        if (i==j)
            K(i,j)=k*((pi*j/l)^2)*l/2;
        else
            K(i,j)= k*((pi*j/l)^2)*0.5*(l/pi)*(1/(i-j)*sin(pi*(i-j))-1/(i+j)*sin(pi*(i+j)));
        end
    end
    f(i,1)=Q*l*(1-cos(i*pi/2))/(i*pi);
end

a=inv(K)*f;

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

fit=[fit;fi2b];

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


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

for i=1:npoints
    for j=1:npoints
        if (i==j)
            K(i,j)=k*((pi*j/l)^2)*l/2;
        else
            K(i,j)= k*((pi*j/l)^2)*0.5*(l/pi)*(1/(i-j)*sin(pi*(i-j))-1/(i+j)*sin(pi*(i+j)));
        end
    end
    f(i,1)=Q*l*(1-cos(i*pi/2))/(i*pi);
end

a=inv(K)*f;

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

fit=[fit;fi2c];

figure(6)
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