有限差分法的Matlab程序（椭圆型方程）
function FD_PDE(fun,gun,a,b,c,d)
% 用有限差分法求解矩形域上的Poisson方程
tol=10^(-6); % 误差界
N=1000; % 最大迭代次数
n=20; % x轴方向的网格数
m=20; % y轴方向的网格数
h=(b-a)/n; % x轴方向的步长
l=(d-c)/m; % y轴方向的步长
for i=1:n-1
x(i)=a+i*h;
end % 定义网格点坐标
for j=1:m-1
y(j)=c+j*l;
end % 定义网格点坐标
u=zeros(n-1,m-1); %对u赋初值
% 下面定义几个参数
r=h^2/l^2;
s=2*(1+r);
k=1;
% 应用Gauss-Seidel法求解差分方程
while k<=N
% 对靠近上边界的网格点进行处理
% 对左上角的网格点进行处理
z=(-h^2*fun(x(1),y(m-1))+gun(a,y(m-1))+r*gun(x(1),d)+r*u(1,m-2)+u(2,m-1))/s;
norm=abs(z-u(1,m-1));
u(1,m-1)=z;
% 对靠近上边界的除第一点和最后点外网格点进行处理
for i=2:n-2
z=(-h^2*fun(x(i),y(m-1))+r*gun(x(i),d)+r*u(i,m-2)+u(i+1,m-1)+u(i-1,m-1))/s;
if abs(u(i,m-1)-z)>norm;
norm=abs(u(i,m-1)-z);
end
u(i,m-1)=z;
end
% 对右上角的网格点进行处理
z=(-h^2*fun(x(n-1),y(m-1))+gun(b,y(m-1))+r*gun(x(n-1),d)+r*u(n-1,m-2)+u(n-2,m-1))/s;
if abs(u(n-1,m-1)-z)>norm
norm=abs(u(n-1,m-1)-z);
end
u(n-1,m-1)=z;
% 对不靠近上下边界的网格点进行处理
for j=m-2:-1:2
% 对靠近左边界的网格点进行处理
z=(-h^2*fun(x(1),y(j))+gun(a,y(j))+r*u(1,j+1)+r*u(1,j-1)+u(2,j))/s;
if abs(u(1,j)-z)>norm
norm=abs(u(1,j)-z);
end
u(1,j)=z;
% 对不靠近左右边界的网格点进行处理
for i=2:n-2
z=(-h^2*fun(x(i),y(j))+u(i-1,j)+r*u(i,j+1)+r*u(i,j-1)+u(i+1,j))/s;
if abs(u(i,j)-z)>norm
norm=abs(u(i,j)-z);
end
u(i,j)=z;
end
% 对靠近右边界的网格点进行处理
matlab 下载 z=(-h^2*fun(x(n-1),y(j))+gun(b,y(j))+r*u(n-1,j+1)+r*u(n-1,j-1)+u(n-2,j))/s;
if abs(u(n-1,j)-z)>norm
norm=abs(u(n-1,j)-z);
end
u(n-1,j)=z;
end
% 对靠近下边界的网格点进行处理
% 对左下角的网格点进行处理
z=(-h^2*fun(x(1),y(1))+gun(a,y(1))+r*gun(x(1),c)+r*u(1,2)+u(2,1))/s;
if abs(u(1,1)-z)>norm
norm=abs(u(1,1)-z);
end
u(1,1)=z;
% 对靠近下边界的除第一点和最后点外网格点进行处理
for i=2:n-2
z=(-h^2*fun(x(i),y(1))+r*gun(x(i),c)+r*u(i,2)+u(i+1,1)+u(i-1,1))/s;
if abs(u(i,1)-z)>norm
norm=abs(u(i,1)-z);
end
u(i,1)=z;
end
% 对右下角的网格点进行处理
z=(-h^2*fun(x(n-1),y(1))+gun(b,y(1))+r*gun(x(n-1),c)+r*u(n-1,2)+u(n-2,1))/s;
if abs(u(n-1,1)-z)>norm
norm=abs(u(n-1,1)-z);
end
u(n-1,1)=z;
% 结果输出
if norm<=tol
fid = fopen('', 'wt');
fprintf(fid,'\n********用有限差分法求解矩形域上Poisson方程的输出结果********\n\n');
fprintf(fid,'迭代次数: %d次\n\n',k);
fprintf(fid,' x的值 y的值 u的值 u的真实值 |u-u(x,y)|\n');
for i=1:n-1
for j=1:m-1
fprintf(fid, '%8.3f %8.3f %14.8f %14.8f %14.8f\n', [x(i),y(j),u(i,j),gun(x(i),y(j)),abs(u(i,j)-gun(x(i),y(j)))]);
end
end
fclose(fid);
break; % 用来结束while循环
end
k=k+1;
end
if k==N+1
fid = fopen('', 'wt');
fprintf(fid,'超过最大迭代次数，求解失败！');
fclose(fid);
end
clc
[a1 a2 a3 a4] = textread('F:\aa.txt','%f %f %f %f');
a = [a1 a2 a3];
a=a';
b=a4';
[pa,mina,maxa,pb,minb,maxb]=premnmx(a,b);
net =newrb(pa,pb,0,1.3,24,2);
an =sim(net,pa);
E = an - pb;
m =sse(E)
n = mse(E)
[f1 f2 f3 f4]= textread('F:\bb.txt','%f %f %f %f');
f = [f1 f2 f3];
f=f';
pf = tramnmx(f,mina,maxa);
an2 = sim(net,pf);
g =postmnmx(an2,minb,maxb);
g= g';
E2 = g- f4;
mm =sse(E2)
nn = mse(E2)
欢迎您的下载，
资料仅供参考！
致力为企业和个人提供合同协议，策划案计划书，学习资料等等
打造全网一站式需求
发表评论