matlab计算拉普拉斯计算图解,急求用matlab编写解拉普拉斯方程的程序

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

这么处理狄拉克边界条件，我的精确解和数值解的误差很大部只哪错了

附程序

% 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;