clear
clc
tic
X=linspace(-3.0,3.0,51);
Y=linspace(-3.0,3.0,51);
R=0.2;%支持域的尺寸
F=zeros(51,51);
for k1=1:51*51
i=mod((k1-1),51)+1 ; %行号
j=floor((k1-1)/51)+1; %列号
x=X(i);
y=Y(j);
p=[1;x;y];%基函数
%确定网格点(x,y)的影响区域的大小,确定影响域内的节点
In_Count=0;%支持域内点的个数
Subscript=zeros(50,2);%用来存储支持域内结点的下标
for k=1:51
for l=1:51
Xk=X(k);
Yl=Y(l);
if sqrt((x-Xk)^2+(y-Yl)^2)<=R;
In_Count=In_Count+1;
Subscript(In_Count,1)=k;
Subscript(In_Count,2)=l;
end
end
end
%计算点x的形函数
A=zeros(3,3);
B=zeros(3,In_Count);
Us=zeros(In_Count,1);
for m=1:In_Count
n1=Subscript(m,1);
n2=Subscript(m,2);
xn1=X(n1);
yn2=Y(n2);
pp=[1;xn1;yn2];
C=zeros(1,In_Count);
D=zeros(In_Count,1);
%权函数的选择
ri=sqrt((x-xn1)^2+(y-yn2)^2)/R;
if ri<=1
W=1-6*ri^2+8*ri^3-3*ri^4;
elseif ri>1
W=0;
end
A=W*pp*pp'+A;
C(1,m)=1;
B=W*pp*C+B;
D(m,1)=1;
z=call_xy(xn1,yn2);%影响域内每一个节点对应一个值
Us=z*D+Us;
end
O=p'*A^(-1)*B;
%计算网格点x处的节点值
U=O*Us;
F(i,j)=U;
end
figure
End=surf(X,Y,F);
xlabel('x'),ylabel('y'),zlabel('f');
t2=toc;
display(strcat('parfor串行计算时间:',num2str(t2),'秒'));