suanli3_5_1_modified.m

《偏微分方程数值解法》孙志忠版中算例的一些MATLAB程序

%向后Euler差分格式
%Ut-Uxx=0,0<x<1,0<t<1
%u(x,0)=exp(x)
%u(0,t)=exp(t),u(1,t)=exp(1+t)
%%%%%%%网格剖分%%%%%%%%%
% first set up the attributes
h = 0.1;  				% the step-size in the x dimension%空间步长
r=1;                  %步长比
dt = h^2 * r;			% the step-size in the t dimension%%时间步长
T0 = 1;	  			% the time-interval over which we will calculate u 时间长度
TN = ceil(T0 / dt);	% the number of steps in the time-interval%时间网格
X0 = 1;					% the x-interval over which we will calculate u空间长度
XN = ceil(X0 / h);	% the number of steps in the x-interval%空间网格
U = zeros(XN + 1, TN + 1);	% set up an array of zeros %初始化
% Next we set up the initial conditions%初始值及边界值
x = 0 :h: X0;
for j=1:XN+1
    U(j , 1) = exp( (j-1)*h );
end
for k=1:TN+1
    U(1,k) =exp( (k-1) *dt);U(XN+1,k)=exp(1+ ( k-1 ) *dt);
end
%%%%%%形成系数矩阵%%%%%%%%
A=diag((5/6+r)*ones(XN-1,1))+diag( (1/12-r/2) *ones(XN-2,1),1 )+diag( (1/12-r/2)*ones(XN-2,1),-1 );
B=diag((5/6-r)*ones(XN-1,1))+diag( (1/12+r/2) *ones(XN-2,1),1 )+diag( (1/12+r/2)*ones(XN-2,1),-1 );
%右端项初始化
F = zeros( XN-1,1 );
%迭代求解
for k = 1 : TN
    F(1)=(1/12+r/2)*exp((k-1)*dt)-(1/12-r/2)*exp(k*dt);
    F(XN-1)=(1/12+r/2)*exp((k-1)*dt+1)-(1/12-r/2)*exp(k*dt+1);
    U(2:XN,k+1) = A \ ( B * U(2:XN,k) + F );
end;
%给出精确解
ana = zeros(XN + 1, TN + 1);
for k = 1 : TN+1
	for j = 1 : XN+1
         ana(j, k ) = exp ( (j-1) *h + (k-1)*dt);
   end;
end;
%给出在T=1时的最大误差
max( abs ( U (:,TN+1) -ana (:,TN+1) ) )
%图示化
t=0:dt:T0;
[xx,tt]= meshgrid(x,t);
xx=xx';   tt=tt';
mesh(xx,tt,ana);
title('analytical solution');
ylabel(sprintf('steps in x domain\nfrom 0 to 1'));
xlabel(sprintf('steps in time-domain\nfrom 0 to 1'));
zlabel('U(x,y)');
pause
mesh(xx,tt,U);
title('numerical solution');
ylabel(sprintf('steps in x domain\nfrom 0 to 1'));
xlabel(sprintf('steps in time-domain\nfrom 0 to 1'));
zlabel('U(xi,tk)');
pause
mesh(xx,tt,U-ana);
title('diffeence beween analical and numeical soluion');
ylabel(sprintf('steps in x domain\nfrom 0 to 1'));
xlabel(sprintf('steps in time-domain\nfrom 0 to 1'));
zlabel('U(xi,tk)');
pause
plot(x,U(:,TN+1),'-or',x,ana(:,TN+1))
title('向前Euler格式数值解与精确解的比较')
legend('uh(x,1)','u(x,1)')%Crank-Nicolson schemes