马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
% 求解在区域R={(x,y):0≤x≤a,0≤y≤b}内的uxx(x,y)+uyy(x,y)=0的近似解,而且满足条件u(x,0)=f1(x), u(x,b)=f2(x),
%其中0≤x≤a且u(0,y)=f3(y), u(a,y)=f4(y),其中0≤y≤b.设Δx=Δy=h,而且存在整数n和m,使得a=nh,b=mh.
%function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1)
%Input -f1,f2,f3,f4 are boundary functions input as strings
% -a and b right end points of [0,a] and [0,b]
% -h step size
% -tol is the tolerance
%Output -U solution matrix
%Initialize parameters and U
function U=dirich(20,180,80,0,4,4,0.5,0.001,10)
n=fix(a/h)+1;
m=fix(b/h)+1;
ave=(a*(feval(f1,0)+feval(f2,0))...
+b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b); %feval('fun',x)等价于求fun(x)值。
U=ave*ones(n,m);
%Boundary conditions
U(1,1:m)=feval(f3,0:h:(m-1)*h)';
U(n,1:m)=feval(f4,0:h:(m-1)*h)';
U(1:n,1)=feval(f1,0:h:(n-1)*h)';
U(1:n,m)=feval(f3,0:h:(n-1)*h)';
U(1,1)=(U(1,2)+U(2,1))/2;
U(1,m)=(U(1,m-1)+U(2,m))/2;
U(n,1)=(U(n-1,2)+U(n,2))/2;
U(n,m)=(U(n-1,m)+U(n,m-1))/2;
%SOR parameter
w=4/(2+sqrt(4-(cos(pi(n-1))+cos(pi/(m-1)))^2));
%Refine approximations and sweep operator throughout
%the grid
err=1;
cnt=0;
while((err>tol)&(cnt<=max1))
err=0;
for j=2:m-1
for i=2:n-1
relx=w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4;
U(i,j)=U(i,j)+relx;
if(err<=abs(relx))
err=abs(relx);
end
end
end
cnt=cnt+1;
end
U=flipud(U'); %flipud(X) returns X with columns preserved and rows flipped in the up/down direction.
运行后出现
??? Error: File: D:\doctor\programe\self-consistent\dirich.m Line: 13 Column: 19
"identifier" expected, "numeric value" found.
[ 本帖最后由 ChaChing 于 2009-2-4 10:17 编辑 ] |