MATLAB Answers

hello everybody i am a student master second year i have a problem i have no knowledge about Matlab anyone can help me explain this code and thanks in advance

28 views (last 30 days)
Omar Idriss
Omar Idriss on 17 Apr 2018
Edited: Rena Berman on 6 May 2021 at 17:26
function kansa1_MQ_1D%
clear ALL;
% Local MPS method; Bi-Harmonic equation
% ni-----# of interior points
% nb-----# of boundary points
% c------shape parameter
% scale with maximum radius
%generate evenly distribute grid points on a unit square
tic
ns=3; nn=30; c=50;
n=nn+1 ; nb=2; ni=n-nb; bb=1;
x1=0:bb/nn:bb;
intnode=zeros(ni,1);
j=0;k=0; %generate the interior and boundary nodes evenly
for i=1:n
if(x1(i)==0 || x1(i)==bb)
j=j+1;
bdpt(j,1)=x1(i);
else
k=k+1;
intnode(k,1)=x1(i);
end
end
x=[intnode(:,1); bdpt(:,1)]';
B=localmpsmatrix(x,ni,nb,n,ns,c);%produce sparse matrix
spy(B)
% Exemple 1
% f(1:ni)=exp(x(1:ni)); % for excat solution exp(x)
% f(ni+1:n)=exp(x(ni+1:n));% for excat solution exp(x)
%Exemple 2
% f(1:ni)=- sin(x(1:ni)); % for excat solution sin(x)
% f(ni+1:n)=sin(x(ni+1:n));% for excat solution sin(x)
%
%Exemple 3
% f(1:ni)=- sin(x(1:ni))+exp(x(1:ni)); % for excat solution sin(x)
% f(ni+1:n)=sin(x(ni+1:n))+exp(x(ni+1:n));% for excat solution sin(x)
%Exemple 4
%f(1:ni)=3.*x(1:ni).^2-x(1:ni); % for excat solution x(x-1)
%f(ni+1:n)=0;% for excat solution x(x-1)
f(1:ni)=(x(1:ni).^2/2+x(1:ni)).*exp(x(1:ni))-x(1:ni).^2/2.*sin(x(1:ni))+x(1:ni).*cos(x(1:ni)); % for excat solution x(x-1)
f(ni+1:n)=sin(x(ni+1:n))+exp(x(ni+1:n));% for excat solution x(x-1)
alpha=B\f';
%u=exp(x)'; %exact solution
%u=sin(x)'; %exact solution
%u=(sin(x)+exp(x))'; %exact solution
%u=(x.*(x-1))';
u=(exp(x(1:n))+sin(x(1:n)))';
error1=max(abs((alpha(1:n)-u(1:n)))); %absolute maximu error
rms1=norm((alpha(1:n)-u(1:n)),2)/sqrt(n); % RMS error
toc
fprintf('n = %6d, nb = %4d, ns= %2d, c =%4d, \n',n,nb,ns,c);
fprintf('max error1: %e\nRMS1 error: %e\n', error1,rms1);
%hold on
%semilogy(param,error1,'b'),xlabel 'shape parameter' , ylabel 'error'
return
%%Produce the local sparse matrix
function localmpsmatrix=localmpsmatrix(x,ni,nb,n,ns,c)
t=x'; i1=0;
tree=kdtree_build(t); %search
[id, D]=kdtree_k_nearest_neighbors(tree,t(1,:),ns);
ctrs=x(id(1:ns))'; %center of locale domain
DM= DistanceMatrix(ctrs, ctrs); %Distance matrix
cc=D(1)*c;
c2=cc*cc;
% r2 = D.*D;
% mq=sqrt(DM+c2);
% A=mq;
% V1=sqrt(r2);
% V=c2./(r2 + c2).^(3/2);
%
% coef=A\V;
B1=zeros(1,ns*ni+nb); B2=B1; B3=B1;
for i=1:ni
% [id]=kdtree_k_nearest_neighbors(tree,t(i,:),ns);
[id, D]=kdtree_k_nearest_neighbors(tree,t(i,:),ns);
ctrs=x(id(1:ns))'; %center of locale domain
DM= DistanceMatrix(ctrs, ctrs); %Distance matrix
% cc=D(1)*c;
% c2=cc*cc;
r2 = D.*D;
mq=sqrt(DM+c2);
A=mq;
%%for solving -div(K*grad(u))=f
V1=(t(i,:)-ctrs(:))./(r2 + c2).^(1/2);
V2=c2./(r2 + c2).^(3/2);
V=t(i,:).^2*V2/2.+t(i,:).*V1;
coef=A\V;
for k=1:ns %sparse matrix storage for interior points
i1=i1+1;
B1(i1)=i; B2(i1)=id(k); B3(i1)=coef(k);
end
end
for i=ni+1:n %indentity matrix for boundary points
i1=i1+1;
B1(i1)=i; B2(i1)=i; B3(i1)=1;
end
localmpsmatrix = sparse(B1(1:i1),B2(1:i1),B3(1:i1),n,n);
return
%%DistanceMatrix: r^2
function DM = DistanceMatrix(dsites, ctrs)
[M,s] = size(dsites); N = length(ctrs);
DM = zeros(M,N);
for i=1:s
[dr,cc]=ndgrid(dsites(:,i),ctrs(:,i));
DM = DM +(dr-cc).^2;
end
return
  6 Comments

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 17 Apr 2018
bdpt is a temporary variable used to construct a list of border points. intnode is a temporary variable used to construct a list of interior nodes. The complete list of nodes is constructed as the interior nodes followed by the border nodes.
In a 1D situation, there are exactly two border nodes, one at the beginning and one at the end. The other n-2 nodes are interior nodes.

More Answers (1)

Omar Idriss
Omar Idriss on 17 Apr 2018
Ok thank you very much Walter Roberson, I have another question please I did not understand this writing B\f'; B a matrix and f a vector ? -div(K*grad(u))=f, I did not understand how we resolved this,,,,..??? CORDIALLY
  1 Comment
Walter Roberson
Walter Roberson on 17 Apr 2018
Sorry, I am not familiar with the computation being done.
B\f is like inv(B)*f for the case where B is invertable and f is a column vector -- that is, one of its main uses is to solve simultaneous equations. It is, though, not restricted to that situation: it can do a least-squared fit for over-determined equations as well.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by