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
1 view (last 30 days)
Show older comments
Omar Idriss
on 17 Apr 2018
Commented: Rena Berman
on 21 Mar 2022
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
7 Comments
Accepted Answer
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.
0 Comments
More Answers (1)
Omar Idriss
on 17 Apr 2018
1 Comment
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.
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!