用matlab解拉普拉斯法,急求用matlab编写解拉普拉斯方程的程序

今天药忘吃喽~ 2023-01-16 09:50 179阅读 0赞

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

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

附程序

% 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

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

%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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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;

用matlab解拉普拉斯法,急求用matlab编写解拉普拉斯方程的程序

发表评论取消回复

还没有评论，来说两句吧...

相关阅读

相关 Matlab：解方程的仿真

相关 matlab中.^的用法

相关用MATLAB求状态空间的解,实验二利用MATLAB求取线性系统的状态空间模型的解

相关用matlab解拉普拉斯法,急求用matlab编写解拉普拉斯方程的程序

相关 matlab复数方程的根,matlab解超越方程的复数根

相关用matlab进行拉普拉斯滤波,matlab拉普拉斯滤波器

相关 Matlab求微分方程的数值解

相关 Matlab求微分方程的解析解（dsolve）

相关 eval的用法（matlab）

相关 matlab中@的用法