[MATLAB数学相关] MATLBA有限差分程序分享

[aiyubaozi]

有问题请联系QQ： 291873404
function [U,err,zuidaERR,Ujingque]=tuoyuan1(jingqueu1,f1,phi1,phi2,phi3,phi4,a,b,c,d,m,n,epsilon)
%此程序用于求解差分方程-[D2u(x)+D2u(y)]=f(x,y)的解
%u(a,y)=phi1,u(b,y)=phi2,u(x,c)=phi3,u(x,d)=phi4
%差分剖分x方向m等分，差分剖分y方向n等分
h1=(b-a)/m;
h2=(d-c)/n;
h=2*[1/h1.^2+1/h2.^2];
x=a:h1:b;
y=c:h2:d;
%u为迭代初始矩阵
u=ones(m+1,n+1);U=zeros(m+1,n+1);
u(1,1:1:n+1)=feval(phi1,y); U(1,1:1:n+1)=feval(phi1,y);
u(1+m,1:1:n+1)=feval(phi2,y); U(1+m,1:1:n+1)=feval(phi2,y);
u(1:m+1,1)=[feval(phi3,x)]'; U(1:m+1,1)=[feval(phi3,x)]';
u(1:m+1,n+1)=[feval(phi4,x)]';U(1:m+1,n+1)=[feval(phi4,x)]';
fanshu=1;
while (fanshu>epsilon)
for i=2:m
j=2:n;
U(i,j)=[feval(f1,a+(i-1)*h1,c+h2:h2:d-h2)+(U(i,j-1)+u(i,j+1))/h2.^2+(U(i-1,j)+u(i+1,j))/h1.^2]/h;
end
fanshu=norm(u-U,inf);
u=U;
end
for i=1:m+1
j=1:n+1;
Ujingque(i,j)=feval(jingqueu1,a+(i-1)*h1,y);
end
err=abs(U-Ujingque);
zuidaERR=max(max(err));
zuidaERR=vpa(zuidaERR,4);
%U=[U(m/4+1,n/4+1),U(m/2+1,n/4+1),U(m*3/4+1,n/4+1),U(m/4+1,n/2+1),U(m/2+1,n/2+1)];
%U=vpa(U,7);
%ERR=[err(m/4+1,n/4+1),err(m/2+1,n/4+1),err(m*3/4+1,n/4+1),err(m/4+1,n/2+1),err(m/2+1,n/2+1)];
%ERR=vpa(ERR,4);