# I want to run this code?

3 views (last 30 days)
Akhtar Jan on 9 Mar 2022
Commented: Akhtar Jan on 9 Mar 2022
shape = 2.5;
N = 30;
Np = N^2;
[X,Y] = meshgrid(linespace(0,1,N),linespace(0,1,N));
x = reshape(X,N^2,1);
y = reshape(Y,N^2,1);
f = -2*(2*y.^3-3*y.^2+1)+6*(1-x.^2).*(2*y-1);
u = (1-x.^2).*(2*y.^3-3*y.^2+1); % exact solution
H = zeros(Np,Np);
rx = zeros(Np,Np);
ry = zeros(Np,Np);
r = zeros(Np,Np);
dirichletBCs = find(x==0 | x==1); % identify boundry and interior centers
neumannBCs = find ((y==0 | y==1)& ~(x==0 | x==1));
interior = find('x~=0' & 'x~ = 1' & 'y~ = 0' & 'y~ = 1' );
f (dirichletBCs)=u(dirichletBCs);
f (neumannBCs)= 0;
for i =1:Np
for j =1:Np
rx ( i , j ) = x ( i )-x ( j ) ;
ry ( i , j ) = y ( i )-y ( j ) ;
r ( i , j ) = sqrt ( rx ( i , j )^2 + ry ( i , j )^2 );
end
end
H(dirichletBCs,:)=mq(r(dirichletBCs,:),shape); % e v a l u a t i on matrix
H(neumannBCs,:)= mqDerivatives(r(neumannBCs,:),ry(neumannBCs,:),shape,1);
Hxx = mqDerivatives(r(interior,:),rx(interior,:),shape,2);
Hyy = mqDerivatives(r(interior,:),ry(interior,:),shape,2);
H(interior,:)=Hxx+Hyy;
kappa=cond(H);
alpha=H\f;
B = mq( r , shape ) ; % system matrix
uh = B*alpha;
error = norm(uh-u,inf);
t = delaunay(x,y);trisurf(t,x,y,abs(uh-u));
xlabel 'x' , ylabel 'y' colormap('Summer')
Walter Roberson on 9 Mar 2022
What problems are you having?

DGM on 9 Mar 2022
I can't test anything, since there are missing functions, but:
% there is no linespace()
%[X,Y] = meshgrid(linespace(0,1,N),linespace(0,1,N));
[X,Y] = meshgrid(linspace(0,1,N),linspace(0,1,N));
% this is an inexplicable misuse of char vectors and invalid operator spacing
%interior = find('x~=0' & 'x~ = 1' & 'y~ = 0' & 'y~ = 1' );
interior = find(x~=0 & x~=1 & y~=0 & y~=1 );
% using command syntax for xlabel,ylabel is technically valid
% but the line continuation will cause an error at colormap()
%xlabel 'x' , ylabel 'y' colormap('Summer')
xlabel('x')
ylabel('y')
colormap('Summer')
Akhtar Jan on 9 Mar 2022
let me check it....

R2015a

### Community Treasure Hunt

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

Start Hunting!

Translated by