clear
Re=100; %式(8)
L=5; %计算域，可以随意设置
u0=1; %内边界条件，u0可随意
miu=u0*L/Re; %扩散系数
N=41; %网格节点数
t=2; %时间步长(松弛因子)
h=L/(N-1); %网格步长
for i=1:N %网格划分
x(i)=(i-1)*h;
end
%---------------------G-S型迎风半隐格式计算-----------------------
u=zeros(N,1); %迭代初值
u(1)=u0; %内边界条件
u1=u;
tol=1;p=0;
while tol>=1e-5 %收敛精度，0.00001
for i=2:N-1
R(i)=abs(u(i))*h/2/miu;
u1(i)=(-t*h*(u(i)*(u(i+1)-u1(i-1)))+2*t*miu*(1+R(i))*(u(i+1)+u1(i-1))+2*h^2*u(i))/(2*h^2+4*t*miu*(1+R(i)));
end
tol=max(abs(u-u1)); %计算两个迭代层的误差
u=u1;
p=p+1;
if p>=5000
disp('G-S型迎风半隐格式不收敛！')
break
end
end
%------------------------G-S型Samarskii型半隐格式计算---------------
tol=1;p=0;
v=zeros(N,1); %迭代初值
v(1)=u0; %内边界条件
v1=v;
while tol>=1e-5 %收敛精度，0.00001
for i=2:N-1
R(i)=abs(v(i))*h/2/miu;
v1(i)=(-t*h*(v(i)*(v(i+1)-v1(i-1)))+2*t*miu*(1/(1+R(i))+R(i))*(v(i+1)+v1(i-1))+2*h^2*v(i))/(2*h^2+4*t*miu*(1/(1+R(i))+R(i)));
end
tol=max(abs(v1-v)); %计算量迭代层误差
v=v1;
p=p+1;
if p>=5000
disp('G-S型Samarskii型半隐格式不收敛！')
break
end
end
%-----------------------牛顿迭代法计算 --------------------
tol=1;U=2; %误差；迭代初值
while tol>=1e-10 %收敛精度0.0000000001
f=(U-1)/(U+1)-exp(-U*Re); %式(9)建立的函数
fd=1/(U+1)-(U-1)/(U+1)^2+Re/exp(U*Re); %式(9)建立的函数的导数
U1=U-f/fd;
tol=abs(U-U1); %计算误差
U=U1;
end
%--------------------------计算精确解-----------------------------------
for i=1:N
ue(i)=u0*U*((1-exp(U*Re*(x(i)/L-1)))/(1+exp(U*Re*(x(i)/L-1))));
end
plot(x,ue,'-*',x,u,'-or',x,v,'-<g') %作图
评论5