该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
这么处理狄拉克边界条件,我的精确解和数值解的误差很大部只哪错了
附程序
% examp1a2.m
clear all
clc
N=5;M=5;n=N*M;
L=2;h=2;
tic
c=sqrt((L/(N-1))^2+(h/(M-1))^2);
dm=2.5*c;
mu=0.25;
E=1e5;
beta=1*10^5;
for j=1:M
for i=1:N
x(i+(j-1)*N)=L/(N-1)*(i-1);
y(i+(j-1)*N)=h/(M-1)*(j-(M+1)/2);
end
end
t=[-0.8611363;-0.339880;0.339880;0.8611363];
G=[0.3478548;0.6521452;0.6521452;0.3478548];
for j=1:(M-1)
for i=1:(N-1)
for ky=1:4
for kx=1:4
xg(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=(x(i+1+(j-1)*N)-x(i+(j-1)*N))/2*t(kx)+(x(i+1+(j-1)*N)+x(i+(j-1)*N))/2;
yg(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=(y(i+j*N)-y(i+(j-1)*N))/2*t(ky)+(y(i+j*N)+y(i+(j-1)*N))/2;
Gx(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=G(kx)*(x(i+1+(j-1)*N)-x(i+(j-1)*N))/2;
Gy(16*(i-1)+16*(j-1)*(N-1)+4*(ky-1)+kx)=G(ky)*(y(i+j*N)-y(i+(j-1)*N))/2;
end
end
end
end
%inputing P
for i=1:n;
P(i,1)=1;
P(i,2)=x(i);
P(i,3)=y(i);
P(i,4)=(x(i))^2;
P(i,5)=x(i)*y(i);
P(i,6)=(y(i))^2;
end
%calculating K
K=0;
for i=1:(4*4*(N-1)*(M-1))
for j=1:n
d(j)=sqrt((xg(i)-x(j))^2+(yg(i)-y(j))^2);
if d(j)<=dm
W(j,j)=(exp(-(d(j)/c)^2)-exp(-6.25))/(1-exp(-6.25));
Wx(j,j)=-exp(-(d(j)/c)^2)*2*(xg(i)-x(j))/c/c/(1-exp(-6.25));
Wy(j,j)=-exp(-(d(j)/c)^2)*2*(yg(i)-y(j))/c/c/(1-exp(-6.25));
else
W(j,j)=0;
Wx(j,j)=0;
Wy(j,j)=0;
end
end
p=[1;xg(i);yg(i);(xg(i))^2;xg(i)*yg(i);(yg(i))^2];
px=[0,1,0,2*xg(i),yg(i),0];
py=[0,0,1,0,xg(i),2*yg(i)];
A=P'*W*P;
Ax=P'*Wx*P;
Ay=P'*Wy*P;
B=P'*W;
Bx=P'*Wx;
By=P'*Wy;
ph=p'*inv(A)*B;
phx=px*inv(A)*B-p'*inv(A)*Ax*inv(A)*B+p'*inv(A)*Bx;
phy=py*inv(A)*B-p'*inv(A)*Ay*inv(A)*B+p'*inv(A)*By;
for j=1:n
phi(1,2*j-1)=ph(j);
phi(1,2*j)=0;
phi(2,2*j-1)=0;
phi(2,2*j)=ph(j);
Ca(1,2*j-1)=phx(j);
Ca(1,2*j)=0;
Ca(2,2*j-1)=0;
Ca(2,2*j)=phx(j);
Cb(1,2*j-1)=phy(j);
Cb(1,2*j)=0;
Cb(2,2*j-1)=0;
Cb(2,2*j)=phy(j);
end
K=K+Gx(i)*Gy(i)*(Ca'*Ca+Cb'*Cb);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j=1:(M-1)
for ik=1:4
xgb(4*(j-1)+ik)=L;
ygb(4*(j-1)+ik)=(y((j+1)*N)-y(j*N))/2*t(ik)+(y((j+1)*N)+y(j*N))/2;
Gyb(4*(j-1)+ik)=G(ik)*(y((j+1)*N)-y(j*N))/2;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:(4*(M-1))
for j=1:n
db(j)=sqrt((xgb(i)-x(j))^2+(ygb(i)-y(j))^2);
if db(j)<=dm
Wb(j,j)=(exp(-(db(j)/c)^2)-exp(-6.25))/(1-exp(-6.25));
else
Wb(j,j)=0;
end
end
pb=[1;xgb(i);ygb(i);(xgb(i))^2;xgb(i)*ygb(i);(ygb(i))^2];
Ab=P'*Wb*P;
Bb=P'*Wb;
phb=pb'*inv(Ab)*Bb;
for j=1:n
phib(1,2*j-1)=phb(j);
phib(1,2*j)=0;
phib(2,2*j-1)=0;
phib(2,2*j)=phb(j);
end
K=K+beta*Gyb(i)*phib'*phib;